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

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

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

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

目录

引言

以正弦信号为例

频域表示

绘制DFT生值图

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

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

FFT求功率谱密度

结论

FFT是FastFourierTransform，快速傅里叶变换的缩写。

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

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

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

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

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

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

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

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

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

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

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

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

%产生正弦波形f=10;%正弦波频率，单位Hz|frequencyofsinewaveoverSampRate=30;%采样倍率|oversamplingratefs=overSampRate*f;%采样频率|samplingfrequencyphase=1/3*pi;%初始相位|desiredphaseshiftinradiansnCyl=5;%显示的周期个数|togeneratefivecyclesofsinewavet=0:1/fs:nCyl*1/f;%时间向量|timebasex=sin(2*pi*f*t+phase);%这里也可以替换为cos函数|replacewithcosifacosinewaveisdesiredplot(t,x);title(['正弦波f=',num2str(f),'Hz']);xlabel('时间(s)');ylabel('幅度');

我们在频域上也可以表示给定信号。就是通过快速傅里叶变换(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的数值计算，所以有很多方法来绘制结果。

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

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

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

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

NFFT=1024;%NFFT点DFT|NFFT-pointDFTX=fft(x,NFFT);%使用FFT计算FFT|computeDFTusingFFTnVals=0:NFFT-1;%DFT样本点|DFTSamplepointsplot(nVals,abs(X));title('双边FFT-无移动|DoubleSidedFFT-withoutFFTShift');xlabel('样本点|Samplepoints(N-pointDFT)')ylabel('DFT值|DFTValues');

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

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

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

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

NFFT=1024;%NFFT点DFT|NFFT-pointDFTX=fft(x,NFFT);%使用FFT计算FFT|computeDFTusingFFTnVals=(0:NFFT-1)/NFFT;%归一化DFT样本点|NormalizedDFTSamplepointsplot(nVals,abs(X));title('双边FFT-无移动|DoubleSidedFFT-withoutFFTShift');xlabel('归一化频率|NormalizedFrequency')ylabel('DFT幅度值|DFTValues');

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

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

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

从-0.5到0.5的x轴。

fftshift-Shiftzero-frequencycomponenttocenterofspectrum

ThisMATLABfunctionrearrangesaFouriertransformXbyshiftingthe

zero-frequencycomponenttothecenterofthearray.

Y=fftshift(X)

Y=fftshift(X,dim)

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

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

figureNFFT=1024;%%NFFT点DFT|NFFT-pointDFTX=fftshift(fft(x,NFFT));%使用FFT计算FFT|computeDFTusingFFTfVals=(-NFFT/2:NFFT/2-1)/NFFT;%归一化DFT样本点|DFTSamplepointsplot(fVals,abs(X));title('双边FFT|DoubleSidedFFT-withFFTShift');xlabel('归一化频率|NormalizedFrequency')ylabel('DFT幅度值|DFTValues');

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

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

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

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

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

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

图5将横坐标转换为频率

figureNFFT=1024;X=fftshift(fft(x,NFFT));fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;plot(fVals,abs(X),'b');title('DoubleSidedFFT-withFFTShift');xlabel('Frequency(Hz)')ylabel('|DFTValues|');

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

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

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

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

图6功率谱密度

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

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

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

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

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

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

图8单边PSD

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

[1]Mathuranathan,"HowtoplotFFTusingMatlab–FFTofbasicsignals:SineandCosinewaves",July16,2014.

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

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

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

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