频谱分析？FFT函数正确使用从零开始

大家好，欢迎来到IT知识分享网。

班长聊通信原理 | 系列文章关注”班长”→主页→”文章”查看

肺炎疫情让年味儿多了不少焦虑，走亲访友都成了奢望。与其冒着被感染的风险出门或聚会，倒不如宅在家里最踏实。

热爱学习的同学们，宅在家里，学点通信原理吧。

本文主要介绍MATLAB中FFT函数的使用，对于FFT|DFT的原理未有提及，需要的同学去文末自取。

目录

1. 引言
2. 以正弦信号为例
3. 频域表示
4. 绘制DFT生值图
5. 根据归一化频率轴绘制原始值
6. 归一化频率轴、调整为正负频率
7. FFT求功率谱密度
8. 结论

01引言

FFT是Fast Fourier Transform，快速傅里叶变换的缩写。

傅里叶变换可以让我们，由时域分析转入变换域(频域)分析。

毫无疑问，我们的实验对象是常见的标准信号。

教材中使用的标准信号有：正弦、余弦、高斯脉冲、方波、孤立矩形脉冲、指数衰减、Chirp信号。

在频域分析之前，我们需要思考如何在MATLAB中产生这些标准信号。

02以正弦信号为例

为了在MATLAB中产生正弦波，第一步是固定正弦波的频率f。

例如，我打算生成一个f = 10Hz的正弦波，其最小振幅和最大振幅分别为-1v和1v。

因为MATLAB是一种处理数字比特的软件，所以既然已经确定了正弦波的频率，下一步就是确定采样率。

根据Nyquist Shannon定理，为了产生/绘制平滑的正弦波，采样率必须远远高于规定的最小所需采样率，该采样率至少是频率f的两倍。

这里选择了30的过采样倍数，这是为了绘制一个平滑的连续的正弦波。

因此，采样率为f_s=30倍f=300Hz。

如果正弦波需要相移，那就单独定义一个初相位。

图1 正弦信号，5个周期，频率为10Hz

% 产生正弦波形 f = 10; %正弦波频率，单位Hz | frequency of sine wave overSampRate = 30; %采样倍率 | oversampling rate fs = overSampRate*f; %采样频率 | sampling frequency phase = 1/3*pi; %初始相位 | desired phase shift in radians nCyl = 5; %显示的周期个数 | to generate five cycles of sine wave t = 0:1/fs:nCyl*1/f; %时间向量 | time base x = sin(2*pi*f*t+phase); %这里也可以替换为cos函数 | replace with cos if a cosine wave is desired plot(t,x); title(['正弦波 f=', num2str(f), 'Hz']); xlabel('时间(s)'); ylabel('幅度');

03频域表示

我们在频域上也可以表示给定信号。就是通过快速傅里叶变换(FFT)来实现的。

MATLAB中FFT(x,n)可以计算n点的DFT。

一般在DFT计算中的点数n取为2的幂，以方便FFT的有效计算。

这里选择n=1024的值。

fft – 快速傅里叶变换

此 MATLAB 函数 用快速傅里叶变换 (fft) 算法计算 X 的离散傅里叶变换 (DFT)。 如果 X 是向量，则 fft(X) 返回该向量的傅里叶变换。

如果 X 是矩阵，则 fft(X) 将 X 的各列视为向量，并返回每列的傅里叶变换。 如果 X 是一个多维数组，则 fft(X) 将沿大小不等于 1

的第一个数组维度的值视为向量，并返回每个向量的傅里叶变换。

Y = fft(X)

Y = fft(X,n)

Y = fft(X,n,dim)

由于FFT只是n点DFT的数值计算，所以有很多方法来绘制结果。

04绘制DFT生值图

所谓的生值图2，就是DFT的原始值：x轴从0到n-1，表示n个样本值。

由于DFT值是复数，所以DFT abs(X)的大小被绘制在y轴上。

从这个图2中，频率轴(X轴)就是简单的0到1023，共计1024个点。看起来和”频率”没啥关系，就是1024个点计算DFT。

图2 FFT函数输出值，直接画图，注意横坐标为样本点

NFFT=1024; % NFFT点DFT | NFFT-point DFT X=fft(x,NFFT); %使用FFT计算FFT | compute DFT using FFT nVals=0:NFFT-1; %DFT样本点 | DFT Sample points plot(nVals,abs(X)); title('双边FFT-无移动 | Double Sided FFT - without FFTShift'); xlabel('样本点 | Sample points (N-point DFT)') ylabel('DFT值 | DFT Values');

05根据归一化频率轴绘制原始值

将频率轴(x轴)归一化。

只需将x轴上的样本索引除以FFT的长度n，这使x轴相对于采样率f_s规范化。

然而，我们仍然不能从图3中找出正弦的频率。

图3 把频率归一化，就是除以样本点总数

