企业官网建设流程全解析

1. 从信号到频谱：为什么我们需要FFT？

如果你正在处理声音、振动、通信信号或者任何随时间变化的物理量，那么你迟早会碰到一个核心问题：这个信号里到底包含了哪些频率成分？它的“音色”是怎样的？是单一频率的纯音，还是多种频率混杂的噪音？要回答这些问题，时域里那条上下波动的曲线就显得力不从心了。这时，我们就需要一把“频率的尺子”，把信号从时间的世界转换到频率的世界去观察。这把尺子，就是傅里叶变换。

傅里叶变换的理论很美，但计算起来很麻烦，尤其是对计算机处理的离散数字信号。直到快速傅里叶变换算法的出现，它才真正成为工程师和科学家手中的利器。FFT不是一种新的变换，而是计算离散傅里叶变换的一种高效算法，它能把计算复杂度从 O(N²) 降到 O(N log N)。想象一下，你要分析一段1秒钟、采样率44100Hz的音频，数据点有44100个。用原始的DFT方法，计算量是天文数字；而用FFT，可能眨眼间就完成了。这就是为什么FFT无处不在，从手机里的音乐播放器到雷达的信号处理芯片，背后都有它的身影。

在MATLAB这个工程计算的神器中，FFT功能被封装得极其友好，几乎是一行代码就能完成从时域到频域的华丽转身。但“友好”不等于“简单”，更不等于“用了就懂”。很多新手，甚至是有一定经验的使用者，常常在几个关键环节上栽跟头：频谱的横坐标到底对应什么物理频率？幅度谱为什么要除以N或者乘以2？单边谱和双边谱有什么区别？相位信息怎么提取才靠谱？这些看似基础的问题，恰恰是区分“会敲代码”和“真懂原理”的关键。

这篇文章，我就结合自己多年在信号处理项目中的实战经验，带你彻底搞懂MATLAB中FFT的使用。我们不只讲fft(x)这一行命令，更要深挖它背后的每一个参数、输出的每一个数组元素的物理意义，以及如何避免那些教科书里不提、但实际工作中一定会遇到的坑。无论你是正在做课程设计的学生，还是需要快速验证算法的工程师，相信这些从实际项目中沉淀下来的细节，都能让你少走弯路。

2. 核心概念前置：采样、点数与频率分辨率

在动手写代码之前，我们必须把几个基础概念夯实地基。这些概念直接决定了你FFT结果的正确性。

2.1 采样定理与奈奎斯特频率

我们的计算机无法处理连续的模拟信号，必须每隔一段时间（采样间隔Ts）对信号进行一次“拍照”，得到一系列离散的数据点，这个过程就是采样。采样率 Fs 就是每秒拍照的次数，单位是Hz。根据奈奎斯特-香农采样定理，为了能够从采样后的数字信号中无失真地恢复原始模拟信号，采样频率 Fs 必须至少是信号中最高频率成分（Fmax）的两倍，即 Fs > 2 * Fmax。

这里引出一个关键概念：奈奎斯特频率（Fs/2）。它是FFT所能分析的最高频率。任何高于 Fs/2 的频率成分，在采样后都会“伪装”成低于 Fs/2 的频率，这种现象称为“混叠”。所以，在采样前，通常需要通过一个抗混叠滤波器，把信号中高于 Fs/2 的频率成分滤除掉。

假设我们以 Fs = 1000 Hz 的速率对信号采样了1秒钟，那么我们就得到了 N = 1000 个数据点。在MATLAB中，这通常是一个长度为1000的行向量或列向量。

2.2 频率分辨率：你能看清多近的两根“谱线”？

频率分辨率 Δf 指的是FFT结果中相邻两个频率点之间的间隔。它决定了你能区分开两个频率多么接近的正弦波。计算公式非常简单：

Δf = Fs / N

其中，N是参与FFT运算的数据点数。注意，这个N不一定等于你采集到的总数据长度，后面我们会讲到“补零”操作。

