做信号处理的人几乎都绕不开傅里叶变换。我刚接触这个领域时第一件事就是按教材公式在 MATLAB 里写了一个“极简版”的离散傅里叶变换DFT然后拿它和内置的fft函数对比结果同一个 4096 点的序列我的函数跑了将近一秒钟fft几乎是瞬间出结果。后来我才明白DFT 是一套数学定义而 FFT 只是让 DFT 变得可实际使用的高效算法。这篇博客我就用 MATLAB 把 DFT、IDFT、FFT 这三件事完整串一遍先手写最直观的 DFT 和 IDFT验证它们和 MATLAB 内置函数的结果一致性再拆解 FFT 的加速思路和调用细节最后用频谱分析、频域滤波、FFT 卷积三个实际例子收尾让初学者既能看懂公式也能直接上手做项目。1. 从“离散采样”到“频谱世界”DFT到底在算什么1.1 计算机里没有连续信号所以必须离散化傅里叶变换的原始定义建立在连续时间函数上它把时域信号拆成一系列复指数分量的叠加。但在计算机里我们拿到的永远是一串采样点比如音频采集卡的输出、示波器抓到的波形、传感器每隔固定时间记录的数据根本没有一个连续函数可以直接积分。这就需要一个“能放进计算机”的傅里叶变换版本把积分变成求和把连续变量变成离散序号DFT 就是干这个的。DFT 的输入是一个长度为 N 的序列 x(0), x(1), ..., x(N-1)输出是同样长度为 N 的频谱序列 X(0), X(1), ..., X(N-1)。每个 X(k) 代表这个信号在频率 k 上的复振幅。公式非常简洁X(k) Σ_{n0}^{N-1} x(n) * W_N^{kn}其中 W_N e^{-j2π/N}x(n) (1/N) * Σ_{k0}^{N-1} X(k) * W_N^{-kn}这里 W_N 被称为旋转因子是整个 DFT 的核心。它本质上是一个落在单位圆上的复数随着 n 和 k 的变化它沿着单位圆转圈。当 k 比较小时它转得慢代表低频k 接近 N/2 时它转得快代表高频。你可以把 DFT 理解成用一组不同转速的“探针”去试探信号每个探针和信号相乘再累加如果信号里有和这个探针同频的分量累加结果就很大如果没有累加结果就接近于零。1.2 DTFT、DFT、DTFS别再傻傻分不清网上经常有人把 DTFT 和 DFT 混着说MATLAB 帮助文档里偶尔也能看到让人困惑的表述。这里用一张表把容易混淆的三种变换关系理清楚名称时域频域用途DTFT离散时间傅里叶变换无限长离散序列连续频谱理论分析DFT离散傅里叶变换有限长离散序列有限长离散频谱计算机计算DFS离散傅里叶级数周期离散序列周期离散频谱周期信号分析DTFT 是无限长序列的傅里叶变换频谱是连续的没法直接放进计算机。DFT 相当于对 DTFT 的频谱做等间隔采样只取 N 个点所以它天然隐含了“时域周期化”的假设——你在计算机里对一段有限长信号做 DFT实际上是在分析这段信号当做一个无限周期信号的一个周期。这也是后面频谱泄漏的根源之一第四节会详细说。1.3 核心结论DFT 是可计算、可逆、可验证的DFT 最重要的性质是“可逆”正变换把时域变到频域逆变换把频域变回时域中间只差一个 1/N 和指数符号的翻转。这意味着你完全可以用 DFT 做滤波先把信号变换到频域把某些频率成分乘 0 或变小再变换回时域。这里提前给一个结论方便后面验证代码对一个长度为 N 的序列 xifft(fft(x))的结果应该和 x 完全一致误差在机器精度级别大约 10^-14 量级。如果你的误差很大那一定是归一化或者符号写错了。这是检查 DFT/IDFT 实现是否正确的黄金标准。2. 手写DFT与IDFT两段MATLAB代码验证数学公式2.1 用旋转因子矩阵实现DFT一次矩阵乘法搞定全部频率点最直观的 DFT 写法是按公式写双重循环外层循环遍历频率下标 k内层循环遍历时间下标 n复杂度是 N^2。对于 N1024 这种小规模数据还好到了 N8192 就会等很久而且代码啰嗦。MATLAB 更擅长矩阵运算所以我建议直接把旋转因子变成一个矩阵用一次矩阵乘法完成全部计算。function X myDFT(x) % 输入: x 为长度为 N 的列向量或行向量 % 输出: X 为长度为 N 的 DFT 结果列向量 N length(x); n (0:N-1).; % 列向量: N x 1 k 0:N-1; % 行向量: 1 x N W exp(-1j * 2 * pi * n * k / N); % 旋转因子矩阵: N x N X W * x(:); % 矩阵乘列向量一行代码完成双重求和 end这里面的精髓在于n * k这个外积。n 是列向量k 是行向量两者相乘得到一个 N×N 矩阵第 (n, k) 个元素正好是 n*k再放进指数函数里就得到了公式需要的旋转因子矩阵 W。每一行对应一个时间索引 n每一列对应一个频率索引 k矩阵乘法本质上是在同时完成 N 组求和。相比之下平时做频谱分析更常用的是fft也就是 MATLAB 内置的快速傅里叶变换函数。它的输出和myDFT是完全一样的数学对象只是算法速度更快。验证方式很简单fs 1000; % 采样率 1000 Hz t (0:1023) / fs; % 1 秒时长1024 个点 x sin(2*pi*50*t) 0.5*cos(2*pi*120*t pi/6); X1 myDFT(x); X2 fft(x); max(abs(X1(:) - X2(:))) % 应该约为 1e-13 甚至更小实测下来这个误差通常在 10^-13 量级甚至更小。如果你看到误差是 10^-3 量级那说明旋转因子矩阵构建时有下标错误最常见的问题是把 0 到 N-1 写成了 1 到 N。2.2 IDFT实现指数符号和归一化系数一个都不能错逆变换是把频域变回时域公式里的指数从负变成正外面乘了 1/N。代码结构和正变换完全对称function x myIDFT(X) % 输入: X 为长度为 N 的 DFT 结果列向量 % 输出: x 为恢复出的时域序列列向量 N length(X); n (0:N-1).; k 0:N-1; W exp(1j * 2 * pi * n * k / N); % 注意这里指数是正号 x (1 / N) * W * X(:); end有一个细节值得单独说明为什么正变换不乘 1/N逆变换才乘这是 DFT 的一种常见归一化约定。MATLAB 的fft和ifft也是这样配套的正变换不乘逆变换乘 1/N两者复合正好回到原始信号。如果你看到别的教材把 1/N 放在正变换里那也合理只是约定不同但正逆变换必须配套否则来回变换之后幅度会差 N 倍。我建议始终以 MATLAB 的约定为准因为实际工程中你大概率会混用内置函数和手写代码。验证时我通常会同时检查两件事x_rec myIDFT(X2); % 1. 恢复信号的虚部应该接近 0 max(abs(imag(x_rec))) % 2. 恢复信号和原始信号一致 max(abs(real(x_rec) - x(:)))第一个检查很关键原始信号是实数时经过正变换得到复数频谱再经逆变换恢复出来的信号的虚部应该在 10^-15 量级。如果虚部很大说明 IDFT 的指数符号或归一化写错了。2.3 用Parseval定理做交叉验证除了“正逆变换互相还原”这个检查我还习惯用帕塞瓦尔定理做第二次验证。它说的是信号在时域的总能量等于频域的总能量。对 DFT 而言具体形式是sum(abs(x).^2) (1/N) * sum(abs(X).^2)随手写一个验证energy_time sum(abs(x).^2); energy_freq sum(abs(X2).^2) / length(x);这两个值理论上完全相等实际误差同样在机器精度级别。如果不等请检查是否有数据被截断、混叠或者归一化出错。帕塞瓦尔定理的好处是它只依赖总能量不依赖某个具体频率点的细节所以它能发现那种“某些点大一点、另一些点小一点”的隐性错误而时域还原检查更容易发现整体性错误两者互补。3. FFT的加速原理与MATLAB调用要领3.1 从 N² 到 Nlog₂N蝶形运算到底做了什么DFT 的矩阵实现虽然优雅但复杂度是 N²。对一个 N4096 的序列矩阵法要做 4096×4096≈1678 万次复数乘法。就算 MATLAB 矩阵运算已经高度优化这个规模也会明显卡顿。FFT 的突破在于利用旋转因子的两个性质对称性 W_N^{kN/2} -W_N^k周期性 W_N^{kN} W_N^k把一个大 DFT 拆成多个小 DFT。最经典的是基2时间抽取算法。设 N 是 2 的整数次幂先把序列按奇偶下标分成两组x_even 和 x_odd那么 N 点 DFT 可以写成两个 N/2 点 DFT 的组合。这个过程不断递归直到变成 2 点 DFT。每一步递归中的合并操作叫做“蝶形运算”因为它的计算图看起来像一只蝴蝶两个输入、两个输出一个加法一个减法再乘一个旋转因子。总共有 log₂N 级每级 N/2 个蝶形所以总计算量大约是 (N/2)log₂N 次复数乘法。用具体数字对比更直观N直接DFT复数乘法次数FFT复数乘法次数加速比1024约 104 万约 5 千约 200 倍4096约 1678 万约 2.5 万约 680 倍65536约 43 亿约 52 万约 8200 倍数据越长FFT 的优势越夸张。这也是为什么实时频谱分析、音频处理、雷达信号处理这些场景里FFT 是绝对的基础设施。3.2 MATLAB中fft、ifft、fftshift、ifftshift的正确打开方式MATLAB 里调用 FFT 极简单X fft(x)就完成了ifft(X)做逆变换。但实际项目中有几个容易忽略的细节第一fft(x, Nfft)这种带第二参数的写法非常常见。如果Nfft大于length(x)MATLAB 会自动在序列尾部补零如果小于length(x)会自动截断。零填充在后面第四节会重点讨论这里记住它不会改变信号本身只改变频谱采样的密度。第二fft输出的频率顺序是 0 到 Fs然后才是负频率。具体来说X(1)对应直流分量 0HzX(2)对应 Fs/N直到X(N)对应 -Fs/N。如果你想画一个以 0Hz 为中心的对称频谱图就需要fftshiftf_shift (-N/2 : N/2 - 1) * (fs / N); X_shift fftshift(X); plot(f_shift, abs(X_shift));第三ifftshift和fftshift并不是完全可互换的。对于偶数长度的序列两者效果一样对于奇数长度的序列ifftshift是fftshift的逆运算方向相反。做频域滤波时如果你先把频域结果fftshift处理完再ifft之前应该用ifftshift把它搬回原来的顺序而不是再用一次fftshift。奇数长度时用错了会直接导致相位错乱时域波形变得面目全非。3.3 实测对比矩阵版DFT vs 内置FFT上点实测数据。我用一个 N4096 的随机信号对比myDFT和fft的耗时N 4096; x randn(N, 1); tic; X1 myDFT(x); t_dft toc; tic; X2 fft(x); t_fft toc; fprintf(矩阵版DFT耗时: %.4f 秒\n, t_dft); fprintf(内置FFT耗时: %.6f 秒\n, t_fft); fprintf(加速比: %.0f 倍\n, t_dft / t_fft);在我自己的笔记本上myDFT大约需要 0.3 到 0.5 秒fft在 1 毫秒以内加速比上百倍。这个差距正好对应 N² 和 Nlog₂N 的复杂度差异。不过要提一句这种手写矩阵版比真正的双重循环版已经快了很多因为 MATLAB 的矩阵乘法底层调用了高度优化的线性代数库。如果你用双重循环逐个元素乘慢得会更明显。所以我的建议是教学和理解用矩阵版生产环境永远用fft。不要把“自己实现 FFT”和“自己造轮子”划等号除非你是为了学习或者特殊硬件需要否则直接用内置函数是更专业、更稳妥的选择。4. 工程中最容易翻车的四个细节4.1 频率轴标定横坐标算错频谱白看fft给出的结果本身不含任何“频率”含义它只有下标 k你必须自己换算成物理频率。换算公式是f(k) k * Fs / Nk 0, 1, ..., N-1频率分辨率相邻两个频率采样点之间的间隔是 Fs/N。举个例子采样率 Fs1000HzN1024 点那么每个频率点间隔约 0.9766Hz。如果你不除以 N直接用 k 当横坐标那 50Hz 的信号会出现在第 50 根线或者后面的位置和真实情况对不上整个分析就失去了物理意义。画单边频谱时通常只取前 N/21 个点对应频率从 0Hz 到 Nyquist 频率 Fs/2。画双边频谱时用fftshift后取整个 N 个点。两者的横轴向量不同千万别混用。我把这两个向量的构建方式写在下面方便直接抄% 单边谱频率轴 f_one_sided (0:N/2) * (Fs / N); % 双边谱频率轴偶数 N f_two_sided (-N/2 : N/2 - 1) * (Fs / N);4.2 幅值归一化为什么峰值不是信号的真实幅值这是初学者最容易困惑的问题。对一个幅值为 A 的正弦信号 sin(2πft)它的 DFT 结果在对应频率点上的幅度是 A*N/2而不是 A。为什么多出 N/2 这个系数因为正弦信号在频域被拆成了两个镜像分量正频率和负频率各占一半能量而每个分量又被 DFT 的求和过程“放大”了 N 倍N 个样本累加。所以如果你关心的是信号的幅度谱要恢复真实幅值需要把对应频率点的幅值乘以 2/N。直流分量k0比较特殊它本身就是真实平均值放大 N 倍所以只除以 N 就行。我一般这样封装X fft(x); % 单边幅度谱 A abs(X(1:N/21)) / N; A(2:end-1) A(2:end-1) * 2; % 除直流外都乘 2 f (0:N/2) * (Fs / N);相位谱也有类似的问题。angle(X)给出的是弧度制的相位对于实信号正频率对应的相位就是该频率分量在时域波形中的初相但要记得把负频率对称点的影响排除在外。还要注意当某频率点幅值接近零时相位是无意义的随机数别画出来吓唬自己。4.3 零填充与频谱泄漏提高“颜值”不等于提高分辨率零填充Zero Padding是fft(x, Nfft)最常见的应用场景。它的作用是让频谱看起来更平滑因为补零后 FFT 的输出点数变多了频谱采样更密集。但要注意一个关键区别谱线间隔变密不等于频率分辨率变高。频率分辨率真正取决于信号的实际观测时长 T N/Fs分辨率的下限是 1/T。你补了一堆零观测时长并没有变长所以并不能把两个靠得很近的频率成分分开。比如 100Hz 和 100.5Hz 的两个正弦波观测 1 秒时分辨率刚好 1Hz勉强能看出是两个峰补零到 2048 点谱线更密了但 100Hz 和 100.5Hz 的两个峰还是糊在一起。频谱泄漏和零填充是两个不同的问题。当你截取一段有限长信号做 DFT 时相当于在时域乘了一个矩形窗。如果截取长度不是信号周期的整数倍这个矩形窗的频谱会有很高的旁瓣导致原来的单根谱线变成一片“拖尾”能量泄漏到附近的频率。缓解办法是加窗比如hann、hamming、blackman。实际项目中我的习惯是先做时域分析判断信号里可能存在哪些频率成分再决定是否需要零填充、选什么窗。加窗虽然能降低泄漏但也会加宽主瓣降低频率分辨率这是不可避免的折中。4.4 相位信息与unwrap别被angle骗了很多初学者只画幅度谱完全忽略相位谱这在某些场景下会出大问题。相位信息在频域滤波、系统辨识、时延估计里都极其重要。比如你要构造一个逆滤波器如果忽略相位只保留幅度滤波后的波形会完全对不上。angle返回的相位范围是 -π 到 π也就是说它被包裹在 [-π, π] 区间内。如果真实相位是 3πangle显示的是 -π 或者接近 -π直接拿去计算会出问题。这时候用unwrap把相位展开成连续曲线phase unwrap(angle(X));还有一个容易误解的地方相位谱只有在信号是窄带、且频率分量幅值较大时才有意义。对宽带噪声或者幅值接近零的频率点相位是随机的不要尝试从里面提取信息。这里整理一个易错点对照表是我自己踩过坑后总结出来的问题表现原因处理办法横坐标是下标而不是频率峰值位置对不上没有乘 Fs/N用 f k*Fs/N峰值是 A*N/2幅度谱偏高没做归一化除 N再乘 2除直流补零后“分辨率提高了”两个近频峰还是重叠零填充不改变物理分辨率延长观测时间相位跳跃剧烈相位谱锯齿状相位被包裹用 unwrapifft前用fftshift而不是ifftshift时域波形错乱奇数长度时两者不同用 ifftshift 恢复顺序5. 完整实例从频谱分析到频域滤波再到FFT卷积5.1 合成一段多分量信号并做频谱分析纸上谈兵没意思我们来做一个完整的闭环练习。假设一段采集到的信号采样率 Fs1000Hz持续 1 秒里面包含一个 50Hz 的正弦波、一个 120Hz 的余弦波以及一定的随机噪声。我们用 MATLAB 生成并分析它Fs 1000; T 1; N Fs * T; t (0:N-1) / Fs; x 1.5 * sin(2*pi*50*t) 0.8 * cos(2*pi*120*t pi/4) 0.1 * randn(1, N); X fft(x); f (0:N/2) * (Fs / N); A abs(X(1:N/21)) / N; A(2:end-1) A(2:end-1) * 2; % 找峰值及其对应频率 [peak_vals, peak_idx] sort(A, descend); for i 1:4 fprintf(峰值 %.3f %.1f Hz\n, peak_vals(i), f(peak_idx(i))); end运行这段代码你应该能看到两个明显的峰值50Hz 附近接近 1.5120Hz 附近接近 0.8噪声成分出现在整个频带里。这个练习虽然简单但它验证了 DFT 的幅值恢复、频率轴标定和信号模型之间的全部逻辑链。5.2 频域做低通滤波再ifft回到时域现在的目标是滤掉 120Hz 成分只保留 50Hz。最直接的办法是在频域构造理想低通滤波器把大于某个截止频率的部分直接置零然后ifft回时域。注意卷积定理在这里的运用时域卷积等效于频域相乘所以滤波就是“频谱逐点相乘”。% 构造理想低通滤波器截止频率 80 Hz H zeros(1, N); H(abs(freq_axis) 80) 1; % freq_axis 需要完整双边的频率轴 Y X .* H; % 频域相乘 y_filtered ifft(Y); % 逆变换回时域如果H是理想矩形ifft后的信号会看到吉布斯效应——在信号起始和结束附近出现明显的振铃振荡这是因为矩形滤波器在时域的脉冲响应是 sinc 函数和有限长信号卷积时在边界处剧烈波动。实际工程里通常用过渡带平滑的滤波器形状比如高斯型、巴特沃斯近似等来缓解这个问题。这里我不展开讲窗函数设计但请牢记理想滤波器在数字信号处理里只是理论工具直接使用会有边界伪影。这个例子还演示了一个重要概念你在频域的任何修改最终都会体现在时域波形上而频域的突变越剧烈时域产生的伪影越明显。这就是为什么频域滤波不能“想怎么改就怎么改”。5.3 用FFT快速计算线性卷积补零长度是关键最后来看 FFT 的另一个经典应用用 FFT 加速线性卷积。时域线性卷积的复杂度是 O(N×M)如果两个序列都很长计算量非常大。但根据卷积定理线性卷积可以通过 FFT 在频域完成复杂度是 O(Nlog₂N)。关键问题是直接fft(x, N) .* fft(h, N)得到的是循环卷积不是线性卷积两者在时域上会有混叠差异。要让圆卷积等于线性卷积必须对两个序列都补零到长度 L ≥ NM-1。这里用一个简单的例子验证x randn(1, 64); h randn(1, 32); % 直接时域卷积基准结果 y_conv conv(x, h); % 用FFT做线性卷积 L length(x) length(h) - 1; y_fft ifft(fft(x, L) .* fft(h, L)); % 对比误差 max(abs(y_conv - y_fft))max(abs(...))的结果应该在 10^-14 量级这证明 FFT 卷积和时域卷积在数值上一致。如果 L 取小了也就是小于 NM-1你会看到误差变得很大那是因为循环卷积把尾部的样本“卷”回了头部发生了时域混叠。FFT 卷积在实际工程中广泛用于大尺寸滤波器的快速计算、图像处理中的卷积神经网络前向计算、以及通信系统中的信道估计等场景。一旦序列长度超过 512 点FFT 方法通常就比直接卷积更快。这里的经验法则是当 N 和 M 都大于 64 时FFT 卷积开始有优势当长度到几千甚至几万时优势是数量级的。我在实际工作中最常犯的一个错误是在做完fftshift后忘记在ifft前ifftshift回来尤其是处理奇数长度信号时结果波形完全对不上。后来我给自己定了一条铁律凡是fft之前用了fftshift的ifft之前就必须用ifftshift并且每次拿到ifft结果后先看虚部和实部是否自洽再进入下一步分析。DFT、IDFT、FFT 这三个名词看起来很接近但只有你亲手写一遍正变换、逆变换再拿内置函数做对比、做滤波、做卷积你才会真正感觉到它们是如何作为一个整体在信号处理里运转的。