数字信号处理研讨DFT近似计算信号频谱.docx
《数字信号处理研讨DFT近似计算信号频谱.docx》由会员分享，可在线阅读，更多相关《数字信号处理研讨DFT近似计算信号频谱.docx（18页珍藏版）》请在冰豆网上搜索。
数字信号处理研讨DFT近似计算信号频谱
BeijingJiaotongUniversity
数字信号处理研讨
DFT近似计算信号频谱
学院：
电子信息工程学院
小组成员：
指导教师：
时间：
2013.06.24
DFT近似计算信号频谱专题研讨
【目的】
（1）掌握利用DFT近似计算不同类型信号频谱的原理和方法。
（2）理解误差产生的原因及减小误差的方法。
（3）培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。
【研讨题目】基本题
1.已知一离散序列为
（1）用L=32点DFT计算该序列的频谱，求出频谱中谱峰的频率；
（2）对序列进行补零，然后分别用L=64、128、256、512点DFT计算该序列的频谱，求出频谱中谱峰的频率；
（3）讨论所获得的结果，给出你的结论。
该结论对序列的频谱计算有何指导意义？
【题目分析】
本题讨论补零对离散序列频谱计算的影响。
【温磬提示】
在计算离散非周期序列频谱时常用Ω/π作为横坐标，称Ω/π为归一化频率（normalizedfrequency）。
在画频谱时需给出横坐标。
每幅图下都需给出简要的文字说明。
由于离散非周期序列频谱是周期的，所以在计算时不必用fftshift函数对fft计算的结果进行重新排列。
【序列频谱计算的基本方法】
连续信号，通过其抽样的离散信号，和离散信号的DFT变换存在如下关系：
通过如上关系，我们就可以通过DFT来求信号的频谱。
【仿真结果】
fm=
0.9375
fm=
0.9219
fm=
0.9141
fm=
0.9102
fm=
0.9082
【结果分析】
增加DFT的点数可以使频谱更容易观察，即减轻了栅栏效应带来的影响。
频谱的横坐标为归一化频率，所以原信号的峰值第一次应该出现在0.2处，随着DFT点数的增大，频谱表示也越来越精确。
【自主学习内容】
1.归一化频率相关知识。
2.通过matlab计算DFT和matlab的绘图操作。
【阅读文献】
1．数字信号处理
2．补零对有限长序列频谱及DFT的影响。
【发现问题】（专题研讨或相关知识点学习中发现的问题）：
DFT点数的增多是否能提高频谱分辨率？
【问题探究】
虽然DFT点数增加使图像更加细致，但是因为不论DFT点数是多少，抽样点数都是相同的，所以每个频谱所包含的信息相同，频谱分辨率只与抽样点数相关，与DFT点数无关，频谱分辨率相同。
所以不能通过增大DFT点数而减少信息损失。
DFT点数的增多不能提过频率分辨率。
【仿真程序】
k=0:
31;
N=16;
x=sin（0.2*pi*k）;
fori=1:
5
N=2*N;
X=fft（x,N）;
k=0:
N-1;
subplot（5,1,i）;
plot（k/N,abs（X））;
xlabel（'f/HZ'）;
ylabel（'Magnitude'）;
fm=find（X==max（X））/N;
fm
end
2已知一离散序列为x[k]=AcosΩ0k+Bcos（（Ω0+∆Ω）k）。
用长度N=64的哈明窗对信号截短后近似计算其频谱。
试用不同的A和B的值（如A和B近似相等，A和B差距较大），确定用哈明窗能分辩的最小的谱峰间隔
中c的值。
【题目分析】
本题讨论用哈明窗计算序列频谱时的频率分辨率。
【仿真结果】
【结果分析】
第一排是c=4，第二排c=3，远大于教材定义的标准c=2。
无论A/B的比值是多少，区分谱峰毫无压力。
第3排中，c=2是教材定义的标准，A与B接近的时候，区分较为容易，当A/B为1.6的时候，几乎是无法区分，而A/B=4的时候，完全无法区分。
所以说在设计使用哈明窗的时候，如果不同频谱的信号所占比例相差较大，哈明窗的c要比2大。
从第一列可以看出，如果A/B象近的时候，即使c稍小于2，也是可以区分的。
每一列中，A/B的比值是相同的，随着c的减小，谱峰的区别能力很容易看出是降低的。
【自主学习内容】
不同窗函数的屏幕分辨率。
【阅读文献】
数字信号处理
【发现问题】（专题研讨或相关知识点学习中发现的问题）：
哈明窗所起的作用？
哈明窗与矩形窗谁的分辨率高？
在仿真的时候c该选用怎样的值？
【问题探究】
哈明窗主瓣较宽，分辨率比矩形窗差，但是哈明窗可以使能量集中在主瓣处。
所以哈明窗可以起到减小泄漏的作用。
在仿真的时候c该选用怎样的值？
在开始的时候选用的c为124816，对比效果不明显，发现c的取值过于大，反复尝试确定在2附近多取值，其他值也不可取得过于大，因为c稍大实验现象基本相同，所以没有意义。
【仿真程序】
N=64;
k=0:
N-1;
c_set=[432.221.91.81.5];
A_set=[11.21.64];
B=1;
w0=0.5*pi;
forc_index=1:
length（c_set）
c=c_set（c_index）;
dw=c*2*pi/N;
forA_index=1:
length（A_set）
A=A_set（A_index）;
x=A*cos（w0*k）+B*cos（（w0+dw）*k）;
Xh=fft（x.*hamming（N）'）;
subplot（length（c_set）,length（A_set）,length（A_set）*（c_index-1）+A_index）;
plot（k,abs（Xh））;
holdon;
title（['c='num2str（c）';A/B='num2str（A）]）;
axis（[1225025]）;
end
end
3已知一离散序列为x[k]=cos（Ω0k）+0.75cos（Ω1k）,0£k£63其中Ω0=0.4π,Ω1=Ω0+π/64
（1）对x[k]做64点FFT,画出此时信号的频谱。
（2）如果
（1）中显示的谱不能分辨两个谱峰，是否可对
（1）中的64点信号补零而分辨出两个谱峰。
通过编程进行证实，并解释其原因。
（3）给出一种能分辨出信号中两个谱峰的计算方案，并进行仿真实验。
【题目分析】
分析影响谱峰分辨率的主要因数，进一步认识补零在在频谱计算中的作用。
【仿真结果】
【结果分析】
频谱的横坐标是归一化频率，图示只观察要区别的部分。
每一个列的抽样点相同，每一行DFT点相同。
但看第一列，虽然增加DFT点数，甚至是增加到32768，也不能够区分两个谱峰，所以要想增加频率的分辨率，可以增加抽样点数，第二列和第三列抽样点数分别是128和256，都可以区分两个谱峰。
【自主学习内容】
如何增加频谱分辨率。
【阅读文献】
数字信号处理。
【发现问题】（专题研讨或相关知识点学习中发现的问题）：
如何计算所需抽样点数？
可以通过信号频率的差值计算。
【问题探究】
1、2、3题讨论的是离散信号频谱的计算问题。
与连续信号频谱计算问题相比较，其计算误差有何不同？
连续信号的频率是非周期的，离散信号的频谱是连续信号频谱的周期话，可能会有混叠误差。
【仿真程序】
w0=0.4*pi;
dw=pi/64;
w1=w0+dw;
N_set=[64128256];
L_set=[6412825651232768];
forN_index=1:
length（N_set）;
N=N_set（N_index）;
k=0:
N-1;
x=cos（w0*k）+0.75*cos（w1*k）;
forL_index=1:
length（L_set）
L=L_set（L_index）;
ifLcontinue
end
X=fft（x,L）;
m=（0:
L-1）*2/L;
subplot（length（L_set）,length（N_set）,（L_index-1）*length（N_set）+N_index）;
plot（m,abs（X））;
axis（[0.380.44040]）;
title（['N='num2str（N）'L='num2str（L）]）;
end
end
4试用DFT近似计算高斯信号
的频谱抽样值。
高斯信号频谱的理论值为
通过与理论值比较，讨论信号的时域截取长度和抽样频率对计算误差的影响。
【题目分析】
连续非周期信号频谱计算的基本方法。
计算中出现误差的主要原因及减小误差的方法。
【仿真结果】
【结果分析】
如图所示，横轴坐标为角频率，左斜线上的图截取时间相同，垂线上的抽样频率相同，很容易能否看出，增加截取时间和抽样频率都可以减小误差。
【自主学习内容】
连续非周期信号的频谱计算方法。
【仿真程序】
d=pi;
fsam_set=[124816];
N_set=[1248];
L=64;
forfsam_index=1:
length（fsam_set）;
fsam=fsam_set（fsam_index）;
Tsam=1/fsam;
forN_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））;
holdon;
G=exp（-w.*w/（4*d））;
plot（w,G,'r'）;
title（['cuttime='num2str（cuttime）'fsam='num2str（fsam）]）;
end
end
扩展题
5本题研究连续周期信号频谱的近似计算问题。
周期为T0的连续时间周期信号x（t）可用Fourier级数表示为
其中
X（nω0）称为连续时间周期信号x（t）的频谱函数。
称为信号的基频（基波），
称为信号的谐波。
如果信号x（t）函数表达式已知，则可由积分得出信号的频谱。
如果信号x（t）函数表达式未知，或者x（t）函数表达式非常复杂，则很难由积分得信号的频谱。
本题的目的就是研究如何利用DFT近似计算连续时间周期信号的频谱。
（1）若在信号x（t）的一个周期T0内抽样N个点，即
，T为抽样周期（间隔），可获得序列x[k]
试分析序列x[k]的DFT与连续时间周期信号x（t）的频谱X（nω0）的关系；
（2）由
（1）的结论，给出由DFT近似计算周期信号频谱X（nω0）的方案；
（3）周期信号x（t）的周期T0=1，x（t）在区间[0，1]的表达式为
x（t）=20t2（1-t）4cos（12πt）
（a）试画出信号x（t）在区间[0，1]的波形；
（b）若要用6次以内的谐波近似表示x（t），试给出计算方案，并计算出近似表示的误差。
讨论出现误差的原因及减小误差的方法。
【理论推导】
DFT计算所得结果X[m]与连续周期信号频谱X（nω0）的关系。
由于
于是得
令n=m+rN;m=0,1,2,…Ｎ－１，ｒ为整数，上式可化为
化简可得
对比IDFT
可得
【计算方案】
根据理论推导结果设计近似计算方案。
分析产生误差的主要原因。
DFT是连续信号频谱的周期化，如果频谱信号高频信号较多，就会产生混叠，所以有误差。
如果主要是低频信号，周期话后互补产生影响，就没有误差。
【扩展分析】
如果周期信号x（t）是带限信号，即信号的最高频率分量为Mω0（是正整数），试确定在一个周期内的最少抽样点N，使得在频谱的计算过程当中不存在混叠误差。
与抽样定理给出的结论比较，发表你的看法。
很容易看出，要想不产生混叠，必须N2M+1。
该结论与抽样定理相似。
因为周期信号是一种特殊的连续信号，所以抽样定理对周期信号也是适用的。
【仿真结果】
【结果分析】
在N=16,k=6时，误差最小，最大误差为0.1229，平均误差为0.07187。
无论是减小N和减小k都很大的增大了误差。
为了分别测试N和k对误差的影响，分别设k=6，求不同N时候的误差，N=16，不同k对误差的影响，分别得到了下图（在原程序上稍加改动，所以在仿真程序中不再另写了）
此图可以看出，在N等于13的时候误差有明显的减小。
在k=6的时候，误差有明显的减小。
可见原函数最大的高频分量为6*w0，根据前面推出的理论，N取6*2+1=13点的时侯误差减小。
实验结果符合结论
【自主学习内容】
周期信号的频谱与DFT的关系。
【发现问题】：
原函数最高频分量大约是多少？
【问题探究】
根据前面推理，最高频分量为6*w0。
【仿真程序】
n=1024;
t=linspace（0,1,n）;
N_set=[1216];
H_set=[46];
T0=1;
count=0;
forN_index=1:
length（N_set）
N=N_set（N_index）;
dT=T0/N;
k=（0:
N-1）*dT;
x=20*k.^2.*（1-k）.^4.*cos（12*pi.*k）;
X=fft（x）;
forH_index=1:
length（H_set）
H=H_set（H_index）;
sum=0;
X1=[X（N-H+1:
N）X（1:
H+1）]/N;
fork=-H:
H;
yt=X1（k+H+1）*exp（j*2*pi*k*t/T0）;
sum=sum+yt;
end
figure（（N_index-1）*length（H_set）+H_index）;
subplot（1,2,1）;
xt=20*t.^2.*（1-t）.^4.*cos（12*pi.*t）;
plot（t,xt,'r',t,sum,'b'）;
holdon;
title（['N='num2str（N）'k='num2str（H）]）;
subplot（1,2,2）;
plot（t,xt-sum）;4
de=abs（xt-sum）;
maxde=max（de）;
meande=mean（de）;
title（{['Diviation:
'];['mean:
'num2str（meande,4）'max:
'num2str（maxde,4）]}）;
count=count+1;
deviation（count）=meande;
end
end
电子文件提交说明:
文件名学号＋姓名＋研讨题目（每组在文件名上只需提供一个学号和姓名）
文件格式MicrosoftWord2003或以下版本
建议用Office套装软件Visio进行绘图。
Matlab所绘图形可在Visio环境下ungroup后进行编辑。
升级会员 