从这个公式你可以直观地理解：采样率 Fs 固定时，你分析的数据越长（N越大），频率分辨率就越高（Δf越小），你看频谱就越“精细”。反之，数据越短，频谱就越“粗糙”。例如，Fs=1000Hz，如果取N=1000个点做FFT，那么 Δf = 1 Hz。这意味着频谱图上每间隔1Hz有一个点。如果两个正弦波的频率相差小于1Hz，它们的频谱峰可能会混在一起，无法分辨。

2.3 FFT点数N的选择：2的幂次方与补零

FFT算法对数据长度N有偏好。当N是2的整数次幂（如256， 512， 1024， 2048）时，算法的计算效率最高。MATLAB的fft函数对任意长度的N都能计算，但内部可能会采用不同的优化策略。

如果你的数据长度不是2的幂次方，通常有两种做法：

1. 直接计算：Y = fft(x)，其中length(x)不是2的幂。MATLAB会处理，但速度可能稍慢。
2. 补零：Y = fft(x, NFFT)，其中NFFT是一个大于length(x)的2的幂次数。例如，x有600个点，你可以设置NFFT = 1024。

补零操作需要深刻理解：它并不能提高真实的频率分辨率！因为补零并没有增加原始信号的实际信息。它的主要作用是：

• 使频谱图看起来更平滑：在原有的频率点之间插值，让曲线更连续美观。
• 便于取2的幂次方，提高计算效率。
• 可能使频率峰值的位置看起来更精确，尤其是当信号频率不是Δf的整数倍时。

真正的频率分辨率只由原始数据长度和采样率决定：Δf_true = Fs / N_original。补零后的分辨率 Δf_apparent = Fs / NFFT，这只是“视觉分辨率”，而非“物理分辨率”。

3. MATLAB FFT实战：从向量到有物理意义的频谱

现在，我们用一个完整的例子，把理论变成代码。我们的目标是：生成一个包含多个频率成分的合成信号，然后用FFT分析它，并得到一张横坐标是物理频率（Hz）、纵坐标是真实幅度（与原始信号一致）的频谱图。

3.1 构造一个测试信号

我们构造一个包含50Hz（幅度1.5）、120Hz（幅度1）和两个高频噪声（200Hz， 幅度0.3；310Hz，幅度0.8）的信号。采样率设为1000Hz，采样时长0.5秒。

%% 1. 参数设置与信号生成 Fs = 1000; % 采样率 (Hz) T = 0.5; % 信号时长 (秒) t = 0:1/Fs:T-1/Fs; % 时间向量，注意‘-1/Fs’以确保点数为 Fs*T N = length(t); % 信号点数 % 生成信号成分 comp1 = 1.5 * sin(2*pi*50*t); % 50 Hz comp2 = 1.0 * sin(2*pi*120*t); % 120 Hz comp3 = 0.3 * sin(2*pi*200*t); % 200 Hz comp4 = 0.8 * sin(2*pi*310*t); % 310 Hz % 合成信号，并加入一些随机噪声 x = comp1 + comp2 + comp3 + comp4 + 0.1*randn(size(t)); % 绘制时域信号 figure(‘Position‘， [100, 100, 800, 400]) subplot(2,1,1) plot(t, x) xlabel(‘时间 (秒)‘) ylabel(‘幅度‘) title(‘时域信号 (含噪声)‘) grid on

运行这部分代码，你会看到时域信号是一条复杂的、看似无规律的波形，无法直接看出里面含有50Hz和120Hz的成分。

3.2 执行FFT与计算双边频谱

接下来，我们对信号x做FFT。这里我们直接使用数据原始长度N。

%% 2. 执行FFT X = fft(x); % X是一个复数数组，包含频域信息

X是一个和x长度相同的复数数组。它的第一个元素X(1)对应的是直流分量（0Hz）。在MATLAB中，FFT输出的频率排列顺序是：从0Hz到正频率，再到负频率。 具体来说：X(1)是0Hz，X(2)到X(N/2+1)对应从Δf到Fs/2的正频率，X(N/2+2)到X(N)对应从-Fs/2+Δf到-Δf的负频率（对于实数信号，这部分是正频率部分的共轭对称）。

为了得到每个频率点对应的物理频率值，我们需要构建频率向量。