NFFT=1024; %NFFT点DFT | NFFT-point DFT X=fft(x,NFFT); %使用FFT计算FFT | compute DFT using FFT nVals=(0:NFFT-1)/NFFT; %归一化DFT样本点 | Normalized DFT Sample points plot(nVals,abs(X)); title('双边FFT-无移动 |Double Sided FFT - without FFTShift'); xlabel('归一化频率 | Normalized Frequency') ylabel('DFT幅度值 | DFT Values');

06归一化频率轴、调整为正负频率

在频域中，由于引入了复指数函数，所以会有正、负频率轴。

所以，为了在具有正值和负值的频率轴上绘制DFT值，样本索引0处的DFT值必须作为数组中心，然后“对称”到负频率轴。

这是通过在MATLAB中使用fftshift函数来完成的。

从-0.5到0.5的x轴。

fftshift – Shift zero-frequency component to center of spectrum

This MATLAB function rearranges a Fourier transform X by shifting the

zero-frequency component to the center of the array.

Y = fftshift(X)

Y = fftshift(X,dim)

fftshiftt通过将零频率分量移动到阵列的中心，这个matlab函数重新排列一个傅里叶变换X。

图4 fftshiftt通过将零频率分量移动到阵列的中心

figure NFFT=1024; %%NFFT点DFT | NFFT-point DFT X=fftshift(fft(x,NFFT)); %使用FFT计算FFT | compute DFT using FFT fVals=(-NFFT/2:NFFT/2-1)/NFFT; %归一化DFT样本点 | DFT Sample points plot(fVals,abs(X)); title('双边FFT|Double Sided FFT - with FFTShift'); xlabel('归一化频率 | Normalized Frequency') ylabel('DFT幅度值 | DFT Values');

至此，横坐标都不能算“真正”的频率，而我们恰恰想要的是频率图。

所以，我们用fs/NFFT可以得到频率间隔△f。

从图5中，我们可以确定FFT绝对值的峰值在10Hz和-10Hz处。

因此，产生的正弦波的频率为10Hz。

这是我们只看频谱图5得出的结论，和我们仿真的一致。

在10Hz和-10Hz峰值旁边的小旁瓣是由于谱泄漏造成的(后续会说这个问题)

图5 将横坐标转换为频率

figure NFFT=1024; X=fftshift(fft(x,NFFT)); fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT; plot(fVals,abs(X),'b'); title('Double Sided FFT - with FFTShift'); xlabel('Frequency (Hz)') ylabel('|DFT Values|');

07FFT求功率谱密度

图6绘制了每个频率分量的功率。功率可以在线性标度或对数标度中绘制。

每个频率分量的能量计算为

P_x(f) = x(f) * x^*(f)

其中x(f)是信号x(t)的频域表示，x^*(f)是其共轭。

图6 功率谱密度

% 功率谱密度，频率移动，幅度值的模 figure NFFT=1024; L=length(x); X=fftshift(fft(x,NFFT)); Px=X.*conj(X)/(NFFT*L); %计算每个频率分量的功率 | Power of each freq components fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT; plot(fVals,Px,'r'); title('功率谱密度 | Power Spectral Density'); xlabel('频率 | Frequency (Hz)') ylabel('功率 | Power');

在对数log尺度上绘制PSD图，为信号处理中产生最常见的PSD图7。

% 功率谱密度，LOG表示，频率移动，幅度值的模 figure NFFT=1024; L=length(x); X=fftshift(fft(x,NFFT)); Px=X.*conj(X)/(NFFT*L); %计算每个频率分量的功率| Power of each freq components fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT; plot(fVals,10*log10(Px),'b'); title('功率谱密度 | Power Spectral Density'); xlabel('频率 | Frequency (Hz)') ylabel('功率 | Power');

图7 PSD功率谱密度，对数坐标

在图8中，省略了x轴的负频部分，只绘制了对应于n点DFT的0到n/2个样本点的FFT值。

相应地，归一化频率轴在0到0.5之间运行，绝对频率(x轴)从0到f_s/2。

图8 单边PSD

结论

关于FFT的使用，还有很多有趣的问题：比如谱泄露spectral leakage，采样点对FFT结果的影响等等。班长会在后续过程中与大家讨论，敬请期待。

相关文章参考

[1]Mathuranathan,”How to plot FFT using Matlab – FFT of basic signals : Sine and Cosine waves”,July 16, 2014.

[2]库利-图基FFT算法，4G通信系统中使用了FFT算法进行时域频域切换

[3]OFDM技术：信号的产生为何与FFT算法有关？为什么要串并转换？

[4]想要画出正确的频谱图，不是直接调用MATLAB FFT函数那么简单

[5]傅里叶变换，是如何让我们抓狂的？

看到这里，欢迎为文章点赞，期待您在留言区的精彩评论！