% 产生正弦波形f = 10; %正弦波频率，单位Hz | frequency of sine waveoverSampRate = 30; %采样倍率 | oversampling ratefs = overSampRatef; %采样频率 | sampling frequencyphase = 1/3pi; %初始相位 | desired phase shift in radiansnCyl = 5; %显示的周期个数 | to generate five cycles of sine wavet = 0:1/fs:nCyl1/f; %时间向量 | time basex = sin(2pift+phase); %这里也可以替换为cos函数 | replace with cos if a cosine wave is desiredplot(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通过将零频率分量移动到阵列的中心figureNFFT=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 将横坐标转换为频率figureNFFT=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 功率谱密度% 功率谱密度，频率移动，幅度值的模figureNFFT=1024;L=length(x); X=fftshift(fft(x,NFFT)); Px=X.conj(X)/(NFFTL); %计算每个频率分量的功率 | 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表示，频率移动，幅度值的模figureNFFT=1024; L=length(x); X=fftshift(fft(x,NFFT)); Px=X.conj(X)/(NFFTL); %计算每个频率分量的功率| Power of each freq components fVals=fs(-NFFT/2:NFFT/2-1)/NFFT; plot(fVals,10log10(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]傅里叶变换，是如何让我们抓狂的？看到这里，欢迎为文章点赞，期待您在留言区的精彩评论
（图片来源网络，侵删）
0 评论