%% 3. 构建频率向量 (双边谱) f = (-N/2 : N/2-1) * (Fs/N); % 以0Hz为中心，从 -Fs/2 到 Fs/2-Δf % 注意：这种构建方式要求N是偶数。如果N是奇数，公式需微调。 % 更通用的方法是使用fftshift配合构建单边频率 X_shifted = fftshift(X); % 将零频分量移动到频谱中心 f_shifted = (-N/2 : N/2-1) * (Fs/N); % 与X_shifted对应的频率向量 % 计算双边幅度谱 magnitude_double = abs(X) / N; % 注意这里除以了N magnitude_shifted = abs(X_shifted) / N;

关键点1：为什么要除以N？DFT的定义中包含了求和。对于一个纯正弦波A*sin(2πf t)，其FFT结果在对应频率点上的谱线幅度（在忽略频谱泄漏的理想情况下）大约是A * N / 2。为了从FFT结果X中恢复出原始信号的真实幅度A，我们需要abs(X) * 2 / N（对于非直流分量）。而先abs(X)/N得到的是“单边谱幅度”的一半（即A/2），这在构建单边谱时会很清晰。这里先除以N是一个中间步骤。

3.3 转换为更实用的单边幅度谱

对于实数信号，其频谱是共轭对称的，负频率部分不提供新的信息。因此，我们通常只显示从0Hz到奈奎斯特频率Fs/2的部分，这就是单边谱。

%% 4. 计算单边幅度谱 P2 = abs(X)/N; % 双边谱幅度 (除以N后) P1 = P2(1:N/2+1); % 取前半部分 (0Hz 到 Fs/2) P1(2:end-1) = 2*P1(2:end-1); % 除了直流分量(0Hz)，其他频率分量幅度乘2 % 构建单边谱频率向量 f_single = (0:(N/2)) * (Fs/N); % 从0Hz到Fs/2 % 绘制单边幅度谱 subplot(2,1,2) stem(f_single, P1, ‘LineWidth‘， 1.5) xlabel(‘频率 (Hz)‘) ylabel(‘幅度 |P1(f)|‘) title(‘单边幅度谱‘) xlim([0, Fs/2]) % 通常只显示到Fs/2 grid on

关键点2：为什么单边谱的非直流分量要乘以2？因为FFT计算出的双边谱中，信号的总能量被平均分配在了正负频率两个峰上，每个峰的幅度是真实幅度A的一半。当我们只显示正频率部分时，需要将幅度乘以2，才能代表该频率成分的真实幅度。直流分量（0Hz）的能量只存在于一个点上，所以不需要乘2。

现在观察频谱图，你应该能在50Hz、120Hz、200Hz和310Hz附近看到清晰的谱峰，并且它们的幅度分别接近1.5， 1.0， 0.3和0.8。噪声则会表现为整个频带上的低矮“基底”。

3.4 相位谱的提取与解读

幅度谱告诉我们信号里有什么频率，以及它们的强度。相位谱则告诉我们这些频率成分的“起始位置”关系，这对于信号重建、滤波器设计、通信系统解调等都至关重要。

%% 5. 计算相位谱 phase = angle(X); % angle函数返回复数的相位角，单位弧度，范围[-π， π] % 同样，我们通常关心单边谱的相位 phase_single = phase(1:N/2+1); % 绘制相位谱 figure stem(f_single, phase_single, ‘LineWidth‘， 1.5) xlabel(‘频率 (Hz)‘) ylabel(‘相位 (弧度)‘) title(‘单边相位谱‘) xlim([0, Fs/2]) grid on

对于我们的理想合成信号（不含噪声），在50Hz， 120Hz等频率点上的相位应该是一个固定值。但由于我们信号是由sin函数生成的，其初始相位是0，而sin函数可以看作cos函数相位偏移-π/2。所以理论上，在这些频率点，相位值应该接近 -π/2（约-1.57弧度）。你可以检查一下谱峰处的相位值。

注意：angle函数返回的相位是“包裹”在[-π， π]区间内的。如果真实相位变化超过了这个范围，会发生相位跳变（从π跳到-π）。在分析连续变化的相位时（如振动分析中的相位差），可能需要使用unwrap函数来解开这种包裹，得到连续的相位曲线：phase_unwrapped = unwrap(phase);。

