模拟信号的数字谱分析

1

内容与要求

画出不同采样率、不同持续期的正弦信号波形图，分别用DTFT 、DFT 和FFT 算法计算频谱。通过MATLAB 仿真计算，达到以下目标：

1. 理解数字谱分析中物理分辨力和计算分辨力的概念及两者间的关系；

2. 弄清数字谱分析中误差成因及减小误差的方法。 2

报告正文

（1）数字谱分析中物理分辨力和计算分辨力的概念及两者间的关系： 物理分辨率△fw=c*fs/L 计算分辨率△fc=fs/N

当c ，fs 都不变时，增大N 可以显著提高计算分辨率，但物理分辨率不变，因为有效的数据点数不变；增大L 时，可以提高物理分辨率，能分辨出间隔少的频率，但计算分辨率不变。窗函数频谱的主瓣应尽量窄,能量尽可能集中在主瓣内,从而在谱分析时获得较高的频率分辨力 。应用DFT 进行信号的数字谱分析时，需要确定的参数主要有：截取的信号长度T 、抽样频率fs 、采样点数N 及选择什么样的窗口。如果在x(n)中有两个频率分别为f1和f2的信号，对x(n)用矩形窗截断时，要分辨出这两个频率，N 必须满足[2]：

212fs

f f N

<- 因此，反映频率分辨率能力的Δf 由下式决定：Δf = fs/ N = 1/( N*T)= 1/L 可见，Δf 反比于截取的信号长度L ，即截取的信号长度越长，分辨能力越好 下面用例子证明上面所分析的结论。设有一序列x(n)含有两种频率成分，f1= 3Hz ，f2=3.04Hz ，()cos(21/2)cos(22/)x n nf f nf fs ππ=+。由抽样定理，我们采用fs =10Hz ，即 T= 0.1s 对其抽样，不应发生频谱混叠问题。这里采用128点和512 点DFT ，对两种情况下的频率分辨率进行比较：设128 点的分辨率为Df1 ， 512 点的分辨率为Df2,则df1=10/128=0.078125Hz,Df2=10/512= 0.01953125Hz 。计算结果表明，128点DFT 未能将两个靠的很近的频率分辨出来，而512点DFT 却达到了分辨频率的要求。

（2）数字谱分析中误差成因及减小误差的方法： ①频谱泄漏：

泄漏现象的产生主要是由于窗函数的突然截断信号在谱函数中产生旁瓣而引起的，谱函数中主瓣的幅值与旁瓣幅度的比值越大则泄漏现象对信号频谱分析的影响就越小。而离散时间序列通常是无限长序列，为了能计算和物理实现，通常都会截取一段来进行计算，也就是加窗。但是加窗后的频谱会发生泄漏，就像频率大小为fs 的冲激经过加窗后会变成以fs 为中心的频带，这就会造成频谱泄漏。 频谱泄漏的危害：

1. 频谱有皱纹，其中较大的旁瓣可能淹没较小谱峰。

2. 产生新的谱峰，造成观察者对信号的误读。（物理分辨率下降） 解决方法：可以通过补零加长窗的长度来减少频谱的泄漏。 ②栅栏效应：

对一函数实行采样就像是通过栅栏看外面一样，就是用栅栏看到的景物来代替全部，这就是栅栏效应。但是这就可能略去重要的部分，使采样失去意义。所以信号的采样要满足采样定理：N>L 。

解决方法：可以通过提高采样点数来减少栅栏效应。 ③混淆：

采样的过程不符合采样定理，即不符合max 2s f f ≥。这就会造成由于采样频率过低而产生的混淆从而不能无失真的恢复。 解决方法：令max 2s f f ≥。 (3)程序： 程序1：

clf

clear all;

L=80;%---采样长度 n=0:L-1;%---采样间隔 N=512;%---采样点数 fs=6000;%---采样频率 f1=100; %---模拟频率

w1=2*pi*f1/fs;%---数字角频率 t=0:0.0001:L;

xt=sin(w1*t);%---原始正弦信号

subplot(5,1,1)

plot(t,xt)

title('正弦信号')

xn=sin(w1*n);%---采样后的信号

subplot(5,1,2)

stem(n,xn,'k.')

title(['采样信号fs=',num2str(fs)]);

%---DTFT

k=0:N-1;

w = (0:N-1) * 2 * pi / N;

Y1=DTFT(xn,w);

subplot(5,1,3)

plot(w/(2*pi),abs(Y1),'r')

title(['DTFT L=',num2str(L),',N=',num2str(N),',fs=',num2str(fs)]) %---DFT

Y2=DFT(xn,N);

subplot(5,1,4)

stem(k*fs/N,abs(Y2),'b.')

title(['DFT L=',num2str(L),',N=',num2str(N),',fs=',num2str(fs)]) %---FFT

Y3=fft(xn,N);

subplot(5,1,5)

stem(k*fs/N,abs(Y3),'g*')

title(['FFT L=',num2str(L),',N=',num2str(N),',fs=',num2str(fs)])

程序2：

%--DTFT函数

function X=DTFT(x,w)

L=length(x);

z=exp(-1j*w);

X=0;

for n=L-1:-1:0

X=x(n+1)+z.*X;

End

程序3：

%---DFT函数

function X=DFT(x,N)

% N为采样点数

L=length(x);%时间序列长度

for k=0:N-1

n=0:L-1;

omega_k=2*pi*k/N;