MATLAB综合实例：时域信号的频谱分析

---

前言

第二天，小刘吃完午饭正趴在桌子上睡觉，做着暴富后，把苦茶子老板辞退的美梦。就在这时，老板从办公室走了出来，敲了敲小刘的桌子。小刘揉了揉眼睛，起来刚要起难，发现是他苦茶子老板，赶紧闭上了要张开的嘴，谄笑的：“老板有何吩咐啊？”
苦茶子老板，拍了拍小刘的肩膀说道：“小刘辛苦了啊！听说你是电子信息工程专业的？”
“是啊，是啊。”
苦茶子老板笑着说道：“这不是巧了吗？我儿子也是电子信息工程专业的，他跟我说他们信号与系统要做什么时域信号的频谱分析。你说我一个搞事业的，哪里懂那些嘛！这个就交给你啦！要求已经发你手机上了，一个小时后发给我啊！”
小刘打开手机【假设数据采样频率为1000Hz，一个信号包含两个正弦波，频率为50Hz、120Hz，振幅为0.7、1，噪声为零平均值的随机噪声。试采用傅里叶变换分析其频谱】。
接着小刘打开了CSDN搜索到：“如何用MATLAB制作时域信号的频谱分析”

---

一、快速傅里叶变换

快速傅里叶变换（FFT）是离散傅里叶变换的快速算法，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅里叶变换的算法进行改进获得的。

快速傅里叶变换命令的调用格式：

命令	意义	命令调用格式
fft	一维快速傅里叶变换	1.Y=fft(X)，计算对向量X的快速傅里叶变换。如果X是矩阵，fft返回对每一列的快速傅里叶变换；2.Y＝fft（X，n），计算向量X的n点FFT。当X的长度小于n时，系统将在X的尾部补零，以构成n点数据；当X的长度大于n时，系统进行截尾；3.Y＝fft（X，［］，dim）或Y＝fft（X，n，dim），计算对指定的第dim维的快速傅里叶。
ff2	二维快速傅里叶变换	1.Y＝fft2（X），计算对X的二维快速傅里叶变换。结果Y与X的维数相同；2.Y＝fft2（X，m，n），计算结果为mxn阶，系统将视情况对X进行截尾或者以0来补齐。
fftshift	将快速傅里叶变换（fft、fft2）的DC分量转移至谱中心	1.Y＝fftshift（X），将DC分量转移至谱中心；2.Y＝fftshift（X，dim），将DC分量转移至dim维谱中心，若dim为1，则上下转移；若dim为2，则左右转移。
ifft	一维逆快速傅里叶变换	1.y＝ifft（X），计算X的逆快速傅里叶变换；2.y＝ifft（X，n），计算向量X的n点逆FFT；3.y＝ifft（···，symflag），计算对指定Y的对称性的逆FFT；4.y＝ifft（X，n，dim），计算对dim维的逆FFT
ifft2	二维逆快速傅里叶变换	1.y＝ifft2（X），计算X的二维逆快速傅里叶变换；2.y＝ifft2（X，m，n），计算向量X的mxn维逆快速傅里叶变换。
ifftn	多维逆快速傅里叶变换	1.y＝ifftn（X），计算X的n维逆快速傅里叶变换；2.y＝ifftn（X，size），系统将视情况对X进行截尾或者以0来补齐。
ifftshift	逆fft平移	1.Y＝ifftshift（X），同时转移行与列；2.Y＝ifftshift（X，dim），若dim为1，则行转移；若dim为2，则列转移

利用快速傅里叶变换实现快速卷积。

>> A=magic(4);%生成4阶魔方矩阵
>> B=ones(3);%生成3阶全1矩阵
>> A(6,6)=0;%将A用零补全为（4+3-1）*（4+3-1）维
>> B(6,6)=0;%将B用零补全为（4+3-1）*（4+3-1）维
>> C=ifft2(fft2(A).*fft2(B));%对A、B进行二维快速傅里叶变换，并将结果相乘，对乘积。
>> %括号内进行二维逆快速傅里叶变换，得到卷积。
>> C

C =

  1 至 5 列

      16             18             21             18             16       
      21             34             47             47             34       
      30             50             69             72             52       
      18             50             81             84             52       
      13             34             55             55             34       
       4             18             33             30             16       

  6 列

      13       
      21       
      33       
      21       
      13       
       1       

>> %进行演算
>> A=magic(4);
>> B=ones(3);
>> D=conv2(A,B)

D =

  1 至 5 列

      16             18             21             18             16       
      21             34             47             47             34       
      30             50             69             72             52       
      18             50             81             84             52       
      13             34             55             55             34       
       4             18             33             30             16       

  6 列

      13       
      21       
      33       
      21       
      13       
       1       

二、时域信号频谱分解操作

1.零平均值的随机噪声

代码如下

>> Fs=1000;
>> T=1/Fs;
>> L=1000;%信号长度
>> t=(0:L-1)*T;%时间向量
>> x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t);%正弦信号的表达式。
>> y=x+2*randn(size(t));%噪声正弦信号和x正弦信号。
>> plot(Fs*t(1:50),y(1:50))%绘制随机噪声信号波
>> title('零平均值噪声信号');
>> xlabel('time(milliseconds)')%标注x轴

2.单边振幅频谱

代码如下

>> NFFT=2^nextpow2(L);%传递给fft的信号长度
>> Y=fft(y,NFFT)/L;%对信号进行快速傅里叶变换，将时域信号转变为频谱
>> f=Fs/2*linspace(0,1,NFFT/2);%FFT变换后的频率
>> plot(f,2*abs(Y(1:NFFT/2))，‘r’)%绘制单边振幅频谱
>> title('y(t)单边振幅频谱')
>> xlabel('Frequency(Hz)')
>> ylabel('|Y(f)|')

---

总结

小刘将做出来的时域信号的频谱分析打印出来发送给了老板，老板看后表示，不愧是电子信息工程专业的高材生呀！随后将文件发送给了自己的儿子，表示要给小刘安排更多的活，让他加速成长！