数字信号处理算法课程设计

一、目的与要求掌握《数字信号处理基础》课程的基本理论；掌握应用MATLAB进行数字信号处理的程序设计方法。二、主要内容设计题目及设计要求：已知低通数字滤波器的性能指标如下：，，，要求：选择合适的窗函数，设计满足上述指标的数字线性相位FIR低通滤波器。用一个图形窗口，包括四个子图，分析显示滤波器的单位冲激响应、相频响应、幅频响应和以dB为纵坐标的幅频响应曲线。用双线性变换法，设计满足上述指标的数字ChebyshevI型低通滤波器。用一个图形窗口，包括三个子图，分析显示滤波器的幅频响应、以dB为纵坐标的幅频响应和相频响应。已知模拟信号其中，，，，取采样频率。要求：以对进行取样，得到。用一个图形窗口，包括两个子图，分别显示以及()的波形；用FFT对进行谱分析，要求频率分辨率不超过5Hz。求出一个记录长度中的最少点数，并用一个图形窗口，包括两个子图，分别显示以及的幅值；用要求1中设计的线性相位低通数字滤波器对进行滤波，求出滤波器的输出，并用FFT对进行谱分析，要求频率分辨率不超过5Hz。求出一个记录长度中的最少点数，并用一个图形窗口，包括四个子图，分别显示（）、、（）和的幅值；用要求2中设计的Chebyshev低通数字滤波器对进行滤波，求出滤波器的输出，并用FFT对进行谱分析，要求频率分辨率不超过5Hz。求出一个记录长度中的最少点数，并用一个图形窗口，包括四个子图，分别显示（）、、（）和的幅值。备注：(1)要求编写一个主程序，完成上述问题。(2)要求各结果图均标出图题以及横纵轴的名称。(3)要求给主要语句加上标注。三、设计思路1）窗函数设计数字线性相位FIR低通滤波器1、根据设计给出的技术指标，计算允许的过渡带宽度及阻带衰减，初步选择窗函数和N值(可选用hamming窗，Blackman窗和Kaiser窗,这里选用Hamming窗)。2、若选用理想低通逼近，则计算积分，求出理想低通的冲激响应hd(n);理想低通的截止频率选择为wc=(ws+wp)/2,hd(n)=sin[wc(n-a)]/pi*n-a);3、将hd(n)与窗函数相乘的FIR数字滤波器的冲激响应h(n)：h(n)=hd(n)w(n)4、计算FIR数字滤波器的频率响应，并验证是否达到所要求的指标。2）双线性变换法设计的数字ChebyshevI型低通滤波器1、根据设计给出的数字滤波器的技术指标，计算模拟原型低通的指标2、模拟切比雪夫1型滤波器的计算，利用function[b,a]=afd_cheb1(wp,ws,Rp,As)进行切比雪夫1型模拟滤波器的设计；利用function[b,a]=cheb1ap_o(N,Rp,omega)进行非归一化的切比雪夫1型模拟原型低通滤波器的设计3、双线性变换法将Ha(s)转换为H(z)4、直接型转换为级联型，计算数字滤波器的频率响应，并验证是否达到所要求的指标。3）模拟信号的波形产生和采样及FFT频谱分析和FIR滤波及IIR滤波N：一个记录长度中的取样数F：频率分辨率fs：采样周期最小记录长度的选择：由tp=NT=1/F可知，给定抽样频率和频率分辨率后，就可以得出最小记录长度N=fs/Ffft(x)：计算N点的DFT,N是序列的长度fft(x，L)：计算L点的DFT,若N>L，则将原序列截短为L点的序列，再计算其L点的DFT；若N<L，则将原序列补零为L点的序列，再计算其L点的DFT；ifft(x)：计算N点的IDFT,N是频域序列的长度ifft(x，L)：计算L点的IDFT,若N>L，则将原序列截短为L点的序列，再计算其L点的IDFT；若N<L，则将原序列补零为L点的序列，再计算其L点的IDFT；为了提高fft和ifft的计算效率，应尽量使序列长度N为2的幂次方，或将序列补零使其长度为2的幂次方。filter(b,a,x)：可采用此函数给出滤波器的输出b=H(z)的分子多项式系数（对FIRb=h)a=H(z)的分母多项式系数（对FIRa=[1])三、程序代码1、主程序代码%一、数字线性相位FIR低通滤波器wp=0.26*pi;ws=0.41*pi;As=50;%数字滤波器的指标tr_width=ws-wp;%(过渡)带宽N=ceil(6.6*pi/tr_width)+1;%经验公式,求采样点数n=[0:1:N-1];wc=(ws+wp)/2;%理想低通的截止频率hd=ideal_lp(wc,N);%理想低通的冲激响应w_ham=(hamming(N))';%哈明窗的窗函数h=hd.*w_ham;%FIR滤波器的实际冲激响应[db,mag,pha,grd,w]=freqz_m(h,[1]);%计算数字滤波器的频率响应delta_w=2*pi/1000;%取1000点Rp=-(min(db(1:1:wp/delta_w+1)));%实际的通带衰减,间隔1dbAs=-round(max(db(ws/delta_w+1:1:501)));%实际的最小阻带衰减figure(1);%分别显示滤波器的单位冲激响应、相频响应、幅频响应和以dB为纵坐标的幅频响应曲线。subplot(2,2,1);stem(n,hd);ylabel('hd(n)');xlabel('n');title('(理想)单位冲激响应');grid;axis([045-0.10.5]);subplot(2,2,2);plot(w/pi,pha);ylabel('相位');xlabel('以\pi为单位的频率');title('相位响应');grid;axis([01-44]);subplot(2,2,3);plot(w/pi,mag);ylabel('|幅度|');xlabel('以\pi为单位的频率');title('幅度响应');grid;axis([01.5-0.51.5]);subplot(2,2,4);plot(w/pi,db);ylabel('对数幅度/db');xlabel('以\pi为单位的频率');title('幅度响应(db)');grid;axis([01.5-1600]);%二、双线性变换法设计数字切比雪夫I型低通滤波器%数字滤波器指标wp=0.26*pi;ws=0.41*pi;Rp=0.75;As=50;%转换为模拟滤波器的指标（模拟原型低通指标）T=1;Fs=1/T;omegap=(2/T)*tan(wp/2);omegas=(2/T)*tan(ws/2);[cs,ds]=afd_cheb1(omegap,omegas,Rp,As);%模拟切比雪夫I型低通滤波器的计算,求解拉普拉斯变换的分子与分母多项式的系数数组[b,a]=bilinear(cs,ds,T);%双线性变换法[C,B,A]=dir2cas(b,a)%IIRDF直接型到级联型的转换函数[db,mag,pha,grd,w]=freqz_m(b,a);%计算数字滤波器的频率响应,频域上采样并记录幅频和相频响应%分别显示数字切比雪夫I型低通滤波器的幅频响应、以dB为纵坐标的幅频响应和相频响应figure(2)subplot(3,1,1);plot(w/pi,mag);ylabel('|幅度|');xlabel('以\pi为单位的频率');title('幅度响应');grid;axis([00.8-0.51.1]);subplot(3,1,2);plot(w/pi,db);title('以db为单位的幅度响应');grid;ylabel('对数幅度/db');xlabel('以\pi为单位的频率');axis([0,1.5-40010]);subplot(3,1,3);plot(w/pi,pha);ylabel('相位');xlabel('以\pi为单位的频率');title('相位响应');grid;axis([01.2-44]);%三（1）、以对xt进行抽样，得到Xnn=0:1:511;%做512点DFTfs=10000;%采样频率t=0:0.000001:10;xt=2*sin(2*pi*120*t)+5*sin(2*pi*4980*t)+8*cos(2*pi*3250*t)+7.5*cos(2*pi*1150*t);T=1/fs;%采样周期t0=n*T;%时间离散化Xn=2*sin(240*pi*t0)+5*sin(2*4980*pi*t0)+8*cos(2*pi*3250*t0)+7.5*cos(2*pi*1150*t0);figure(3)%显示x(n)(0<=n<=511)以及x(t)的波形；subplot(2,1,1);stem(n,Xn);title('Xn');grid;ylabel('Xn');xlabel('n');axis([01000-2525])subplot(2,1,2);plot(t,xt);title('xt曲线');grid;ylabel('xt');xlabel('t');axis([00.01-2525]);%三（2）、对Xn（0<=n<=Nx-1）进行FFT变换求Xk，对Xn进行谱分析F=5;Nx=fs/F;%取样数(因为此时fs远大于x(t)中的最高频率,应该用fs才能分辨出四个谱峰f1~f4),计算得Nx=2000点T=1/fs;%采样周期n=0:Nx-1;%做2000点DFTt0=n*T;%时间离散化Xn=2*sin(240*pi*t0)+5*sin(2*4980*pi*t0)+8*cos(2*pi*3250*t0)+7.5*cos(2*pi*1150*t0);%fprintf('\n***直接FFT时一个记录长度最小点数Nx=%2.0f\n',Nx);k=0:Nx-1;Xk=fft(Xn);magX=abs(Xk);figure(4)%显示x(n)以及X(k)的幅值；subplot(2,1,1);stem(n,Xn);title('Xn');grid;ylabel('Xn');xlabel('n');axis([02000-2525]);subplot(2,1,2);stem(k,magX);title('Xk');grid;ylabel('magX');xlabel('k');axis([02048010000]);%三（3）、用线性相位低通滤波器对Xn进行滤波得到y1n，对y1n进行FFT变换Fy1=5;%频率分辨率f0=1150;Ny1=2*f0/Fy1;%取样数(x(n)通过低通滤波器之后，高频成分f2,f3被滤除,仅剩下f1,f4,计算得Ny1=2*f4/Fy1=460点)%fprintf('\n**y1n用FFT进行谱分析时一个记录长度最小点数Ny1=%2.0f\n',Ny1);y1n=filter(h,[1],Xn);%利用递归滤波器(IIR)或非递归滤波器(FIR)对数据进行滤波(对FIRb=h,a=[1]),Xn经过用哈明窗设计的线性相位低通滤波器滤波得到y1nYk1=fft(y1n);%对y1n进行FFT变换magYk1=abs(Yk1);figure(5)%显示x(n)(0<=n<=Nx-1)、X(k)、y1(n)(0<=n<=Ny1-1)和Y1(k)的幅值；subplot(2,2,1);stem(n,Xn);title('Xn');grid;ylabel('Xn');xlabel('n');axis([02000-2525]);subplot(2,2,2);stem(k,magX);title('Xk');grid;ylabel('magX');xlabel('k');axis([02048010000]);subplot(2,2,3);stem(n,y1n);title('y1n');grid;ylabel('y1n');xlabel('n');axis([0200-2020]);subplot(2,2,4);stem(k,magYk1);title('Yk1');grid;ylabel('magY1');xlabel('k');axis([0204808000]);%三（4）、用由双线性变换法设计的数字切比雪夫I型低通滤波器对Xn进行滤波得到y2n，对y2n进行FFT变换y2n=filter(b,a,Xn);%利用递归滤波器(IIR)或非递归滤波器(FIR)对数据进行滤波(b和a值在要求2中已求出，直接调用)%fprintf('\n***y2n用FFT进行谱分析时一个记录长度最小点数Ny2=%2.0f\n',Ny2);Ny2=Ny1;%%取样数(x(n)通过低通滤波器之后，高频成分f2,f3被滤除,仅剩下f1,f4,计算得Ny2=Ny1=460点)XK2=fft(Xn);%对Xn进行DFT运算AM1=abs(XK2);%取模值|XK|YK2=fft(y2n);%对y2n进行DFT运算AM2=abs(YK2);%取模值|YK|figure(6)%显示x(n)(0<=n<=Nx-1)、X(k)、y2(n)(0<=n<=Ny2-1)和Y2(k)的幅值subplot(2,2,1);stem(n,Xn);grid;xlabel('n');ylabel('x(n)');title('取样信号x(n)');axis([0100-4040]);subplot(2,2,2);stem(n,AM1);grid;xlabel('k');ylabel('X(k)');title('X(k)的幅值');subplot(2,2,3);stem(Ny2,y2n);grid;xlabel('n');ylabel('y2(n)');title('滤波器的输出y2(n)');axis([0100-1010]);subplot(2,2,4);stem(Ny2,AM2);grid;xlabel('k');ylabel('Y2(k)');title('Y2(k)的幅值');输出要求的参数结果如下：***FIR低通滤波器的阶数=46As=53***切比雪夫滤波器的阶