DSP题目：FFT算法的Matlab实现及其应用研究

DSP 题目FFT算法的Matlab实现及应用研究最近帮室友调毕设的信号处理部分他拿了个麦克风录的杂音想把背景的50Hz工频噪音去掉上来就问我“为啥我fft出来的峰不对”——害这问题我刚学DSP的时候也踩过无数坑今天就来唠唠FFT的Matlab实现和那些容易踩的坑。先从最直观的例子开始分析正弦信号频谱先搞个最简单的场景生成一个100Hz的纯净正弦波用Matlab快速算它的频谱% 基础参数设置 Fs 1000; % 采样频率要至少是信号最高频率的2倍这里留了余量 T 1/Fs; % 单次采样间隔 L 1000; % 总采样点数对应1秒时长 t (0:L-1)*T; % 生成时间轴 % 生成100Hz正弦信号 x sin(2*pi*100*t); % 计算FFT并处理频谱 Y fft(x); % 先算双侧频谱再转成我们常用的单侧频谱 P2 abs(Y/L); P1 P2(1:L/21); % 除了直流和奈奎斯特频率点其他点都要翻倍因为能量被分到正负频率两侧了 P1(2:end-1) 2*P1(2:end-1); % 生成对应的频率轴 f Fs*(0:(L/2))/L; % 画图看结果 figure(Name,纯净正弦信号频谱) plot(f,P1) title(100Hz正弦波单侧频谱) xlabel(频率(Hz)) ylabel(幅值) grid on这里很多新手会卡壳为啥算出来的峰值是0.5而不是1就是没做那个归一化和单侧转换——Matlab内置的fft输出的是包含正负频率的完整频谱直接取绝对值的话幅值会翻倍而且我们只需要看0到Fs/2的部分就够了。实战场景带噪信号滤波就说室友遇到的问题他录的信号里混了50Hz工频噪和高斯白噪我们来模拟一下% 生成带噪信号100Hz信号 50Hz工频噪 随机白噪 x_noisy sin(2*pi*100*t) 0.3*sin(2*pi*50*t) 0.2*randn(size(t)); % 同样做频谱分析 Y_noisy fft(x_noisy); P2_noisy abs(Y_noisy/L); P1_noisy P2_noisy(1:L/21); P1_noisy(2:end-1) 2*P1_noisy(2:end-1); figure(Name,带噪信号频谱) plot(f,P1_noisy) title(带噪信号单侧频谱) xlabel(频率(Hz)) ylabel(幅值) grid on这时候能清晰看到50Hz和100Hz的两个峰还有一堆杂乱的小尖峰——就是高斯白噪的频谱。接下来我们把50Hz的峰删掉% 滤除50Hz噪音 Y_filtered Y_noisy; % 找到50Hz对应的索引区间留一点余量防止频偏 idx 48:52; % 实信号的FFT是共轭对称的正负频率位置都要置0 Y_filtered(idx) 0; Y_filtered(L - idx 2) 0; % 逆FFT恢复干净信号 x_filtered ifft(Y_filtered); % 画图对比效果 figure(Name滤波前后信号对比) subplot(2,1,1) plot(t, real(x_noisy)) title(原始带噪信号) xlabel(时间(s)) subplot(2,1,2) plot(t, real(x_filtered)) title(滤除50Hz噪后的信号) xlabel(时间(s))这里要注意ifft输出的是复数结果因为我们的原始信号都是实信号所以只需要取real部分就行虚部都是浮点运算的微小误差可以直接忽略。自己写一个极简版FFT搞懂原理光用内置函数总觉得隔着一层我们来写个最基础的基2FFT也就是课本里的Cooley-Tukey算法原理就是把N点DFT拆成两个N/2点的DFT再用蝶形运算拼起来复杂度从O(N²)降到O(NlogN)function y my_fft(x) N length(x); % 递归终止条件单点DFT就是它本身 if N 1 y x; return; end % 奇偶项拆分 x_even x(1:2:end); x_odd x(2:2:end); % 递归计算两个子DFT y_even my_fft(x_even); y_odd my_fft(x_odd); % 蝶形运算组合结果 y zeros(1, N); for k 1:N/2 % 旋转因子 twiddle exp(-2*pi*1j*(k-1)/N); y(k) y_even(k) twiddle * y_odd(k); y(k N/2) y_even(k) - twiddle * y_odd(k); end end我们来测试一下和内置fft的差距% 测试自定义FFT x_test sin(2*pi*100*t(1:64)); % 取64个点刚好是2的整数次幂 Y_mine my_fft(x_test); Y_builtin fft(x_test); % 计算最大误差基本在1e-10级别都是浮点精度问题 err max(abs(Y_mine - Y_builtin)); disp([自定义FFT和内置FFT的最大误差, num2str(err)]);跑出来的误差极小说明我们自己写的函数是对的就是递归版本效率不高实际工程里还是用内置的fft更快。避坑指南几个新手最容易踩的点补零≠提升分辨率很多人以为补零能让频谱更精细但其实真实的频谱分辨率是Fs/LL是原始采样点数补零只是在频域做插值让图像看起来更平滑而已。比如我们只采了100个点分辨率就是10Hz两个间隔小于10Hz的信号峰是没法分辨的补零到1024点也没用。matlab% 补零效果演示xshort sin(2pi100t(1:100)); % 仅100个采样点Lshort length(xshort);% 不补零的频谱Y1 fft(xshort);P11 2abs(Y1(1:Lshort/21))/Lshort;f1 Fs(0:(Lshort/2))/Lshort;% 补零到1024点的频谱Y2 fft(xshort, 1024);P12 2abs(Y2(1:1024/21))/Lshort;f2 Fs*(0:(1024/2))/1024;figure(Name,补零对比)subplot(2,1,1)plot(f1,P1_1)title(未补零频谱)xlim([90 110])DSP 题目FFT算法的Matlab实现及应用研究subplot(2,1,2)plot(f2,P1_2)title(补零到1024点频谱)xlim([90 110])别忘了去直流分量如果你的信号有直流偏移频谱在0Hz的地方会出现一个很高的尖峰不需要的话可以先把信号减去均值x x - mean(x);频率轴别算错别写成f (0:L-1)*Fs/L;然后直接用单侧频谱只需要取前半部分而且要对应好索引。最后说一句FFT的应用远不止滤波图像处理里的二维FFT、音频压缩、雷达信号分析都离不开它先把基础的坑踩完后面玩起来就顺多了。