4. 频谱泄漏与加窗：如何让谱峰更“瘦”更准？

在上一节的理想例子中，我们的信号频率（50， 120， 200， 310）恰好是频率分辨率Δf（= Fs/N = 1000/500 = 2Hz）的整数倍。这种情况下，信号能量完美地集中在单一的频率点上，谱线又“瘦”又高。但在现实中，信号频率很少正好是Δf的整数倍。

4.1 频谱泄漏现象

让我们修改一下信号，把50Hz改成51Hz（不是2Hz的整数倍）。

%% 演示频谱泄漏 Fs = 1000; T = 0.5; t = 0:1/Fs:T-1/Fs; N = length(t); % 信号频率不是频率分辨率的整数倍 f_signal = 51; % Hz x_leak = sin(2*pi*f_signal*t); % 做FFT X_leak = fft(x_leak); P2_leak = abs(X_leak)/N; P1_leak = P2_leak(1:N/2+1); P1_leak(2:end-1) = 2*P1_leak(2:end-1); f_single = (0:(N/2)) * (Fs/N); figure subplot(2,1,1) stem(f_single, P1_leak, ‘LineWidth‘， 1.5) title([‘频谱泄漏示例：信号频率 = ‘， num2str(f_signal)， ‘Hz‘]) xlabel(‘频率 (Hz)‘) ylabel(‘幅度‘) xlim([40， 70]) grid on

你会发现，51Hz信号的频谱不再是一个干净的尖峰，而是像一座“小山”，能量“泄漏”到了旁边的频率点上。主峰变胖、变矮，旁边出现了许多不应该有的旁瓣。这会带来两个问题：1) 频率估计不精确；2) 强信号的小旁瓣可能会淹没附近弱信号的主峰，导致无法检测。

4.2 加窗函数的作用与选择

频谱泄漏的根本原因在于，我们对信号进行了“截断”。我们分析的是一段有限长的信号，这相当于用一个矩形窗去乘一个无限长的信号。矩形窗在时域是突然开始、突然结束的，其频谱有很高的旁瓣。这种时域的不连续性导致了频域的严重泄漏。

加窗，就是在做FFT之前，用一个窗函数（通常两端小、中间大）去乘原始信号，让信号的起始和结束部分平滑地过渡到零，从而减少截断带来的频谱泄漏。代价是主峰会进一步加宽（频率分辨率轻微下降），并且信号幅度会有微小衰减（需要修正）。

MATLAB提供了丰富的窗函数，如汉宁窗（hann）、汉明窗（hamming）、布莱克曼窗（blackman）等。

%% 应用汉宁窗 win = hann(N)‘; % 生成汉宁窗，转置成行向量 x_windowed = x_leak .* win; % 加窗 % 计算加窗信号的FFT X_win = fft(x_windowed); P2_win = abs(X_win)/N; % 注意：加窗后，幅度需要除以窗函数的相干增益进行补偿。 % 对于汉宁窗，相干增益约为0.5。更准确的做法是除以窗函数的能量（范数）。 ENBW = norm(win， 2)^2 / N; % 等效噪声带宽的一种计算 P2_win_corrected = abs(X_win) / (sum(win)); % 常用幅度修正方法1 % 或 P2_win_corrected = abs(X_win) / sqrt(mean(win.^2)*N); % 方法2 P1_win = P2_win_corrected(1:N/2+1); P1_win(2:end-1) = 2*P1_win(2:end-1); subplot(2,1,2) stem(f_single, P1_win, ‘LineWidth‘， 1.5， ‘Color‘， ‘r‘) title(‘加汉宁窗后的频谱‘) xlabel(‘频率 (Hz)‘) ylabel(‘修正后幅度‘) xlim([40， 70]) grid on

加窗后，虽然主峰更宽了，但旁瓣被显著抑制，频谱看起来更“干净”。这对于分析含有多个频率成分，尤其是强弱信号并存的场景非常有用。

如何选择窗函数？这是一个权衡：

