内容与要求
利用DFT对多种信号(例如由多个正弦信号组成的信号)进行频谱分析,并研究不同数据长度、补零、加窗等对频率分辨率的影响。
方法原理
1、引入
当数字计算机对信号进行频谱分析时,要求信号必须以离散值作为输入,而计算机输出所得的频谱值自然也是离散的。因此,要使信号是时间的连续函数、频谱是频率的连续函数或者信号及频谱二者都是变量的连续函数这三种形式的信号能用数字计算机进行计算,必须针对每一种形式的具体情况,或者在时域与频域上取样,或者在时域上取样,或者在频域上取样。信号在时域上取样导致频率的周期函数,在频域上取样导致时域的周期函数,最后都将使原时间函数和频率函数二者都成为周期离散的函数。我们采用DFT(离散傅里叶变换)来对连续时间信号的傅里叶变换进行逼近,进而分析连续时间信号的频谱。
离散傅里叶变换是有限长序列的傅里叶变换,它相当于把信号的傅里叶变换进行等频率间隔采样,并且有限长序列的离散傅里叶变换和周期序列的离散傅里叶级数本质是一样的。
2、推导
离散傅里叶级数定义为 将上式两端乘以并对n在0~N-1求和可得
因为
所以
这样
用k代替m得
令,则DFS
IDFS
其中都是周期为N的周期序列,DFS[·]表示离散傅里叶级数正变换,IDFS[·]表示离散傅里叶级数反变换。习惯上,对于长为N的周期序列,把0nN-1区间称为主值区,把称为的主值序列,同样也称为的主值序列。
由于,对于周期序列仅有N个样值,对于任何一个周期进行研究就可以得到它的全部信息。在主值区研究与是等价的,因此在主值区计算DFS和DFT是相等的,所以DFT计算公式形式与DFS基本相同。
其关系为,
所以离散傅里叶正变换 0kN-1
3、定义
DFT:设有限长序列x (n) 长为N(0nN-1),其离散傅里叶变换是一个长为N的频率有限长序列(0kN-1),其正变换为
0kN-1
DFT的分辨率:指其能够分辨的最小频率间隔。
频率分辨率主要由数据截断的长度决定,即时间长度的倒数。也可以说由时间窗函数的傅里叶变换,即谱窗的主瓣宽度决定。不同的谱窗的主瓣宽度不同。矩形窗的主瓣宽度最窄,但其副瓣最高(不利于对频率相邻弱信号的分辨),其它常用的窗函数的主瓣宽度与其副瓣高度近似存在反比关系。主瓣窄,副瓣高,有利于相邻强信号的分辨,但不利于相邻弱信号的分辨。主瓣宽,副瓣低不利于相邻强信号的分辨,但可能有利于相邻弱信号的分辨。
4、实质
把有限长序列当做周期序列的主值序列进行DFS变换,x(n)、X(k)的长度均为N,都是N个值,因此二者具有的信息量是相等的。已知x(n)可以唯一确定X(k),已知X(k)可以唯一确定x(n)。
虽然离散傅里叶变换是两个有限长序列之间的变化,但它们是利用DFS关系推导出来的,因而隐含着周期性
作业内容
1、几种信号的频谱分析
(1)自定义DFT函数
function xk=dft_1(xn)
N=length(xn);
WN=exp(-1i*2*pi/N);
n=0:1:N-1;
k=0:1:N-1;
nk=k'*n;
WNnk=WN.^(nk);
xk=xn*WNnk;
| end |
N=input('N=');
n=0:1:N-1;
xn=input('xn=');
Xk=dft_1(xn1,N);
subplot(3,1,1)
stem(n,xn,'.k');
xlabel('n');
axis([0,N,-2.5,2.5]);
w=2*pi*(0:1:2047)/2048;
Xw=xn*exp(-1i*n'*w);
subplot(3,1,2);
plot(w/pi,abs(Xw));
xlabel('w');
axis([0,1,0,N]);
subplot(3,1,3)
k1=0:1:N-1;
w1=2*pi/N*k1;
stem(w1/pi,abs(Xk),'.k');
xlabel('w');
| axis([0,1,0,N]); |
①xn=cos(0.4*pi*n)+sin(0.6*pi*n)
②xn=0.02*n
③xn=heaviside(n)
2、几种因素对频率分辨率的影响
(1)数据长度(取xn=cos(0.4*pi*n)+sin(0.6*pi*n),wn=boxcar(N)矩形窗)
①N=10
②N=20
③N=100
结论:由图可见,数据长度的增长改变了频谱混叠作用,提高了物理分辨率。
(2)补零(取xn=cos(0.4*pi*n)+sin(0.6*pi*n),N=100, wn=boxcar(N) 矩形窗)
①不补零
②补零至N=300
结论:由图可见,补零只改变了Xk的密度,截断函数的频谱混叠作用没有改变。这说明,补零仅仅是提高了计算分辨率,得到的是高密度频谱,而得不到高分辨率谱。
(3)加窗(取xn=cos(0.*4*pi*n)+ sin(0.6*pi*n), N=100)
①矩形窗wn=boxcar(N)
i)窗长度40
ii)窗长度100
②三角形窗wn= triang(N)
i)窗长度40
ii)窗长度100
③汉宁窗wn=hanning(N)
i)窗长度40
ii)窗长度100
④海明窗wn=hamming(N)
i)窗长度40
ii)窗长度100
⑤布拉克曼窗wn=blackman(N)
i)窗长度40
ii)窗长度100
结论:不同的谱窗的主瓣宽度不同。矩形窗的主瓣宽度最窄。窗长度为数据长度时分辨率最高。
3.与MATLAB自带函数对比
(1)幅频特性
①代码
N=100;
n=0:1:N-1;
xn=cos(0.4*pi*n)+sin(0.6*pi*n);
subplot(3,1,1)
stem(n,xn,'.k');
title('时域序列图xn');
xlabel('n');
axis([0,N,-2.5,2.5]);
Xw=fft(xn,N);
subplot(3,1,2)
k1=0:1:N-1;
w1=2*pi/N*k1;
stem(w1/pi,abs(Xw),'.k');
title('fft幅频特性');
xlabel('频率');
axis([0,1,0,N]);
Xk=dft_1(xn,N);
subplot(3,1,3)
k2=0:1:N-1;
w2=2*pi/N*k2;
stem(w2/pi,abs(Xk),'.k');
title('dft幅频特性');
xlabel('频率');
| axis([0,1,0,N]); |
(2)相频特性
①代码
N=100;
n=0:1:N-1;
xn=cos(0.4*pi*n)+sin(0.6*pi*n);
subplot(3,1,1)
stem(n,xn,'.k');
title('时域序列图xn');
xlabel('n');
axis([0,N,-2.5,2.5]);
Xw=fft(xn,N);
subplot(3,1,2)
k2=0:1:N-1;
w2=2*pi/N*k2;
stem(w2/pi,angle(Xw),'.k');
title('fft相频特性');
xlabel('频率');
axis([0,1,0,N]);
Xk=dft_1(xn,N);
subplot(3,1,3)
k1=0:1:N-1;
w1=2*pi/N*k1;
stem(w1/pi,angle(Xk),'.k');
title('dft相频特性');
xlabel('频率');
| axis([0,1,0,N]); |
结论:dft与MATLAB自带函数相比,相频特性不同,幅频特性无显著区别。下载本文