数字信号处理第二次研学.doc

数字信号处理课程研究性学习报告姓名 学号 同组成员 指导教师 时间 DFT近似计算信号频谱专题研讨【目的】(1) 掌握利用DFT近似计算不同类型信号频谱的原理和方法。(2) 理解误差产生的原因及减小误差的方法。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】 基本题 1利用DFT分析x(t)=Acos(2pf1t)+Bcos(2pf2t)的频谱，其中f1=100Hz，f2=120Hz。(1)A=B=1; (2)A=1,B=0.2。【题目分析】分析题目，给出合适的DFT参数由取样定理知，要使信号频谱不混叠，则抽样频率不小于最高频率的两倍。而要满足信号分辨率的要求，抽样点数Nfsam/f。在对信号做DFT时，由于对信号进行截短，因此会产生频谱泄漏，要想从频谱中很好的分辨出个频率分量，需要考虑时域抽样频率，所加的窗函数，窗函数的长度，以及DFT的点数等参数对结果的影响。(1)A=B=1，即x(t)=cos(2pf1t)+cos(2pf2t)矩形窗1：条件：fsam=240Hz；N=20；L=512矩形窗2：条件：fsam=600Hz；N=40；L=512矩形窗3：fsam=1200Hz；N=80；L=512Hamming窗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(2pf1t)+0.2cos(2pf2t)矩形窗：N=100；L=512；fs=600Hamming窗：N=100；L=512；fs=600【仿真结果】【结果分析】对实验结果进行分析比较，回答：加窗对谱分析有何影响？如何选择合适的窗函数？选择合适DFT参数的原则？在（1）中进行矩形窗仿真时，我们选择了不同的fsam，分别为240,600,1200它们均满足抽样定理，但是我们在实验中却发现，在240hz时出现了混叠现象。所以，在实际应用中抽样频率应大于最低抽样频率3-5倍才能有更好的结果。进行hamming窗仿真时，在保证抽样频率相同的条件下，取不同的长度也40,120。其中，N=40不满足N=60的要求，我们可以看到出现了混叠，利用Nfsam/f，我们可以算出当fs=600时N为60时恰好可以分出频谱，而实际中N=60时无法分出两个频率分量，而N=120时，仿真效果良好在（2）的条件下进行仿真时，我们选取了相同的N、L、fsam值，但是分别使用了矩形窗和哈明窗。使用矩形窗时，幅度较小的峰值与旁瓣的幅度接近，甚至难以区分，效果不理想。使用哈明窗后，泄露现象被有效遏制，所以可以清楚区分主瓣、旁瓣。所以，在选择参数进行DFT变换时，应该保证抽样频率满足抽样定理，并且能大于最小抽样值3-5倍。长度选择保证Nfsam/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(jomega);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(jomega);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(jomega);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(jomega);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时的DFTL=32、64、128、256、512时用 DFT求出的频谱中谱峰对应的频率（注意：谱峰对应的频率的单位是Hz）【结果分析】1. zero padding对我们“看到的”频谱是有影响的，但对频谱本身是没有影响。它没有改变频谱和序列一一对应的关系，但减小了频谱的最小间隔，平衡了栅栏现象，使我们能够得到更多的信息。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(Normalized Frequency);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(Normalized Frequency);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(Normalized Frequency);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(Normalized Frequency);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(Normalized Frequency);ylabel(Amplitude);%以下程序是为求频谱中谱峰的频率Y1max=max(abs(Y1);%求出谱峰的幅度x1=find(abs(Y1)=Y1max)-1;%找到谱峰对应的点m1=x1*(2/L1);%找到谱峰对应的频率，为便于观察，输出的结果应乘以pidisp(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=AcosW0k+Bcos ( (W0+DW)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;while i=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);grid on; title(i) xlabel(Normalized Frequency);ylabel(Normalized Amplitude); hold on; i=i-0.1;end【研讨题目】 基本题 5已知一离散序列为 xk=cos(W0k)+0.75cos(W1k), 0 k 63 其中W0=0.4p, W1=W0+p/64(1) 对xk做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(Normalized frequency);ylabel(Magnitude);title(128点);【研讨内容】基本题6. 语音信号的频率分析(1) 采集若干wav格式的男女生话音信号、女高音男低音演唱信号。(2) 分析所采集信号的频谱分析，给出男生和女生话音信号的频率范围，女高音的最高频率是多少，男低音的最低频率是多少？【题目分析】题目要求采集wav信号。主要用Wavread函数分析【仿真结果】男生朗诵：女生朗诵：男生清唱：女生清唱：【结果分析】在朗诵时男女生的频率范围差不多主要在1100hz左右，但在清唱是女生频率要高于男生，女高音最高可达20000，男低音可达100hz【自主学习内容】音频函数的读取【发现问题】fft转换后要进行shiftfft【仿真程序】y,fs,nbits=wavread (C:UsersAdministratorDesktop男生朗诵.wav);sound(y,fs,nbits);n = length (y) ; Y=fftshift(fft(y,n); subplot(2,1,1);plot(y);title(原始信号波形);subplot(2,1,2);plot(abs(Y);title(原始信号频谱)【研讨题目】 扩展题7本题研究连续周期信号频谱的近似计算问题。 周期为T0的连续时间周期信号x(t)可用Fourier级数表示为其中X(nw0)称为连续时间周期信号x(t)的频谱函数。称为信号的基频（基波），称为信号的谐波。如果信号x(t)函数表达式已知，则可由积分得出信号的频谱。如果信号x(t)函数表达式未知，或者x(t)函数表达式非常复杂，则很难由积分得信号的频谱。本题的目的就是研究如何利用DFT近似计算连续时间周期信号的频谱。(1)若在信号x(t)的一个周期T0内抽样N个点，即， T为抽样周期(间隔)，可获得序列xk试分析序列xk的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)若要用10次以内的谐波近似表示x(t)，试给出计算方案，并计算出近似表示的误差。讨论出现误差的原因及减小误差的方法。 【理论推导】DFT计算所得结果Xm与连续周期信号频谱X(nw0)的关系。【计算方案】根据理论推导结果设计近似计算方案。分析产生误差的主要原因。Xm是连续周期信号频谱X(n0)的周期化，在乘上比例因子N所得，由于对于一般连续周期信号的频谱，当n时，X（n0）0,所以可以用Xm来近似较为低次的谐波分量，但仍会有一些高频谐波分量在周期化是发生混叠。如果信号是带限信号，但抽样时不满足抽样定理，也同样会产生混叠误差。用矩形窗截取了一个周期长度上的信号抽样值，因此会产生频率泄露。如果信号不是带限信号，在分析其频谱的时候，要截断成时限的谱，此时会产生截断误差。在对信号