《数字信号处理》课程研究性学习报告
姓名
学号
同组成员
指导教师
时间
DFT近似计算信号频谱专题研讨
【目的】
(1) 掌握利用DFT近似计算不同类型信号频谱的原理和方法。
(2) 理解误差产生的原因及减小误差的方法。
(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。
【研讨题目】基本题
1．利用DFT分析x(t)=A cos(2πf1t)+B cos(2πf2t)的频谱，其中f1=100Hz，f2=120Hz。
(1)A=B=1; (2)A=1,B=0.2。要求选择不同的窗函数。
【题目分析】
1.对于第一小问，A=B=1，抽样频率应大于最高频率的2倍，才能避免频率混叠，另外由于
信号无限长，所以采用矩形窗进行截短，要想分辨f1,f2两个频率，应满足N≧f sam/△f
2.第二问中f2频率信号比较弱，如果也采用矩形窗，会使得频率泄漏比较大，无法检测到f2频率分量，因此应选择旁瓣较小的Hamming窗
【仿真结果】
【结果分析】
对实验结果进行比较，总结出选择合适DFT参数的原则。【自主学习内容】
【阅读文献】
【发现问题】(专题研讨或相关知识点学习中发现的问题)：
【问题探究】
【仿真程序】
N=30;
L=512;
f1=100;f2=120;fs=500;
T=1/fs;
ws=2*pi*fs;
t=(0:N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
X=fftshift(fft(x,L));
w=(-ws/2+(0:L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('幅度谱')
N=30;
L=512;
f1=100;f2=120;fs=500;
T=1/fs;
ws=2*pi*fs;
t=(0:N-1)*T;
x=cos(2*pi*f1*t)+0.15*cos(2*pi*f2*t);
wh=(hamming(N))';
x=x.*wh;
X=fftshift(fft(x,L));
w=(-ws/2+(0:L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('幅度谱')
【研讨题目】基本题
2．已知一离散序列为
=
=k
[
k
x
k
]
sin(
,
31
2.0
,1,0
),
π
(1)用L=32点DFT计算该序列的频谱，求出频谱中谱峰的频率；
(2)对序列进行补零，然后分别用L=64、128、256、512点DFT计算该序列的频谱，求出频谱中谱峰的频率；
(3)讨论所获得的结果，给出你的结论。该结论对序列的频谱计算有何指导意义？
【题目分析】
本题讨论补零对离散序列频谱计算的影响。
【温磬提示】
在计算离散非周期序列频谱时常用Ω/π作为横坐标，称Ω/π为归一化频率(normalized frequency)。在画频谱时需给出横坐标。每幅图下都需给出简要的文字说明。
由于离散非周期序列频谱是周期的，所以在计算时不需要用fftshift 函数对fft计算的结果进行重新排列。
【序列频谱计算的基本方法】
【仿真结果】
谱峰
14.5561
谱峰的频率
3
DFT点数:L=64 谱峰
14.5561
谱峰的频率
6
DFT点数:L=128 谱峰
15.2671
谱峰的频率
13
DFT点数:L=256 谱峰
15.2671
谱峰的频率
26
DFT点数:L=512
谱峰
15.3563
谱峰的频率
51
【结果分析】
【自主学习内容】
【阅读文献】
【发现问题】(专题研讨或相关知识点学习中发现的问题)：【问题探究】
【仿真程序】
1．
N=32;
k=0:N-1;
xk=sin(0.2*pi*k);
Xm=abs(fft(xk,N));
stem(2*k/N,Xm);
xlabel(' Normalized frequency of 32 DFT');
ylabel('幅度谱');
grid on;
[a1 b1]=max(Xm);
disp('谱峰');disp(a1');
disp('谱峰的频率');disp(k(b1)');
谱峰
14.5561
谱峰的频率
3
2.
N=32
k=0:N-1;
xk=sin(0.2*pi*k);
L=input('DFT点数:L=');
Xm=abs(fft(xk,L));
m=0:L-1;
stem(2*m/L,Xm);
xlabel(['normalized frequency of L=',num2str(L)]);hold on; ylabel('幅度谱');
grid on;
[a1 b1]=max(Xm);
disp('谱峰');disp(a1');
disp('谱峰的频率');disp(m(b1)');
谱峰
14.5561
谱峰的频率
3
DFT点数:L=64
谱峰
14.5561
谱峰的频率 6
DFT 点数:L=128 谱峰 15.2671
谱峰的频率 13
DFT 点数:L=256 谱峰 15.2671
谱峰的频率 26
DFT 点数:L=512 谱峰 15.3563
谱峰的频率 51
（3）.DFT 点数越多，分辨率越高，补零只能使序列的频谱变得细致,但不能提高序列的频率分辨率,
只有采集更多的有效数据,才能得到序列的高分辨率频谱。
【研讨题目】 基本题
3．已知一离散序列为 x [k ]=A cos Ω0k +B cos ( (Ω0+?Ω)k )。用长度N =64的哈明窗对信号截短后近似计
算其频谱。试用不同的A 和B 的取值，确定用哈明窗能分辩的最小的谱峰间隔N
c π2Δw =Ω中
c 的值。（M2-3）
【题目分析】
本题讨论用哈明窗计算序列频谱时的频率分辨率。
【仿真结果】
【结果分析】
当两个相近频率幅度相差不大时，分辨效果比较好。C从2到1.8，1.6，1.4变化，可知c>1.6时能分辨开，并且c>2时较为理想
【自主学习内容】
【阅读文献】
【发现问题】(专题研讨或相关知识点学习中发现的问题)：
【问题探究】
在离散序列频谱计算中为何要用窗函数？用不同的窗函数对计算结果有何影响？与矩形窗相比哈明窗有何特点?如何选择窗函数？
由于计算机存储空间有限，因此我们需要利用窗函数来截取，得到短的序列近似计算频谱
哈明窗旁瓣带走能量少，不同频率之间影响小，但其主瓣宽，因而影响分辨率。所以c的选取根据分辨率要求来选择。
【仿真程序】
w0=0.2*pi;
dw1=pi/32;
dw2=pi/35.6;
dw3=pi/40;
dw4=pi/45.7;
k=0:63;
A=1;
B=1;
L=1028;
x1=(A*cos(w0*k)+B*cos((w0+dw1)*k)).*(0.54-0.46*(cos(2*pi/63*k)));
x2=(A*cos(w0*k)+B*cos((w0+dw2)*k)).*(0.54-0.46*(cos(2*pi/63*k)));
x3=(A*cos(w0*k)+B*cos((w0+dw3)*k)).*(0.54-0.46*(cos(2*pi/63*k)));
x4=(A*cos(w0*k)+B*cos((w0+dw4)*k)).*(0.54-0.46*(cos(2*pi/63*k)));
X1=fftshift(fft(x1,L));
X2=fftshift(fft(x2,L));
X3=fftshift(fft(x3,L));
X4=fftshift(fft(x4,L));
w=-pi:2*pi/L:pi-2*pi/L;
subplot(2,2,1);
plot(w/pi,abs(X1));
subplot(2,2,2);
plot(w/pi,abs(X2));
subplot(2,2,3);
plot(w/pi,abs(X3));
subplot(2,2,4);
plot(w/pi,abs(X4));
【研讨题目】基本题
4．已知一离散序列为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) 给出一种能分辨出信号中两个谱峰的计算方案，并进行仿真实验。（M2-4）
【题目分析】
分析影响谱峰分辨率的主要因数，进一步认识补零在在频谱计算中的作用。
【仿真结果】
k=0:63;L=512;
【结果分析】
即使在序列后面补零仍不能分辨出两个谱峰。由于在后面继续补零，序列的不会有任何变化。由于L的增加，显示更多的细节，这样能减小栅栏现象的影响，但连续信号在离散化和加窗过程中，由于混叠或泄漏等过程造成的信号频谱中信息的失真，则无论如何补0也无法再恢复已损失的信息。
提高分辨率，可以通过增加N的值实现。
【自主学习内容】
【阅读文献】
【发现问题】(专题研讨或相关知识点学习中发现的问题)：
【问题探究】
【仿真程序】
k=0:63;
L=64;
w0=0.4*pi;
w1=w0+pi/64;
x=cos(w0*k)+0.75*cos(w1*k);
x=[x zeros(1,L-length(x))];
X=fftshift(fft(x));
w=-pi:2*pi/L:pi-2*pi/L; plot(w/pi,abs(X));
k=0:63; L=128;
w0=0.4*pi; w1=w0+pi/64;
x=cos(w0*k)+0.75*cos(w1*k); x=[x zeros(1,L-length(x))]; X=fftshift(fft(x));
w=-pi:2*pi/L:pi-2*pi/L; plot(w/pi,abs(X));
k=0:63; L=512;
w0=0.4*pi; w1=w0+pi/64;
x=cos(w0*k)+0.75*cos(w1*k); x=[x zeros(1,L-length(x))]; X=fftshift(fft(x));
w=-pi:2*pi/L:pi-2*pi/L; plot(w/pi,abs(X));
k=0:127; L=128;
w0=0.4*pi; w1=w0+pi/64;
x=cos(w0*k)+0.75*cos(w1*k); x=[x zeros(1,L-length(x))]; X=fftshift(fft(x));
w=-pi:2*pi/L:pi-2*pi/L; plot(w/pi,abs(X));
【研讨题目】 基本题
5．试用DFT 近似计算高斯信号)exp()(2
dt t g -=的频谱抽样值。高斯信号频谱的理论值为
)4ex p(π)j (2
d
d G ω
ω-
=
通过与理论值比较，讨论信号的时域截取长度和抽样频率对计算误差的影响。（M2-6）
【题目分析】
连续非周期信号频谱计算的基本方法。计算中出现误差的主要原因及减小误差的方法。【仿真结果】
【结果分析】
由于信号及频谱都有理论表达式，在进行误差分析时希望给出一些定量的结果。
抽样频率的影响：
非带限信号会发生频谱混叠。又理论值可以看出，随频率的平方的增加而指数衰减。对比实验1，2.可知，提高抽样频率有助于结果更加精确，但同时对计算机也提出了更高的要求。
时域截取长度的影响：
时域截取会发生泄漏现象，本题使用矩形窗截取，泄漏现象其实比较明显。对比实验2，3可知，时域截取长度越长，结果与理论值越接近。
【自主学习内容】
【阅读文献】
【发现问题】(专题研讨或相关知识点学习中发现的问题)：
d值的选取问题。由理论计算公式可知，d值选为pi更为方便。
同时，公式可知，进行DFT变换之前函数还要乘T，这样才能与G（jw）对应。
【问题探究】
【仿真程序】
N=2;
d=pi;
L=64;
fs=1;
T=1/fs;
t=(-N/2:N/2-1)*T;
g=exp(-d*(t.^2));
w=-pi/T:2*pi/L/T:(pi-2*pi/L)/T; G1=T*fftshift(fft(g,L));
G2=sqrt(pi/d)*exp(-(w.^2)/4/d); G3=abs(G1-G2);
subplot(3,1,1)
plot(w,abs(G1));
subplot(3,1,2)
plot(w,abs(G2));
subplot(3,1,3)
plot(w,abs(G3));
N=2;
d=pi;
L=64;
fs=2;
T=1/fs;
t=(-N/2:N/2-1)*T;
g=exp(-d*(t.^2));
ws=2*pi*fs;
w=-pi/T:2*pi/L/T:(pi-2*pi/L)/T; G1=T*fftshift(fft(g,L));
G2=sqrt(pi/d)*exp(-(w.^2)/4/d); G3=abs(G1-G2);
subplot(3,1,1)
plot(w,abs(G1));
subplot(3,1,2)
plot(w,abs(G2));
subplot(3,1,3)
plot(w,abs(G3));
N=4;
d=pi;
L=64;
fs=2;
T=1/fs;
t=(-N/2:N/2-1)*T;
g=exp(-d*(t.^2));
ws=2*pi*fs;
w=-pi/T:2*pi/L/T:(pi-2*pi/L)/T; G1=T*fftshift(fft(g,L));
G2=sqrt(pi/d)*exp(-(w.^2)/4/d); G3=abs(G1-G2);
subplot(3,1,1)
plot(w,abs(G1));
plot(w,abs(G2)); subplot(3,1,3) plot(w,abs(G3));