• 矩形窗：频率分辨率最高（主瓣最窄），但旁瓣最高，泄漏最严重。适用于瞬态信号或精确已知周期的情况。
• 汉宁窗：旁瓣衰减好，频率分辨率中等。是最常用的通用窗，适合大多数频谱分析。
• 汉明窗：与汉宁窗类似，但第一个旁瓣更低，旁瓣衰减速度稍慢。
• 布莱克曼窗：旁瓣抑制最好，但主瓣最宽，频率分辨率最低。适用于需要极低旁瓣的场合。

实操心得：对于一般的频谱分析，我通常默认使用汉宁窗。除非有特殊理由（比如追求最高的频率分辨率，或者已知信号是同步采样的完整周期），否则不要用矩形窗。加窗后，一定要记得对幅度进行修正，否则所有频率分量的幅度都会偏低。MATLAB的信号处理工具箱（Signal Processing Toolbox）中的pwelch（功率谱密度估计）等函数已经内置了加窗和修正逻辑，在需要做谱估计时直接使用这些高级函数会更省心、更准确。

5. 功率谱密度：从幅度到能量视角

在很多工程应用，特别是噪声分析、振动测试、通信系统中，我们更关心信号功率在频域的分布，而不是单个频率点的幅度。这时就需要计算功率谱密度。

5.1 周期图法

最简单的方法是直接对幅度谱求平方，并考虑单边谱和系数。

%% 基于FFT的周期图法计算功率谱 x = comp1 + comp2 + 0.5*randn(size(t)); % 用带噪声的信号示例 X = fft(x); Pxx_raw = (abs(X).^2) / (N*Fs); % 双边功率谱密度估计 Pxx_single = Pxx_raw(1:N/2+1); Pxx_single(2:end-1) = 2 * Pxx_single(2:end-1); % 转换为单边 f_psd = (0:(N/2)) * (Fs/N); figure plot(f_psd, 10*log10(Pxx_single)) % 用dB表示 xlabel(‘频率 (Hz)‘) ylabel(‘功率/频率 (dB/Hz)‘) title(‘使用周期图法估计的单边功率谱密度‘) grid on

这种方法称为周期图法。但它估计的方差很大，曲线非常“毛糙”，不稳定。

5.2 韦尔奇方法：更优的PSD估计

韦尔奇方法是实际应用中的标准做法。它将长信号分成重叠的若干段，对每一段加窗并计算周期图，最后对所有段的周期图进行平均。这大大降低了估计的方差，得到了更平滑、更稳定的功率谱估计。

MATLAB中可以直接使用pwelch函数：

%% 使用pwelch函数（推荐） % 参数设置 window = hann(N/4); % 窗函数，段长度设为N/4 noverlap = length(window)/2; % 50%重叠 nfft = max(256， 2^nextpow2(length(window))); % FFT点数，至少256 [Pxx_welch， f_welch] = pwelch(x， window， noverlap， nfft， Fs); figure plot(f_welch, 10*log10(Pxx_welch)) xlabel(‘频率 (Hz)‘) ylabel(‘功率/频率 (dB/Hz)‘) title(‘使用Welch方法估计的功率谱密度‘) grid on

pwelch函数自动处理了加窗、重叠、平均、幅度修正等所有细节，返回的Pxx_welch就是估计的单边功率谱密度，其单位是原信号单位的平方每Hz（如 V²/Hz）。转换为dB单位后，可以更清晰地观察不同频率成分的相对强度。

注意事项：pwelch函数输出的频率向量f_welch只包含正频率部分（单边谱）。window的长度和noverlap的选择会影响结果：窗口越长，频率分辨率越高，但方差越大，曲线越不平滑；重叠越多，用于平均的段数越多，方差越小，曲线越平滑，但计算量也越大。通常选择50%的重叠是一个很好的折中。

6. 实战中的高频问题与调试技巧

掌握了基本原理和标准流程后，在实际项目中你还会遇到一些更具体、更棘手的问题。这里分享几个我踩过的坑和对应的解决方案。

6.1 频谱图中频率轴对不上？

这是最常见的问题之一。症状：明明输入一个100Hz的信号，谱峰却出现在200Hz或者别的莫名其妙的位置。

