ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

离散傅里叶变换DFT/IDFT:从数学公式到频谱分析工程实践

离散傅里叶变换DFT/IDFT:从数学公式到频谱分析工程实践 很多人在学“信号处理”时第一次被劝退就是在离散傅里叶变换这一节。当年我啃这部分的时候公式里的求和符号、复指数、下标k和n每一个都认识放到一起就成了天书。后来在雷达信号处理项目里被频谱分析反复折腾才慢慢把DFT/IDFT从“数学定义”变成了“手上的工具”。这篇东西不是复述教科书而是把DFT和IDFT从公式到代码、从理论到工程那层窗户纸捅破顺便把那些网上说法各异、教材里又不爱写的坑也一并填了。这篇文章适合刚接触数字信号处理的学生、做嵌入式或通信的工程师也包括那些要用频谱分析但不想被教程绕晕的硬件开发者。读完你至少能明白三件事DFT到底在做怎样的数学运算、IDFT为什么能完美复原信号、以及在实际代码里怎么避开那些让你结果“看起来不对劲”的常见陷阱。1. 从连续傅里叶到离散DFT计算机眼中的频谱1.1 连续傅里叶变换的理想与现实理论上一个连续时间信号x(t)的傅里叶变换长这样X(f)\int_{-\infty}^{\infty} x(t)e^{-j2\pi f t} dt积分上下限全体实数被积函数连续这在推导公式时很完美。但拿到现实里计算机没法处理连续积分也没法保存无限长的波形。示波器采回来的信号永远是一串有限长度的离散点比如采样率8kHz、采了1024个点这就是你手上全部信息。于是问题变成这1024个点构成的有限长序列能不能得到一个类似“频谱”的东西DFT就是干这个的。1.2 从DTFT到DFT频域也被采样了教科书还会提到离散时间傅里叶变换DTFTX(e^{j\omega})\sum_{n-\infty}^{\infty} x[n]e^{-j\omega n}DTFT的输入是离散序列但频率变量ω连续输出的频谱还是实变函数。计算机同样没法保存一条连续频域曲线所以要在频域轴上等间隔打点只取N个频率点来观察。这一步就是频域采样。频域采样之后时域会发生什么一句话时域序列会被周期延拓。DFT正是“时域有限长频域有限长”的组合两个域上有限长也意味着两个域都在隐含地周期性重复。这解释了后面很多工程现象比如循环卷积。1.3 从连续到DFT的直观流程图很多人一上来就被公式砸晕其实不看数学也能理解DFT干了什么拿一个长度为N的离散序列把它看成是某个周期信号的一个周期用N个等间隔的频率点去试探这个序列每个频率点算出一个复数表示“这个频率的成分有多大、相位是什么”这就像用一把有N根筛齿的梳子把信号里不同频率的分量梳出来。每一根齿对应一个频率。频率梳的齿间距跟信号长度N有关齿的位置和采样率有关。2. DFT公式拆解一个求和到底算出了什么2.1 公式里每个符号都是什么意思先上标准定义X[k]\sum_{n0}^{N-1} x[n]e^{-j 2\pi kn/N},\quad k0,1,\dots,N-1这里x[n]是时域第n个采样点X[k]是频域第k个频率点N是序列长度。e^{-j2πkn/N}是一个复指数可以看成单位圆上旋转的向量。k固定时它像一个“频率梳子”n是变量表示在时域上一步步旋转采样。换个角度X[k]等于序列x[n]与复指数e^{-j2πkn/N}的“点积”。如果x[n]里恰好含有与这个复指数同频率的成分点积结果就很大如果不含有点积结果接近0。所以DFT的本质是一组“模板匹配”用N个不同频率的模板逐个去和信号做内积。这跟我之前做匹配滤波器的思路一模一样只不过匹配滤波器是主动设计的模板DFT的模板就是一组等间隔频率的复指数。2.2 用欧拉公式理解正负频率再看公式里的复指数根据欧拉公式展开e^{-j 2\pi kn/N}\cos(2\pi kn/N)-j\sin(2\pi kn/N)所以X[k]是序列x[n]分别与余弦、正弦做相关运算后的组合余弦相关对应实部正弦相关对应虚部。这也就解释了为什么实数信号的双边频谱是共轭对称的因为k和N-k对应的余弦相同、正弦相反两个X[k]自然就成了共轭。理解这层之后就不会再问“为什么负频率也有值”这种问题了。负频率在数学上存在在物理上可以理解成旋转方向相反的复指数分量。2.3 拿个具体序列走一遍假设有一个4点序列x [1, 2, 3, 4]N4。我们算X[0]X[0]\sum_{n0}^{3}x[n]123410这就是直流分量0Hz的大小。再算X[1]X[1]1*e^{-j0} 2*e^{-j\pi/2} 3*e^{-j\pi} 4*e^{-j3\pi/2}1 2*(-j) 3*(-1) 4*(j) -2 2jX[2]X[2]1 2*e^{-j\pi} 3*e^{-j2\pi} 4*e^{-j3\pi}1 - 2 3 - 4 -2X[3]X[3]1 2*e^{-j3\pi/2} 3*e^{-j3\pi} 4*e^{-j9\pi/2}1 2j - 3 - 4j -2 - 2j可以看到X[1]与X[3]共轭这正是实数序列的频谱特性。这套手算过程虽然慢但能帮你把求和符号从“抽象符号”变成“运算动作”。3. IDFT为什么除以N怎么完美复原3.1 逆变换公式与DFT的对称性DFT把时域N个点变成频域N个点信息量没有丢失所以理论上可以用N个X[k]恢复出原来的N个x[n]。逆变换IDFT定义为x[n]\frac{1}{N}\sum_{k0}^{N-1} X[k]e^{j 2\pi kn/N},\quad n0,1,\dots,N-1注意指数符号变成了正号前面多了个1/N。这个1/N是整个复原的“加权平均”系数。为什么必须有它简单说DFT相当于把原始信号“投影”到N个正交基上每个基的模长都是√N。为了还原原坐标需要把投影结果除以基的模长平方也就是N。3.2 用正交性理解复原过程更严谨地说复指数序列e^{j2πkn/N}在k0~N-1这N个基向量两两正交。学过线性代数都知道正交基下坐标还原就是每个基上的投影除以基的内积模长平方。对任意两个不同的频率k1和k2\sum_{n0}^{N-1} e^{j2\pi k1 n/N} e^{-j2\pi k2 n/N}\sum_{n0}^{N-1} e^{j2\pi (k1-k2)n/N}当k1≠k2时这个等比数列求和等于0当k1k2时等于N。这就是正交性的数学表述。所以把X[k]代回IDFT时k那一项只会从x[n]中抽出对应频率的分量其他分量全部抵消最后剩下N×x[n]再除以N就得到原序列。我当时学到这里才恍然大悟DFT/IDFT不是两个分离的算法而是一对变换对频域不过是对同一份信息换了个坐标系描述。3.3 一个快速验证用2点序列试x [2, 5]N2。直接按公式X[0]257 X[1]2*e^0 5*e^{-j\pi} 2 - 5 -3再IDFTx[0](7*1 (-3)*1)/2 2 x[1](7*1 (-3)*(-1))/2 5两次变换原序列原样回来。很多迷糊点只要亲自推一遍这种微型例子就能消除。4. DFT在信号处理里的经典坑泄漏、栅栏与卷积混淆4.1 频谱泄漏当被截断的正弦不再是“整周期”实际处理信号时你拿到的永远是有限长的一块数据相当于对无限长信号做了一个矩形窗截断。如果截断长度不是信号周期的整数倍频谱就会“糊”掉原本一根干净频谱线变成一坨展宽的主瓣和一堆旁瓣。这就是频谱泄漏。举个例子一个50Hz正弦波采样率1000Hz采样点数100恰好是5个完整周期DFT后会在50Hz处得到一个干净的冲击峰。如果采样点数改成97第97个点不在周期结束位置频谱里就会出现很多不该有的低频分量。原因就是矩形窗在边界处制造了不连续性等价于给信号附加了高频成分。4.2 补零能“插值”但提高不了分辨率很多人发现频谱“不够平滑”第一反应是补零把N从1000补到4000。这样DFT频域点数也变成了4000曲线确实更细腻。但注意补零不会让物理分辨率变高它只是把原本1000个频率点的频谱用sinc插值到4000个点。真实的分辨率由信号的有效时长决定近似等于fs/N。如果两个频率差异小于1/T补多少零都没用。所以工程上要分辨两个邻近频率正确手段是采集更多时间的数据而不是在尾部加零。我做雷达多目标测距时就吃过补零的亏。以为补零能让峰值更“尖”结果只是曲线变光滑两个目标的峰值照样重合。后来延长了积分时间问题才解决。4.3 窗函数不是玄学是工程刚需为了抑制频谱泄漏可以在做DFT之前给序列乘一个窗函数。矩形窗两边突然切断频率响应旁瓣高汉宁窗、汉明窗、布莱克曼窗等把两端压到接近0小来截断突变换来更低的旁瓣代价是主瓣变宽。选择哪个窗要根据任务取舍窗类型主瓣宽度旁瓣衰减典型用途矩形最窄-13dB瞬态测量、整周期采样汉宁较宽-31dB一般频谱分析汉明较宽-43dB语音信号布莱克曼最宽-58dB精度要求高的频谱分析4.4 循环卷积与线性卷积你以为的卷积不是DFT做的那个DFT有一个重要性质时域循环卷积对应频域乘积。注意是“循环卷积”不是我们在滤波器里更常用的“线性卷积”。如果你直接对两个长度都为N的序列做DFT频域乘积再IDFT得到的是循环卷积结果和线性卷积不一样除非补零到足够长度。要求线性卷积时需要把两个序列分别补零到至少N1N2-1的长度再做DFT。这个知识点在我后来的FIR滤波器中用得非常频繁如果你用FFT做快速卷积忘了这步会得到莫名其妙的“环绕”结果。5. 用Python手写DFT/IDFT再和FFT对照5.1 最朴素的O(N²)实现先按定义实现逻辑清晰适合验证理解import numpy as np def my_dft(x): N len(x) X np.zeros(N, dtypecomplex) for k in range(N): for n in range(N): X[k] x[n] * np.exp(-2j * np.pi * k * n / N) return X def my_idft(X): N len(X) x np.zeros(N, dtypecomplex) for n in range(N): for k in range(N): x[n] X[k] * np.exp(2j * np.pi * k * n / N) return x / N测试一下x np.array([1.0, 2.0, 3.0, 4.0]) X my_dft(x) x_hat my_idft(X) print(X) print(x_hat.real)输出结果和前面手算完全一致IDFT后能恢复出[1,2,3,4]虚部是浮点零级的小数。5.2 向量化加速别写双层循环直接用NumPy矩阵运算或者广播机制能快不少def dft_matrix(N): n np.arange(N)[:, None] k np.arange(N)[None, :] return np.exp(-2j * np.pi * k * n / N) def my_dft_fast(x): N len(x) W dft_matrix(N) return W x这个思路在理解DFT是线性变换时也很有用一个N×N矩阵作用于一个N维向量矩阵的每一行就是一个频率模板。5.3 和NumPy官方FFT对照用上面代码和np.fft.fft对比x np.random.randn(64) X1 my_dft(x) X2 np.fft.fft(x) np.testing.assert_allclose(X1, X2, atol1e-10)如果一切正常会静默通过。再试试IDFTx_hat np.fft.ifft(X2) np.testing.assert_allclose(x_hat, x, atol1e-10)实际工程中没人会手写O(N²)的DFT都是用FFT但看懂朴素实现能帮你避开很多“为什么结果和书上不一样”的疑惑。5.4 复杂度对比N1024时的差距朴素的DFT需要N²次复数运算也就是1048576次FFT只需要约N*log2(N)/2大约5120次。N越大差距越恐怖。所以我在嵌入式项目里从来不用库函数里的直接DFT全部走FFT实现哪怕是简单的8点DFT也要写成定点快拍形式。6. 实测经验频率轴标定、幅度归一化与相位展开6.1 频率轴怎么和采样率对应DFT输出的X[k]没有单位你需要自己把k映射到物理频率。如果采样率是fsDFT点数N那么第k个频点对应的物理频率是f_k \frac{k \cdot f_s}{N}这里k的范围是0到N-1。前N/2个点对应0到fs/2的正频率后N/2个点对应从-fs/2到0的负频率或者按周期重复的角度理解。做频谱图时通常只画前N/2或者用np.fft.fftfreq生成真实的频率刻度。我见过至少三个工程师因为直接忽略频率轴刻度把中心频率当成零点导致雷达测速偏差大。建议规范写法fre np.fft.fftfreq(N, d1/fs)6.2 幅度谱的归一化哪些频谱值要乘2要表示真实的单边幅度谱规则是直流分量k0的幅度 |X[0]| / N正频率分量k1到N/2-1的幅度 2 * |X[k]| / N奈奎斯特频率kN/2的特殊情况如果N是偶数且信号只在正频率有实正弦最后一点可能是分量也要单独处理通常也是|X[N/2]|/N很多人直接用np.abs(X)画频谱看到一正一负两个对称峰值还很大就不知所措。其实就是因为没做归一化和单边化处理。以50Hz正弦幅度A1为例采样率1000HzN1000DFT后X[50]的绝对值约为A*N/2500所以除以N乘以2后才是1。6.3 相位谱的展开与IDFT后的小虚部直接计算np.angle(X)得到的是[-π, π]之间的主值但真实相位可能因为延迟线性递增产生跳跃。需要做相位展开phase unwrap才能看到连续的相位变化。NumPy提供np.unwrap能处理。IDFT后理论上应该是纯实数但浮点运算会造成一些NaN级别的虚部残留比如1e-15j。处理办法很简单如果已知原信号是实数就取x_hat.real。不要直接把虚部扔了或者强行取绝对值那会破坏信号。正确做法是设置一个阈值比如小于1e-10就认为虚部是数值噪声。6.4 采样定理与混叠DFT结果“不干净”的根源DFT只处理离散样本可前提是采样率必须大于信号最高频率的两倍否则高频分量会折叠到低频频谱出现虚假分量。这个坑在采集中频信号时尤其严重。我之前做分布式阵列信号处理前端ADC采样率不够结果每个频点都出现一个幻影峰折腾半天才发现是抗混叠滤波器没焊好。如果你看到某个频谱分量在对折频率附近特别异常先检查采样率再检查前端滤波器别急着怀疑DFT算错了。6.5 给自己留个后路使用FFT时的最佳实践最后总结几条我在工程里保留下来的习惯确定采样率和信号有效带宽先估算需要的DFT点数满足频率分辨率的需求再开算对非整周期截断的数据先加窗再做FFT频谱图同时画幅度谱和相位谱不要只看幅度做频域滤波时把频域修改后直接用ifft转回时域注意滤波器的缓降边缘避免吉布斯现象用np.fft.rfft处理实数信号可以省一半内存和运算时间这些习惯帮我少走了很多弯路。信号处理里DFT/IDFT不是背两个公式就完事它像一把刀用好了能切削出清晰的频谱用不好只会得到一堆毛刺。希望这篇能让你更快上手。
返回列表