DSP研究性学习报告频谱计算.doc

《数字信号处理》课程研究性学习报告

指导教师 薛健
时间 2014年4月
DFT近似计算信号频谱专题研讨
【目的】
(1) 掌握利用DFT近似计算不同类型信号频谱的原理和方法。
(2) 理解误差产生的原因及减小误差的方法。
(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。
【研讨题目】 基本题
1. 已知一离散序列为

(1)用L=32点DFT计算该序列的频谱，求出频谱中谱峰的频率；
(2)对序列进行补零，然后分别用L=64、128、256、512点DFT计算该序列的频谱，求出频谱中谱峰的频率；
(3)讨论所获得的结果，给出你的结论。该结论对序列的频谱计算有何指导意义？

【题目分析】
本题讨论补零对离散序列频谱计算的影响。

【温磬提示】
在计算离散非周期序列频谱时常用W/p作为横坐标，称W/p为归一化频率(normalized frequency)。在画频谱时需给出横坐标。每幅图下都需给出简要的文字说明。
由于离散非周期序列频谱是周期的，所以在计算时不必用fftshift 函数对fft计算的结果进行重新排列。

【序列频谱计算的基本方法】
（1） 周期为N的离散周期信号的频谱为

（2） 离散非周期信号的频谱为
,
（3） 有限长N的序列的频谱为
， m=0,1,…,N-1

【仿真结果】
（1）

（2）

【结果分析】
在有限长序列后补0，不会增加任何信息，补0前后的两序列对应的DTFT完全一致，但补0后的DFT存在明显差别。从信号表示的角度来看，对于长度为N的时域序列x[k],可由N点的DFT对应频域序列X[m]唯一表示，X[m]是序列x[k]的离散时间Fourier变换在一个周期2内的等间隔抽样，由于抽样间隔不同所以X[m]不同。但从信号频谱分析的角度，在序列x[k]后补0，可以在的一个周期2内获得更多的抽样值，获得更多信息。

【自主学习内容】

【阅读文献】

【发现问题】 (专题研讨或相关知识点学习中发现的问题)：

【问题探究】

2 已知一离散序列为 x [k]=AcosW0k+Bcos ( (W0+DW)k)。用长度N=64的哈明窗对信号截短后近似计算其频谱。试用不同的A和B的值(如 A和B近似相等，A和B近差距较大)，确定用哈明窗能分辩的最小的谱峰间隔中c的值。

【题目分析】
本题讨论的是Hamming窗截断信号后信号频谱的频率分辨率。
首先，加窗处理对频谱分析主要有两个方面的影响：频谱中出现多余的高频分量，这是窗函数突然截断而引起的，称为频率泄漏现象；谱线变成了有一定宽度的谱峰，谱峰的宽度与信号的长度成反比，所取信号的长度越长，谱峰的宽度就越窄，当信号中两个不同频率分量的频率差𝛥𝑓小于谱峰宽度时，计算出的频谱可能显示不出两个明显的峰值。为使计算出的频谱能显示出相邻的谱峰，则要求相邻频率分量的频率差𝛥𝑓大于谱峰的有效宽度，即：
𝛥𝑓≥𝛥𝑓𝑤 =𝛥𝛺𝑤/2𝜋𝑇
所以，本题采用Hamming窗截断信号在进行频谱分析。Hamming窗是以降低频率分辨率为代价来减少频率泄漏的，使用Hamming窗时要求能分辨的谱峰的间隔为：
𝛥𝑓≥𝑐𝛥𝑓𝑤 =𝑐∕𝑇𝑃=𝑐𝑓𝑠𝑎𝑚/𝑁
按照题目要求，选取𝛺0=0.5，再分别令A=1、B=1和A=1、B=0.1，不断改变c值，观察频谱的变化，即可确定刚好可以分辨谱峰的c值。

【仿真结果】
当A=1、B=1时，c取1、2、3、4时的频谱如下图所示。其中w=𝛥𝛺𝑤。
可以看出，c=2时两谱峰不能分辨，c=3、c=4时可以分辨。

当A=1、B=0.1时，c取1、2、3、4时的频谱如下图所示。其中w=𝛥𝛺𝑤。

可以看出，c=2时两谱峰不能分辨，c=3、c=4时可以分辨。

【结果分析】
将实验结果与教材中定义的哈明窗有效宽度相比较，发表你的看法。

Hamming窗函数的幅值有中心向两端逐渐减弱，因而其高频分量明显减小，频谱中旁瓣的幅度较小，主瓣峰值与第一个旁瓣峰值相对衰减很大可达40dB，Hamming窗是以降低频率分辨率为代价来减少频率泄漏的。
信号的时域加窗会导致在信号的频谱中产生频率泄漏。对同一个窗函数，增加长度N虽然可以减少主瓣宽度，但不减少旁瓣泄漏。为在主瓣宽度与旁瓣泄漏之间取得良好平衡，应根据实际信号的特性采用合适的窗函数。

【自主学习内容】

【阅读文献】

【发现问题】 (专题研讨或相关知识点学习中发现的问题)：

问题：如果使用矩形窗来对信号进行截短，信号的频谱将如何变化？

解答：当A=1、B=1时，使用矩形窗截短信号的频谱如图。

当A=1、B=0.1时，使用矩形窗截短信号的频谱如图。

可以看出，Haming窗可以有效减少旁瓣泄漏，但会造成分辨率下降。且当存在小信号也就是信号中存在较弱的频率分量时，Haming窗比矩形窗更好。

【问题探究】
在离散序列频谱计算中为何要用窗函数？用不同的窗函数对计算结果有何影响？与矩形窗相比哈明窗有何特点?如何选择窗函数？

如果连续信号无限长，则离散后的序列也无限长，无法适用于DFT分析，需要对其进行加窗截短使之成为有限长序列，既可以对信号进行DFT分析，也可以减少计算机的计算量，提取出我们想要的重点信息。
对于不同的信号需要选择不同的窗函数，否则可能会因为窗函数选择的不合理造成较大的失真，从而无法得到符合真实情况的有用信息。
与矩形窗相比，Haming窗主瓣宽度增加，频率分辨率下降，但是旁瓣的能量也随之下降，减少了频率泄漏。但是Hamming窗是以降低频率分辨率为代价来减少频率泄漏的。
根据本题的研究结果，可知当存在小信号也就是信号中存在较弱的频率分量时，Haming窗比矩形窗更好。

3 已知一离散序列为 x[k]=cos(W0k)+0.75cos(W1k), 0£ k £ 63 其中W0=0.4p, W1=W0+p/64
(1) 对x[k]做64点FFT, 画出此时信号的频谱。
(2) 如果(1)中显示的谱不能分辨两个谱峰，是否可对(1)中的64点信号补零而分辨出两个谱峰。通过编程进行证实，并解释其原因 。
(3) 给出一种能分辨出信号中两个谱峰的计算方案，并进行仿真实验。

【题目分析】
分析影响谱峰分辨率的主要因数，进一步认识补零在在频谱计算中的作用。

本题目具体探究了关于信号分辨率的具体问题。进一步研究了关于信号处理中参数的选择问题。
我们知道在进行信号的频域分析的时候需要一下4个参数
1截取信号长度N
2抽样频率fs
3采样点数L
4选择什么样的窗
本题目中选择的矩形窗，在我们常用的三大类窗中（矩形窗，三角窗，指数窗）是主瓣宽度最小的一种窗。而本题目中给出的信号所需要的分辨频率p/64此时c=0.5<1，而且本题目中两个信号的幅度差不多没有小信号所以如果修改窗函数类型没法达到预期目的。
首先我们按照题目提示进行补零实验。
改变采样点数L只能加大频域的取样点数补零起到对原 X(K ) 做插值的作用一方面克服“栅栏“效应使谱的外观得到平滑另一方面由于对数据截断时所引起的频谱泄露有可能在频谱中出现一些难以确认的谱峰补零后有可能消除这种现象。但是这种方法只能改变频谱的可视分辨率。对于实际分辨率没有贡献。
所以解决这个问题应该从截取信号的长度即窗的宽度下手，改变矩形窗主瓣宽度。

【仿真结果】

【结果分析】

(1) 仿真结果如图1所示，即信号频谱图。
(2) 由（1）结果可知图中显示的谱不能分辨两个谱峰，进行补零后频谱图如图2所示，由图可知并不能使分辨率提高，依然分辨不出两个谱峰。
(3) 解决方案：采用哈明窗，取N=64*2，仿真结果如图3所示，可见此时可清晰分辨出两个谱峰。

【自主学习内容】

【阅读文献】

【发现问题】 (专题研讨或相关知识点学习中发现的问题)：

【问题探究】
1、2、3题讨论的是离散信号频谱的计算问题。与连续信号频谱计算问题相比较，其计算误差有何不同？

4 试用DFT近似计算高斯信号的频谱抽样值。高斯信号频谱的理论值为

通过与理论值比较，讨论信号的时域截取长度和抽样频率对计算误差的影响。

【题目分析】
连续非周期信号频谱计算的基本方法。计算中出现误差的主要原因及减小误差的方法。

连续非周期信号频谱与Fourier的关系：
𝑋(e 𝑗𝑤)=ｓａｍ）］＝
在计算中出现3个问题：
1、 混叠现象：当连续信号不是带限信号时，在连续信号离散化时，就会出现信号的频谱混叠。再用DFT分析连续信号频谱时，信号抽样频率f对分析时影响精度较大，因为其直接影响频谱的混叠程度。
2、 泄漏现象：由于矩形窗在其两个端点突然截断，使得频谱中存在许多高频分量，而增加信号的长度不能减少频率泄漏，而且信号的谱峰宽度与信号的长度成正比。
3、 栅栏现象：由于DFT得出的是离散序列，无法反映连续信号频谱中抽样点的细节，为了改善这种现象，需要增加DET的点数。

【仿真结果】

【结果分析】

由于信号在时域和频谱都有理论表达式，在进行误差分析时希望给出一些定量的结果。
当截取时间增大时，谱线逐渐变为谱峰，而且谱峰宽度随截取时间在一定范围变化时内变窄。
当抽样频率变化时，对ＤＦＴ分析频谱的精度影响较大，以为其直接影响频谱混叠的精度。
　　很容易能看出，增加截取时间和抽样频率都可以减小误差。

【自主学习内容】

连续非周期信号的频谱计算方法。

【阅读文献】

【发现问题】 (专题研讨或相关知识点学习中发现的问题)：

【问题探究】

扩展题
5 本题研究连续周期信号频谱的近似计算问题。
周期为T0的连续时间周期信号x(t)可用Fourier级数表示为

其中

X(nw0)称为连续时间周期信号x(t)的频谱函数。称为信号的基频（基波），称为信号的谐波。
如果信号x(t)函数表达式已知，则可由积分得出信号的频谱。如果信号x(t)函数表达式未知，或者x(t)函数表达式非常复杂，则很难由积分得信号的频谱。本题的目的就是研究如何利用DFT近似计算连续时间周期信号的频谱。
(1)若在信号x(t)的一个周期T0内抽样N个点，即， T为抽样周期(间隔)，可获得序列x[k]

试分析序列x[k]的DFT与连续时间周期信号x(t)的频谱X(nw0)的关系；
(2)由(1)的结论，给出由DFT近似计算周期信号频谱X(nw0)的方案；
(3)周期信号x(t)的周期T0=1，x(t)在区间[0，1]的表达式为
x(t)=20t2(1-t)4cos(12pt)
(a)试画出信号x(t)在区间[0，1]的波形；
(b)若要用4次以内的谐波（包含4次）近似表示x(t)，试给出计算方案，并计算出近似表示的误差。讨论出现误差的原因及减小误差的方法。

【题目分析】

【理论推导】
DFT计算所得结果X[m]与连续周期信号频谱X(nw0)的关系。

【计算方案】
根据理论推导结果设计近似计算方案。分析产生误差的主要原因。

【扩展分析】
如果周期信号x(t)是带限信号，即信号的最高频率分量为Mw0(是正整数)，试确定在一个周期内的最少抽样点N，使得在频谱的计算过程当中不存在混叠误差。与抽样定理给出的结论比较，发表你的看法。

【仿真结果】

【结果分析】
讨论DFT点数对近似计算的影响，讨论所取谐波项的多少对近似计算的影响。误差分析要给出定量的结果，如平均误差，最大误差等。
与连续非周期信号频谱计算过程中存在的误差相比较，连续周期信号频谱的计算计算误差有何异同？
【自主学习内容】

【阅读文献】

【发现问题】 (专题研讨或相关知识点学习中发现的问题)：

【问题探究】

【仿真程序】

1、
（1）
k=0:31;
x=sin(0.2*pi*k);
x_32=fft(x,32);
L=0:31;
plot(L/32,abs(x_32));
xlabel('Normalized frequency');ylabel('Magnitude');
title('L=32');
（2）
k=0:31;
x=sin(0.2.*pi.*k);
X_64=fft(x,64);
L=0:63;
subplot(2,2,1);
plot(L/64,abs(X_64));
xlabel('Normalized frequency');ylabel('Magnitude');
title('L=64');
X_128=fft(x,128);
L=0:127;
subplot(2,2,2);
plot(L/128,abs(X_128));
xlabel('Normalized frequency');ylabel('Magnitude');
title('L=128');
X_256=fft(x,256);
L=0:255;
subplot(2,2,3);
plot(L/256,abs(X_256));
xlabel('Normalized frequency');ylabel('Magnitude');
title('L=256');
X_512=fft(x,512);
L=0:511;
subplot(2,2,4);
plot(L/512,abs(X_512));
xlabel('Normalized frequency');ylabel('Magnitude');
title('L=512');

2、
A=1;
B=1;
k=0:63;%ÊäÈëÐźų¤¶ÈΪ64
l=0:511;
%c=input('c=');ÊäÈëcÖµ
N=64;
L=512;
subplot(221);
c=1;
x=A*cos(0.5*pi*k)+B*cos((0.5*pi+c*2*pi/N)*k);%ÊäÈëÐźÅ
h=(hamming(N))';%³¤¶ÈΪ64µÄHamming´°
x=x.*h; %ÓÃHamming´°½Ø¶ÏÐźÅ
X=fft(x,L);%¿ìËÙFourier±ä»»
plot(l/L,abs(X));%»­³öͼÏñ
xlabel('¹éÒ»»¯ÆµÂÊ');
ylabel('·ùÖµ');
w=2*c*pi/N;
title(strcat('w= ',num2str(w), ' ; c= ',num2str(c)));
subplot(222);
c=2;
x=A*cos(0.5*pi*k)+B*cos((0.5*pi+c*2*pi/N)*k);%ÊäÈëÐźÅ
h=(hamming(N))';%³¤¶ÈΪ64µÄHamming´°
x=x.*h; %ÓÃHamming´°½Ø¶ÏÐźÅ
X=fft(x,L);%¿ìËÙFourier±ä»»
plot(l/L,abs(X));%»­³öͼÏñ
xlabel('¹éÒ»»¯ÆµÂÊ');
ylabel('·ùÖµ');
w=2*c*pi/N;
title(strcat('w= ',num2str(w), ' ; c= ',num2str(c)));
subplot(223);
c=3;
x=A*cos(0.5*pi*k)+B*cos((0.5*pi+c*2*pi/N)*k);%ÊäÈëÐźÅ
h=(hamming(N))';%³¤¶ÈΪ64µÄHamming´°
x=x.*h; %ÓÃHamming´°½Ø¶ÏÐźÅ
X=fft(x,L);%¿ìËÙFourier±ä»»
plot(l/L,abs(X));%»­³öͼÏñ
xlabel('¹éÒ»»¯ÆµÂÊ');
ylabel('·ùÖµ');
w=2*c*pi/N;
title(strcat('w= ',num2str(w), ' ; c= ',num2str(c)));
subplot(224);
c=4;
x=A*cos(0.5*pi*k)+B*cos((0.5*pi+c*2*pi/N)*k);%ÊäÈëÐźÅ
h=(hamming(N))';%³¤¶ÈΪ64µÄHamming´°
x=x.*h; %ÓÃHamming´°½Ø¶ÏÐźÅ
X=fft(x,L);%¿ìËÙFourier±ä»»
plot(l/L,abs(X));%»­³öͼÏñ
xlabel('¹éÒ»»¯ÆµÂÊ');
ylabel('·ùÖµ');
w=2*c*pi/N;
title(strcat('w= ',num2str(w), ' ; c= ',num2str(c)));

3、
(1)
W0=0.4*pi;
W1=W0+pi/64;
k=0:63;
x=cos(W0*k)+0.75*cos(W1*k);
y=fft(x);
yy=fftshift(abs(y));
k=linspace(-pi,pi,length(yy));
stem(k,yy,’.’);
(2)
W0=0.4*pi;
W1=W0+pi/64;
k=0:63;
x=cos(W0*k)+0.75*cos(W1*k);
y=fft(x,1024);
yy=fftshift(abs(y));
k=linspace(-pi,pi,length(yy));
plot(k,yy);
(3)
W0=0.4*pi;
W1=W0+pi/64;
N=64*2;
k=0:N-1;
wh=(hamming(N))';
x=(cos(W0*k)+0.75*cos(W1*k)).*wh;
y=fft(x,1024);
yy=fftshift(abs(y));
k=linspace(-pi,pi,length(yy));
plot(k,yy);

4、
d = pi;
fsam_set = [1 2 4 8 16];
N_set = [1 2 4 8];
L = 64;
for fsam_index = 1:length(fsam_set);
fsam = fsam_set(fsam_index);
Tsam = 1/fsam;
for N_index = 1:length(N_set);
N = N_set(N_index);
cuttime = N*Tsam;
k = -N/2:N/2-1;
t = k*Tsam;
x = exp(-d*t.^2);
X = Tsam*fftshift(fft(x,L));
w = linspace(-pi*fsam,pi*fsam,L);
subplot(length(N_set),length(fsam_set),fsam_index+length(fsam_set)*(N_index-1));
plot(w,abs(X));
hold on;
G = exp(-w.*w/(4*d));
plot(w,G,'r');
title(['cuttime=' num2str(cuttime) ' fsam=' num2str(fsam)]);
end
end

电子文件提交说明:
文件名 学号＋姓名＋研讨题目 (每组在文件名上只需提供一个学号和姓名)
文件格式 Microsoft Word 2003 或以下版本
建议用Office 套装软件 Visio 进行绘图。 Matlab所绘图形可在Visio环境下ungroup 后进行编辑。