信号的频谱分析及MATLAB 实现（实例）
摘自：张登奇,杨慧银.信号的频谱分析及MATLAB 实现[J].湖南理工学院学报(自然科学版),2010,(03)
摘 要：DFT 是在时域和频域上都已离散的傅里叶变换，适于数值计算且有快速算法，是利用计算机实现信号频谱分析的常用数学工具。
文章介绍了利用DFT 分析信号频谱的基本流程，重点阐述了频谱分析过程中误差形成的原因及减小分析误差的主要措施，实例列举了MATLAB 环境下频谱分析的实现程序。
通过与理论分析的对比，解释了利用DFT 分析信号频谱时存在的频谱混叠、频谱泄漏及栅栏效应，并提出了相应的改进方法。
关键词：MATLAB ；频谱分析；离散傅里叶变换；频谱混叠；频谱泄漏；栅栏效应
3 分析实例
对信号进行频谱分析时，由于信号不同，傅里叶分析的频率单位也可能不同，频率轴有不同的定标方式。
为了便于对不同信号的傅里叶分析进行对比，这里统一采用无纲量的归一化频率单位，即模拟频率对采样频率归一化；模拟角频率对采样角频率归一化；数字频率对2π归一化；DFT 的k 值对总点数归一化。
同时，为了便于与理论值进行对比，理解误差的形成和大小，这里以确定信号的幅度谱分析为例进行分析说明。
假设信号为：)()(t u e t x t
-=，分析过程：首先利用CTFT 公式计算其模拟频谱的理论值；然后对其进行等间隔理想采样，得到)(n x 序列，利用DTFT 公式计算采样序列的数字连续频谱理论值，通过与模拟频谱的理论值对比，理解混叠误差形成的原因及减小误差的措施；接下来是对)(n x 序列进行加窗处理，得到有限长加窗序列)(n xw ，再次利用DTFT 公式计算加窗后序列)(n xw 的数字连续频谱，并与加窗前)(n x 的数字连续频谱进行对比，理解截断误差形成的原因及减小误差的措施；最后是对加窗序列进行DFT 运算，得到加窗后序列)(n xw 的DFT 值，它是对)(n xw 数字连续频谱进行等间隔采样的采样值，通过对比，理解栅栏效应及DFT 点数对栅栏效应的影响。
利用MATLAB 实现上述分析过程的程序如下： clc;close all;clear;
%CTFT 程序，以x(t)=exp(-t) t>=0 为例
%利用数值运算计算并绘制连续信号波形
L=4, %定义信号波形显示时间长度
fs=4,T=1/fs; %定义采样频率和采样周期
t_num=linspace(0,L,100);%取若干时点,点数决定作图精度
xt_num=exp(-1*t_num);%计算信号在各时点的数值
subplot(3,2,1);plot(t_num,xt_num),%绘信号波形
xlabel('时间(秒)'),ylabel('x(t)'),%加标签
grid,title('(a) 信号时域波形'),%加网格和标题
%利用符号运算和数值运算计算连续信号幅度谱的理论值
syms t W %定义时间和角频率符号对象
xt=exp(-1*t)*heaviside(t),%连续信号解析式
XW=fourier(xt,t,W),%用完整调用格式计算其傅氏变换
%在0两边取若干归一化频点，点数决定作图精度
w1=[linspace(-0.5,0,50),linspace(0,1.5,150)];
XW_num=subs(XW,W,w1*2*pi*fs);%利用置换函数求频谱数值解mag1=abs(XW_num);%计算各频点频谱的幅值
subplot(3,2,2);plot(w1,mag1),%绘制归一化频率幅值谱线xlabel('频率(*2*pi*fs)rad/s'),ylabel('幅度'),%加标签grid,title('(b) 连续信号幅频理论值'),%加网格和标题
%DTFT程序，以x(n)=exp(-nT) n>=0 为例
%利用数值运算计算并绘制离散信号图形
N=L*fs+1;n_num=0:N-1;%生成信号波形采样点
xn_num=exp(-1*n_num*T);%计算信号理想采样后的序列值subplot(3,2,3);stem(n_num,xn_num,'b.'),%绘序列图形xlabel('n'),ylabel('x(n)'),%加标签
grid,title('(c) 理想采样图形'),%加网格和标题
%利用符号运算和数值运算计算离散信号幅度谱的理论值syms n z w %定义符号对象
xn=exp(-n*T), %定义离散信号
Xz=ztrans(xn,n,z),%用完整调用格式计算其Z变换
%利用复合函数计算序列傅里叶变换的解析解
Z=exp(j*w);Hejw=compose(Xz,Z,z,w);
Hejw_num=subs(Hejw,w,w1*2*pi);%求频谱数值解
mag2=abs(Hejw_num);%计算各频点频谱的幅度
subplot(3,2,4);plot(w1,mag2*T),%绘制频谱幅度曲线xlabel('频率(*2*pi)rad'),ylabel('幅度'),%加标签
grid,title('(d) 离散信号幅频理论值'),%加网格和标题
%序列加窗图示及频谱幅值绘制
%利用数值运算计算并绘制加窗后序列xw(n)的图形
M=8;win=(window(@rectwin,M))';%定义窗点和窗型
xwn=xn_num.*[win,zeros(1,N-M)];%给离散信号加窗
subplot(3,2,5);stem(n_num,xwn,'b.'),%加窗序列图示xlabel('n'),ylabel('xw(n)'),%加标签
grid,title('(e) 加窗序列图形'),%加网格和标题
%利用符号运算和数值运算计算加窗序列的频谱幅值
%先求加窗序列的Z变换,注意表达式长度限制问题
Xwz=0;for n=0:(M-1);Xwz = Xwz+xwn(n+1)*z^(-n); end
%利用复合函数计算加窗序列傅里叶变换的解析解
Zw=exp(j*w);HejwM=compose(Xwz,Zw,z,w);
HejwM_num=subs(HejwM,w,w1*2*pi);%求频谱数值解
mag3=abs(HejwM_num);%计算各频点频谱的幅度
subplot(3,2,6);plot(w1,mag3*T),%绘频谱幅度曲线
%利用DFT计算加窗序列xw(n)的离散谱幅值
Ndft=16, Xk=fft(xwn,Ndft);%定义DFT点数和DFT运算
Xk0=fftshift(Xk)*T;%将DFT值0对称和幅值加权处理
if mod(Ndft,2)==0; N1=Ndft; else N1=Ndft-1; end;
k=[0:(Ndft-1)]-N1/2;wk=k/Ndft;%0对称取值并归一化
hold on;stem(wk,abs(Xk0),'r.'),%绘制DFT图形
legend('幅谱','DFT',0),%加响应图例，位置自动最佳
xlabel('归一化频率'),ylabel('幅度'),%加标签
grid,title('( f ) 加窗序列幅谱及其DFT幅值'),
plot(w1,mag1,'k:'),%与连续信号幅谱的理论值比较
该程序过程清晰、容易理
解，程序运行结果如图2所示。
图2（a）是信号的时域波形，
图2（b）是对应的幅度谱图。
由于在归一化频率为0.5的地
方还有较大幅度，所以对信号进
行理想采样后存在较大的混叠
失真，表现在图2（d）中归一
化频率为-0.5～0.5范围内的波
形与图2（b）的波形明显不同。
图2（e）是理想采样序列加矩
形窗得到的图形，对应的幅度谱
线如图2（f）中实线所示。
该
谱线与图2（d）相比有明显的
不同（如波动现象），这是加窗
图2 频谱分析图解说明
截断的结果。
窗口长度越短截断
效应会越明显。
对加窗序列进行DFT运算，只要DFT点数大于等于窗口长度，算出的DFT值就是对加窗序列的连续频谱在一个周期内进行的等间隔采样的采样值。
在图2（f）中表现为DFT幅值在加窗序列连续幅谱的谱线上，是连续频谱曲线上的有限个数据点，即所谓栅栏效应。
图2（f）中用虚线画出了连续信号的幅谱理论值，与DFT的幅值对比，存在一定误差，但只要采样频率足够高，时窗长度足够长，FFT 点数足够大，得到的DFT值越逼近实际频谱。