• 检查1：采样率Fs赋值是否正确？确保你构建时间向量t和计算频率向量f时使用的是同一个Fs。
• 检查2：频率向量计算公式是否正确？对于单边谱，f = (0:N/2) * (Fs/N)。确保N是FFT的长度（可能是补零后的NFFT）。如果你用了fftshift，频率向量也要相应地从负频率开始。
• 检查3：信号是否是实信号？如果你处理的是复数信号（如通信中的I/Q数据），那么频谱不是共轭对称的，不能简单取前半部分。你需要显示整个-Fs/2到Fs/2的双边谱。

一个可靠的频率向量生成模板：

Fs = your_sample_rate; N = length(your_signal); % 或你指定的NFFT f = (0:N-1)*(Fs/N); % 双边谱频率 (0 到 Fs) f_shift = (-N/2:N/2-1)*(Fs/N); % 用于fftshift后的频率 (-Fs/2 到 Fs/2) f_single = (0:N/2)*(Fs/N); % 单边谱频率 (0 到 Fs/2)

6.2 幅度谱的幅度不对？

谱峰找到了，频率也对，但幅度和信号的实际幅值对不上。

• 根本原因：归一化因子错误。回顾第3节，核心公式要记牢：
  • 双边幅度谱（未修正）：magnitude = abs(X)
  • 恢复真实幅度的双边谱：对于非直流分量，true_magnitude_double = 2 * abs(X) / N
  • 恢复真实幅度的单边谱：true_magnitude_single = 2 * abs(X(1:N/2+1)) / N，并且true_magnitude_single(1)对应直流，不乘2；true_magnitude_single(end)对应奈奎斯特频率，如果N是偶数也不乘2（但通常我们只显示到N/2）。
• 加窗后的修正：如果加了窗，幅度会衰减。修正因子通常是窗函数的能量和或相干增益。对于pwelch等高级函数，内部已做修正。手动修正可以参考第4.2节的代码。

6.3 如何精确测量频率和相位？

当信号频率不是频率分辨率的整数倍时，直接取谱峰对应的频率和相位会有误差。

• 频率插值法：可以通过谱峰附近几个点的幅度，进行抛物线或重心法插值，来估计更精确的频率。MATLAB信号处理工具箱中的findpeaks函数可以结合插值选项使用。
• 相位测量：直接从angle(X(k))读取的相位，对噪声和频谱泄漏非常敏感。一种更稳健的方法是：phase = atan2(imag(X(k))， real(X(k)))。对于高精度需求，可以考虑使用基于解析信号的希尔伯特变换方法，或者专门的正弦波拟合算法。
• 整周期采样：在条件允许的情况下，尽量使采样时长包含信号周期的整数倍。这样可以完全避免频谱泄漏，获得最精确的幅度和相位。这需要事先知道或估计信号的主频率。

6.4 处理大数据量时的性能与内存

当信号长度N非常大（例如上百万点）时，直接fft(x)可能会消耗大量内存和计算时间。

• 分段处理：使用pwelch本身就是一种分段平均。对于其他需要全数据FFT的操作，可以考虑先下采样（如果高频信息不重要），或者使用迭代/分段的方法。
• 使用fft(X, [], dim)指定维度：如果你的数据是多通道的（例如多路传感器数据），确保沿正确的维度（通常是列）进行FFT，避免无谓的循环。
• GPU加速：对于超大规模计算，如果拥有Parallel Computing Toolbox和兼容的GPU，可以使用gpuArray将数据送入GPU，然后使用fft，速度会有数量级的提升。例如：x_gpu = gpuArray(x); X_gpu = fft(x_gpu); X = gather(X_gpu);。

6.5 从仿真工具（如Vivado）导出数据给MATLAB分析

从Vivado Simulation或Xilinx FFT IP核导出的数据，常常需要预处理。

