数字信号处理project2

数字信号处理Project21.实验题目Designalinear-phaseFIRlowpassfilterusingthewindowdesigntechnique.Thespecificationsaregivenasfollows:stopbandedge:5rad/secpassbandedge:3rad/secmaximumpassbandattenuation:0.3dBminimumstopbandbandattenuation:40dBsamplingfrequency:15rad/sec(a).Useeachofthefollowingwindowsforthedesign:Hamming,Hann,Blackman,andKaiser.Showtheidealimpulseresponse,theactualimpulseresponse,andplotthegainresponseofthedesignedfiltersforeachcase.Givethemaindesignsteps.Commentonyourresults.Donotusethefunctionfir1ofMatlab.(b).Repeat(a)exceptfortheidealimpulseresponse,byusingthefunctionfir1ofMatlab.Compareresultsof(a)withthoseof(b).2.实验设计(a)题：1.确定通带截止频率由于题中给定的频率为实际频率，所以先对给定频率进行归一化。1532p，1552s通带的截止频率为：2psc2.确定不同窗下的长度M，N归一化的过渡带宽为ps对于Hann、Hamming、Blackman窗用下面的函数计算M，N。c通过查表得,分别为3.11,3.32,5.56。cM12MN对于Kaiser窗，M、N用下面的公式计算，其中s为最小阻带衰减为：40dB。21NMΔ285.28sN3.计算理想低通的脉冲响应MnMnnnhdc,sin][4.使用窗函数产生窗系数winHamming:)cos(46.054.0][Mnnw，MnMHann:)cos(5.05.0][Mnnw，MnMBlackman:)2cos(08.0)cos(5.042.0][MnMnnw，MnMKaiser:因为4050，所以，)21(07886.0)21(5842.04.0ss)()/(1][020IMnInw，MnM5.对理想低通滤波器的脉冲响应加窗'*.winhdht6.画出频谱图(b)题：将(a)中计算得到的窗函数win，直接带入matlab自带的函数fir1，即可产生时域的加窗低通滤波的冲激响应：filter_t=fir1(N-1,wc/pi,win);3.实验结果与分析(a)题：Hamming:Hann:Blackman:Kaiser:分析：图中标出了spww与。窗的类型最大通带衰减(dB)最小阻带衰减(dB)Hamming-0.006-73.1Hann-0.036-50.58Blackman-0.002-76.3Kaiser-0.176-36.34代码运行结果：由此可以看出，Kaiser窗比较接近指标，同时其阶数（N=17）也是最少的。但就衰减性能而言，其余三个的最小阻带衰减(-73.1dB，-50.58dB，-76.3dB)比Kaiser的-36.34dB要好。Hamming窗阻带波纹较大，Blackman窗阻带衰减较好但波纹较密，Hann窗处于二者之间,Kaiser窗较稳定，同时波纹较大。因为Hamming窗阻带衰减较好，阶数不算太大，阻带波纹较大，所以可以用来设计要求不高的FIR滤波器；因为Blackman阻带衰减较好,且随着频率的增加，波纹衰减增强，但因为阶数比其他要大，所以可以用来设计要求较高且较少考虑成本和复杂度的FIR滤波器；Kaiser窗有最低的阶数，所以可以用来设计较精确的阻带衰减滤波器。(b)题：Hamming:Hann:Blackman:Kaiser:代码运行结果：分析：各窗在两种方式下的图形相似，阶数相同，由此可见结果是一致的。4.代码(a)题：wp=2*pi*3/15;ws=2*pi*5/15;alphap=0.3;%最大通带衰减alphas=40;%最小阻带衰减NUM_Hamming=1;NUM_Hann=2;NUM_Blackman=3;NUM_Kaiser=4;c=[3.32,3.11,5.56]*pi;%各个的窗函数的C值wc=(ws+wp)/2;delt_w=ws-wp;freq_labels={'hamming窗低通滤波器增益响应','hann窗低通滤波器增益响应','Blackman窗低通滤波器增益响应','Kaiser窗低通滤波器增益响应'};ht_labels={'hamming窗真实冲激响应','hann窗真实冲激响应','Blackman窗真实冲激响应','Kaiser窗真实冲激响应'};hd_labels={'理想冲激响应','理想冲激响应','理想冲激响应','理想冲激响应'};forfilter_kind=NUM_Hamming:NUM_Kaiser,switchfilter_kind,caseNUM_Hamming,M=ceil(c(filter_kind)/delt_w);%向上取整N=2*M+1;win=hamming(N);display(['Hamming窗生成的冲激响应系数：','(阶数N=',num2str(N),')']);figurecaseNUM_Hann,M=ceil(c(filter_kind)/delt_w);%向上取整N=2*M+1;win=hann(N);display(['Hann窗生成的冲激响应系数：','(阶数N=',num2str(N),')']);figurecaseNUM_Blackman,M=ceil(c(filter_kind)/delt_w);%向上取整N=2*M+1;win=blackman(N);display(['Blackman窗生成的冲激响应系数：','(阶数N=',num2str(N),')']);figurecaseNUM_Kaiser,N=ceil((alphas-8)/(delt_w*2.285));%向上取整M=(N-1)/2;b=0.5842*(alphas-21)^0.4+0.07886*(alphas-21);win=kaiser(N,b);display(['Kaiser窗生成的冲激响应系数：','(阶数N=',num2str(N),')']);figureotherwisedisp('error');endn=-M:M;hd=sin(wc*n)./(pi*n);hd(find(n==0))=wc*cos(wc*0)/pi;%因为n=0时公式“sin(wc*n)./(pi*n)”中分母为0，计算结果为NaN,所以应该另用%一种方法单独计算n=0时的hd值。这里采用洛必达法则对原hd公式进行分子分母微分%得hd=wc*cos(w%c*n)./pi，n-0;这样可以求得当n=0时原函数值。(也可以将n加上eps%—-系统精度，使n最大程度接近0而不等0）ht=hd.*win';display(['',num2str(ht)]);%打印显示冲激响应系数%wvtool(ht);subplot(1,3,1);plot(n,hd,'.-')%画理想冲激响应系数图title(hd_labels(filter_kind));xlabel('n','FontSize',12);ylabel('hd','fontsize',12);gridonsubplot(1,3,2);plot(n,ht,'.-')%画冲激响应系数图title(ht_labels(filter_kind));xlabel('n','FontSize',12);ylabel('ht','fontsize',12);gridon[h,w]=freqz(ht,1,512);W=w/pi;H=20*log10(abs(h));subplot(1,3,3);%画增益响应holdontitle(freq_labels(filter_kind));plot(W,H);xlabel(['pi(omega/omega','s)']);ylabel('增益(dB)');gridon%为了更直观的看到所设计FIR是否达到指标，可分别在图中标出相关点%从而可以对不同窗函数进行对比。标示过程如下：dotp=round(mean(find(wwp-0.1&wwp+0.1)));%找到点wp的坐标dots=round(mean(find(wws-0.1&wws+0.1)));%找到点ws的坐标plot(W(dotp),H(dotp),'.r','MarkerSize',25);%对应处打上显目的红点plot(W(dots),H(dots),'.r','MarkerSize',25);text(W(dotp),H(dotp),['[','',num2str(round(W(dotp)*1000)/1000),...',',num2str(round(H(dotp)*1000)/1000),'dB]'],'FontSize',15);%标示wp在频响图%中的位置，从中可以%看来相应增益text(W(dots),H(dots),['[','',num2str(round(W(dots)*100)/100),...',',num2str(round(H(dots)*100)/100),'dB]'],'fontsize',15);%标示ws在频响图%的位置holdoffend(b)题：wp=2*pi*3/15;ws=2*pi*5/15;alphap=0.3;%最大通带衰减alphas=40;%最小阻带衰减NUM_Hamming=1;NUM_Hann=2;NUM_Blackman=3;NUM_Kaiser=4;c=[3.32,3.11,5.56]*pi;%各个的窗函数的C值wc=(ws+wp)/2;delt_w=ws-wp;freq_labels={'hamming窗低通滤波器增益响应','hann窗低通滤波器增益响应','Blackman窗低通滤波器增益响应','Kaiser窗低通滤波器增益响应'};ht_labels={'hamming窗真实冲激响应','hann窗真实冲激响应','Blackman窗真实冲激响应','Kaiser窗真实冲激响应'};forfilter_kind=NUM_Hamming:NUM_Kaiser,switchfilter_kind,caseNUM_Hamming,M=ceil(c(filter_kind)/delt_w);%向上取整N=2*M+1;win=hamming(N);display(['hamming窗生成的冲激响应系数：','(阶数N=',num2str(N),')']);figurecaseNUM_Hann,M=ceil(c(filter_kind)/delt_w);%向上取整N=2*M+1;win=hann(N);display(['hann窗生成的冲激响应系数：','(阶数N=',num2str(N),')']);figurecaseNUM_Blackman,M=ceil(c(filter_kind)/delt_w);%向上取整N=2*M+1;win=blackman(N);display(['Blackman窗生成的冲激响应系数：','(阶数N=',num2str(N),')']);figurecaseNUM_Kaiser,N=ceil((alphas-8)/(delt_w*2.285));