数字信号处理第二次研学doc.docx
- 文档编号:11084118
- 上传时间:2023-02-24
- 格式:DOCX
- 页数:38
- 大小:881.10KB
数字信号处理第二次研学doc.docx
《数字信号处理第二次研学doc.docx》由会员分享,可在线阅读,更多相关《数字信号处理第二次研学doc.docx(38页珍藏版)》请在冰豆网上搜索。
数字信号处理第二次研学doc
《数字信号处理》课程研究性学习报告
姓名
学号
同组成员
指导教师
时间
DFT近似计算信号频谱专题研讨
【目的】
(1)掌握利用DFT近似计算不同类型信号频谱的原理和方法。
(2)理解误差产生的原因及减小误差的方法。
(3)培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
【研讨题目】基本题
1.利用DFT分析x(t)=Acos(2πf1t)+Bcos(2πf2t)的频谱,其中f1=100Hz,f2=120Hz。
(1)A=B=1;
(2)A=1,B=0.2。
【题目分析】
分析题目,给出合适的DFT参数
由取样定理知,要使信号频谱不混叠,则抽样频率不小于最高频率的两倍。
而要满足信号分辨率的要求,抽样点数N≧fsam/△f。
在对信号做DFT时,由于对信号进行截短,因此会产生频谱泄漏,要想从频谱中很好的分辨出个频率分量,需要考虑时域抽样频率,所加的窗函数,窗函数的长度,以及DFT的点数等参数对结果的影响。
(1)A=B=1,即x(t)=cos(2πf1t)+cos(2πf2t)
矩形窗1:
条件:
fsam=240Hz;N=20;L=512
矩形窗2:
条件:
fsam=600Hz;N=40;L=512
矩形窗3:
fsam=1200Hz;N=80;L=512
Hamming窗1:
N=40;L=512;fs=600;
Hamming窗2:
N=60;L=512;fs=600;
Hamming窗3:
N=120;L=512;fs=600;
(2)A=1,B=0.2,即x(t)=cos(2πf1t)+0.2cos(2πf2t)
矩形窗:
N=100;L=512;fs=600
Hamming窗:
N=100;L=512;fs=600
【仿真结果】
【结果分析】
对实验结果进行分析比较,回答:
加窗对谱分析有何影响?
如何选择合适的窗函数?
选择合适DFT参数的原则?
在
(1)中进行矩形窗仿真时,我们选择了不同的fsam,分别为240,600,1200它们均满足抽样定理,但是我们在实验中却发现,在240hz时出现了混叠现象。
所以,在实际应用中抽样频率应大于最低抽样频率3-5倍才能有更好的结果。
进行hamming窗仿真时,在保证抽样频率相同的条件下,取不同的长度也40,120。
其中,N=40不满足N>=60的要求,我们可以看到出现了混叠,利用N≧fsam/△f,我们可以算出当fs=600时N为60时恰好可以分出频谱,而实际中N=60时无法分出两个频率分量,而N=120时,仿真效果良好
在
(2)的条件下进行仿真时,我们选取了相同的N、L、fsam值,但是分别使用了矩形窗和哈明窗。
使用矩形窗时,幅度较小的峰值与旁瓣的幅度接近,甚至难以区分,效果不理想。
使用哈明窗后,泄露现象被有效遏制,所以可以清楚区分主瓣、旁瓣。
所以,在选择参数进行DFT变换时,应该保证抽样频率满足抽样定理,并且能大于最小抽样值3-5倍。
长度选择保证N≧fsam/△f,且实际中取最小长度即N=fsam/△f时,会出现混叠。
为防止泄露现象,特别是峰值之间差异较大时,应该选择加特殊的窗,如哈明窗。
【发现问题】
按照理论分析最小抽样频率只需要满足2fmax就可以满足抽样定理,但在仿真中发现该频率无法满足要求,频谱发生严重的混叠。
所以抽样频率应为最小抽样频率3-5倍。
另外,在使用哈明窗作为窗函数时,按照理论分析,当fs=600时N为60时恰好可以分出频谱,而实际中N=60时无法分出两个频率分量,当N=90时则可以分出。
因此在做DFT是窗函数长度应大于最小长度。
【仿真程序】
N=20;
L=512;
f1=100;f2=120;fs=240;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
X=fft(x,L);
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('矩形窗1.1');
N=40;
L=512;
f1=100;f2=120;fs=600;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
X=fft(x,L);
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('矩形窗1.2')
N=80;
L=512;
f1=100;f2=120;fs=1200;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
X=fft(x,L);
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('矩形窗1.3')
N=40;
L=512;
f1=100;f2=120;fs=600;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
wh=(hamming(N))';
x=x.*wh;
X=fft(x,L);
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('哈明窗1.1')
:
N=60;
L=512;
f1=100;f2=120;fs=600;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
wh=(hamming(N))';
x=x.*wh;
X=fft(x,L);
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('哈明窗1.2')
N=120;
L=512;
f1=100;f2=120;fs=600;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
wh=(hamming(N))';
x=x.*wh;
X=fft(x,L);
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('哈明窗1.3')
N=100;
L=512;
f1=100;f2=120;fs=600;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);
X=fft(x,L);
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('矩形窗2.1')
N=100;
L=512;
f1=100;f2=120;fs=600;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);
wh=(hamming(N))';
x=x.*wh;
X=fft(x,L);
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X));
ylabel('哈明窗2.1')
【研讨题目】基本题
2.试用DFT近似计算高斯信号
的频谱抽样值。
高斯信号频谱的理论值为
通过与理论值比较,讨论信号的时域截取长度和抽样频率对计算误差的影响。
(M2-6)
【题目分析】
连续非周期信号频谱计算的基本方法。
计算中出现误差的主要原因及减小误差的方法。
【仿真结果】
【结果分析】
由于信号及频谱都有理论表达式,在进行误差分析时希望给出一些定量的结果。
当时域截取长度相同时,抽样间隔越小时误差越小,当抽样间隔一定时,时域截取长度越长,误差越小。
当取抽样间隔为1S,时域截取长度为2S时,误差较大,绝对误差在0.5左右;当抽样间隔为0,5S,时域截取长度为2S时,误差比间隔为1S时小,绝对误差不大于0.2;当抽样间隔为0.5S时域截取长度为4S时,误差更小,绝对误差不大于0.04。
因为时域截取长度越长,保留下来的原信号中的信息越多,抽样间隔越小,频谱越不容易发生混叠,所以所得频谱与理论值相比,误差更小。
【问题探究】
抽样频率的大小对频谱的影响:
信号抽样会导致频谱的周期化,如果对于非带限信号,则会发生频谱混叠。
理论值是随频率平方的增加指数衰减,所以当抽样频率比较大的时候,混叠不会太严重。
时域截取长度对频谱的影响:
时域截取会产生泄露现象,如果时域截取长度比较短,则会产生比较明显的泄露现象。
【仿真程序】
Ts=0.5;
N=2;
N0=64;
k=(-N/2:
(N/2))*Ts;
x=exp(-pi*(k).^2);
X=Ts*fftshift(fft(x,N0));
w=-pi/Ts:
2*pi/N0/Ts:
(pi-2*pi/N0)/Ts;
XT=(pi/pi)^0.5*exp(-w.^2/4/pi);
subplot(2,1,1)
plot(w/pi,abs(X),'-o',w/pi,XT);
xlabel('\omega/\pi');
ylabel('X(j\omega)');
legend('试验值','理论值');
title(['Ts=',num2str(Ts)'''N=',num2str(N)]);
subplot(2,1,2)
plot(w/pi,abs(X)-XT)
ylabel('实验误差')
xlabel('\omega/\pi');
Ts=1;
N=2;
N0=64;
k=(-N/2:
(N/2))*Ts;
x=exp(-pi*(k).^2);
X=Ts*fftshift(fft(x,N0));
w=-pi/Ts:
2*pi/N0/Ts:
(pi-2*pi/N0)/Ts;
XT=(pi/pi)^0.5*exp(-w.^2/4/pi);
subplot(2,1,1)
plot(w/pi,abs(X),'-o',w/pi,XT);
xlabel('\omega/\pi');
ylabel('X(j\omega)');
legend('试验值','理论值');
title(['Ts=',num2str(Ts)'''N=',num2str(N)]);
subplot(2,1,2)
plot(w/pi,abs(X)-XT)
ylabel('实验误差')
xlabel('\omega/\pi');
Ts=0.5;
N=4;
N0=64;
k=(-N/2:
(N/2))*Ts;
x=exp(-pi*(k).^2);
X=Ts*fftshift(fft(x,N0));
w=-pi/Ts:
2*pi/N0/Ts:
(pi-2*pi/N0)/Ts;
XT=(pi/pi)^0.5*exp(-w.^2/4/pi);
subplot(2,1,1)
plot(w/pi,abs(X),'-o',w/pi,XT);
xlabel('\omega/\pi');
ylabel('X(j\omega)');
legend('试验值','理论值');
title(['Ts=',num2str(Ts)'''N=',num2str(N)]);
subplot(2,1,2)
plot(w/pi,abs(X)-XT)
ylabel('实验误差')
xlabel('\omega/\pi');
Ts=1;
N=4;
N0=64;
k=(-N/2:
(N/2))*Ts;
x=exp(-pi*(k).^2);
X=Ts*fftshift(fft(x,N0));
w=-pi/Ts:
2*pi/N0/Ts:
(pi-2*pi/N0)/Ts;
XT=(pi/pi)^0.5*exp(-w.^2/4/pi);
subplot(2,1,1)
plot(w/pi,abs(X),'-o',w/pi,XT);
xlabel('\omega/\pi');
ylabel('X(j\omega)');
legend('试验值','理论值');
title(['Ts=',num2str(Ts)'''N=',num2str(N)]);
subplot(2,1,2)
plot(w/pi,abs(X)-XT)
ylabel('实验误差')
xlabel('\omega/\pi');
【研讨题目】基本题
3.已知一离散序列为
(1)用L=32点DFT计算该序列的频谱,求出频谱中谱峰的频率;
(2)对序列进行补零,然后分别用L=64、128、256、512点DFT计算该序列的频谱,求出频谱中谱峰的频率;
(3)讨论所获得的结果,给出你的结论。
该结论对序列的频谱计算有何指导意义?
【题目分析】本题讨论补零对离散序列频谱计算的影响。
【序列频谱计算的基本方法】
1.解析法——依据定义严格计算
2.用DFT近似计算
【仿真结果】
这是L=32时的DFT
这是L=64、128,256、512时的DFT
L=32、64、128、256、512时用DFT求出的频谱中谱峰对应的频率
(注意:
谱峰对应的频率的单位是πHz)
【结果分析】
1.zeropadding对我们“看到的”频谱是有影响的,但对频谱本身是没有影响。
它没有改变频谱和序列一一对应的关系,但减小了频谱的最小间隔,平衡了栅栏现象,使我们能够得到更多的信息。
2.L=32,64,128,256,512谱峰对应的频率依次为0.1875π,0.1875π,0.2031π,0.2031π,
0.1992π;而谱峰对应的频率的理论值为0.2π。
L=32,64,128,256,512时对应的误差依次为
0.0125π,0.0125π,0.0031π,0.0031π,0.0008π。
说明L的数值越大、频谱的间隔越小,得到的信息就越准确。
【自主学习内容】
Find函数的使
【仿真程序】
k=0:
31;
xk=sin(0.2*pi*k);
L1=32;
Y1=fft(xk);
figure
(1);
stem(k/length(Y1),abs(Y1),'.');
title('L=32');xlabel('NormalizedFrequency');ylabel('Amplitude');
L2=64;L3=128;L4=256;L5=512;
figure
(2);
Y2=fft(xk,L2);
k2=0:
length(Y2)-1;
subplot(2,2,1);stem(k2/length(Y2),abs(Y2),'.');
title('L=64');xlabel('NormalizedFrequency');ylabel('Amplitude');
Y3=fft(xk,L3);
k3=0:
length(Y3)-1;
subplot(2,2,2);stem(k3/length(Y3),abs(Y3),'.');
title('L=128');xlabel('NormalizedFrequency');ylabel('Amplitude');
Y4=fft(xk,L4);
k4=0:
length(Y4)-1;
subplot(2,2,3);stem(k4/length(Y4),abs(Y4),'.');
title('L=256');xlabel('NormalizedFrequency');ylabel('Amplitude');
Y5=fft(xk,L5);
k5=0:
length(Y5)-1;
subplot(2,2,4);stem(k5/length(Y5),abs(Y5),'.');
title('L=512');xlabel('NormalizedFrequency');ylabel('Amplitude');
%以下程序是为求频谱中谱峰的频率
Y1max=max(abs(Y1));%求出谱峰的幅度
x1=find(abs(Y1)==Y1max)-1;%找到谱峰对应的点
m1=x1*(2/L1);%找到谱峰对应的频率,为便于观察,输出的结果应乘以pi
disp('L=32时幅度最大值');disp(Y1max);%输出谱峰的幅度
disp('L=32时谱峰对应的频率(单位:
πHz)');disp(m1);%输出谱峰对应的频率
Y2max=max(abs(Y2));
x2=find(abs(Y2)==Y2max)-1;
m2=x2*(2/L2);
disp('L=64时幅度最大值');disp(Y2max);
disp('L=64时谱峰对应的频率(单位:
πHz)');disp(m2);
Y3max=max(abs(Y3));
x3=find(abs(Y3)==Y3max)-1;
m3=x3*(2/L3);
disp('L=128时幅度最大值');disp(Y3max);
disp('L=128时谱峰对应的频率(单位:
πHz)');disp(m3);
Y4max=max(abs(Y4));
x4=find(abs(Y4)==Y4max)-1;
m4=x4*(2/L4);
disp('L=256时幅度最大值');disp(Y4max);
disp('L=256时谱峰对应的频率(单位:
πHz)');disp(m4);
Y5max=max(abs(Y5));
x5=find(abs(Y5)==Y5max)-1;
m5=x5*(2/L5);
disp('L=512时幅度最大值');disp(Y5max);
disp('L=512时谱峰对应的频率(单位:
πHz)');disp(m5);
【研讨题目】基本题
4.已知一离散序列为x[k]=AcosΩ0k+Bcos((Ω0+∆Ω)k)。
用长度N=64的哈明窗对信号截短后近似计算其频谱。
试用不同的A和B的取值,确定用哈明窗能分辩的最小的谱峰间隔
中c的值。
(M2-3)
【仿真结果】
1.A=B=1,c从0.5到6变化时分辨率的变化非常生动,我将其做成了一个视频(Hamming.mpg)。
部分仿真结果如下:
(注意:
标题为c的值)
2.A=1,B=0.2时,仿真结果如下:
(注意:
标题为c的取值)
【结果分析】
C的临界值在2左右,收到幅度等其他因素的影响会有所不同
【自主学习内容】
循环语句的使用
【问题探究】
在离散序列频谱计算中为何要用窗函数?
用不同的窗函数对计算结果有何影响?
与矩形窗相比哈明窗有何特点?
如何选择窗函数?
为了将无限长序列截短,使之能够进行DFT;窗函数的截止方式对频谱的泄露现象有巨大影响,进而直接影响到分辨率;根据实际情况合理选择窗函数,如,考虑旁瓣是否会将另一谱峰掩盖掉。
【仿真程序】
N=64;
k=0:
N-1;
A=1;B=1;
W=10;
wh=(hamming(N))';
w=0:
1023;
i=6;n=57;
whilei>=0.5
dW=2*pi*i/N;
x=A*cos(W*k)+B*cos((W+dW)*k);
Y=fft(x.*wh,1024);
n=n-1;
figure(n);plot(w/length(Y),abs(Y)/max(abs(Y)));gridon;
title(i)
xlabel('NormalizedFrequency');ylabel('NormalizedAmplitude');
holdon;
i=i-0.1;
end
【研讨题目】基本题
5.已知一离散序列为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)
【题目分析】
分析影响谱峰分辨率的主要因数,进一步认识补零在在频谱计算中的作用。
【仿真结果】
(3)
【结果分析】
(2)不能通过添加0的方法使分辨率提高,因为在抽取的时候是抽取64个点,此时造成的信息损失是不可逆的,不可以通过添0的办法增加死去的信息。
(3)要想提高分辨率,就要提高抽取的点数,所以可以通过加窗函数(窗的长度大于64)的办法使频谱分开。
【仿真程序】
(1)k=0:
64;omega0=0.4*pi;omega1=0.4*pi+pi/64;
x=cos(omega0*k)+0.75.*cos(omega1*k);
k=0:
64;omega0=0.4*pi;omega1=0.4*pi+pi/64;
x=cos(omega0*k)+0.75.*cos(omega1*k);
Hm=fft(x,64);
w=linspace(-pi,pi,64);
stem(w,abs(Hm),'.');xlabel('频率');ylabel('幅度');title('L=64')
(2)k=0:
64;omega0=0.4*pi;omega1=0.4*pi+pi/64;
x=cos(omega0*k)+0.75.*cos(omega1*k);
k=0:
64;omega0=0.4*pi;omega1=0.4*pi+pi/64;
x=cos(omega0*k)+0.75.*cos(omega1*k);
Hm=fft(x,128);
w=linspace(-pi,pi,128);
stem(w,abs(Hm),'.');xlabel('频率');ylabel('幅度');title('L=128)
k=0:
64;omega0=0.4*pi;omega1=0.4*pi+pi/64;
x=cos(omega0*k)+0.75.*cos(omega1*k);
k=0:
64;omega0=0.4*pi;omega1=0.4*pi+pi/64;
x=cos(omega0*k)+0.75.*cos(omega1*k);
Hm=fft(x,256);
w=linspace(-pi,pi,256);
stem(w,abs(Hm),'.');xlabel('频率');ylabel('幅度');title('L=256');
(3)
k=0:
128;L=128;
w0=0.4*pi;w1=w0+pi/64;
x=cos(w0*k)+0.75*cos(w1*k);
X=fftshift(fft(x,L));
w=(-2*pi/2+(0:
L-1)*2*pi/L)/(2*pi);
plot(w,abs(X));xlabel('Normalizedfrequency');ylabel('Magnitude');title('128点');
【研讨内容】—
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 第二次 doc