• 数据格式：导出的数据可能是二进制、十六进制文本或.csv文件。使用fscanf、textscan或readmatrix读取。
• 复数处理：FFT IP核的输出通常是分离的实部（I）和虚部（Q）数据。你需要将它们组合成复数：data_complex = I + 1j*Q;。
• 位宽与定点数：导出的数据可能是定点数，带有特定的位宽和小数点位。你需要根据IP核的配置，将其转换为MATLAB中的浮点数。例如，如果输出是ap_fixed<16，14>（总共16位，14位整数），在MATLAB中可能需要除以2^(16-14)或进行类似的缩放。
• 顺序问题：如热词中提到的“vivado中fft核输出iq反了”，这通常是因为对输出数据格式的理解有误。仔细阅读IP核文档，确认输出是I+ jQ还是Q + jI，以及输出是自然顺序还是倒位序。Xilinx FFT IP核通常支持多种输出顺序，需要在配置时和读取时保持一致。如果顺序不对，可以使用fftshift或自己编写索引重排代码进行调整。

7. 超越基础：几个进阶应用场景

掌握了单信号分析后，FFT在MATLAB中还能玩出更多花样。

7.1 使用fft2进行二维图像频率分析

FFT可以推广到二维，用于图像处理。图像的二维FFT反映了图像在水平和垂直方向上的空间频率成分。低频对应图像中平缓变化的区域（如背景），高频对应边缘和细节。

%% 图像二维FFT示例 img = imread(‘cameraman.tif‘); % 读取灰度图像 img_double = im2double(img); % 转换为双精度 F = fft2(img_double); % 二维FFT F_shifted = fftshift(F); % 将零频移到中心 magnitude_spectrum = log(1 + abs(F_shifted)); % 对数变换便于显示 phase_spectrum = angle(F_shifted); figure subplot(1,3,1)， imshow(img)， title(‘原图‘) subplot(1,3,2)， imshow(magnitude_spectrum， [])， title(‘幅度谱（对数）‘) subplot(1,3,3)， imshow(phase_spectrum， [-pi pi])， title(‘相位谱‘) colormap gray

通过修改幅度谱或相位谱，再进行逆变换ifft2，可以实现图像滤波、压缩等操作。

7.2 使用goertzel函数进行单频点能量检测

如果你只关心少数几个特定频率（例如DTMF电话拨号音解码），使用完整的FFT计算所有频率点是浪费的。Goertzel算法是一种高效的递归算法，用于计算DFT在单个或多个特定频率点上的值。

%% 使用Goertzel算法检测特定频率 Fs = 8000; t = 0:1/Fs:0.1-1/Fs; x = 0.5*sin(2*pi*697*t) + 0.5*sin(2*pi*1209*t); % DTMF信号“1” % 我们只关心DTMF的行频和列频 target_freqs = [697, 770, 852, 941, 1209, 1336, 1477]; N = length(x); indices = round(target_freqs * N / Fs) + 1; % 对应的DFT索引（从1开始） % 使用Goertzel for k = 1:length(target_freqs) idx = indices(k); % goertzel函数需要信号和频率索引（从0到N-1） det = goertzel(x， idx-1); % idx-1 转换为从0开始的索引 magnitude(k) = abs(det) * 2 / N; % 计算幅度 end figure stem(target_freqs， magnitude) xlabel(‘频率 (Hz)‘) ylabel(‘检测幅度‘) title(‘使用Goertzel算法检测DTMF频率‘) grid on

可以看到，在697Hz和1209Hz处有明显的峰值，对应按键“1”。

7.3 使用spectrogram函数绘制时频谱图

对于非平稳信号（频率随时间变化，如鸟叫声、雷达信号），单纯的FFT会丢失时间信息。短时傅里叶变换通过一个滑动的窗，对信号分段进行FFT，从而得到信号频率随时间变化的图谱，即频谱图。

MATLAB中的spectrogram函数可以一键生成。

%% 生成并分析一个频率线性变化的信号（啁啾信号） Fs = 1000; t = 0:1/Fs:2; x = chirp(t， 0， 1， 250); % 频率从0Hz线性增加到250Hz figure spectrogram(x， 256， 250， 256， Fs， ‘yaxis‘) % 窗长256，重叠250，FFT点数256 title(‘啁啾信号的时频谱图‘) colorbar

图中，颜色代表能量强度，纵轴是频率，横轴是时间。可以清晰地看到一条从低频斜向高频的亮线，这就是频率的变化过程。spectrogram的参数（窗长、重叠、FFT点数）需要根据信号特性调整，以在时间分辨率和频率分辨率之间取得平衡。