数字信号处理课设报告.doc

长安大学数字信号处理课设报告 学 院：信息工程学院 专 业：电子信息工程 班 级：2011240301 姓 名： 学 号：22 指导老师：陈玲目录一、课程设计目的-1二、课程设计原理-11. IIR滤波器设计原理-12.用窗函数法设计FIR滤波器原理-53.用FFT对信号作频分析的原理-6三、课程设计过程与结果-71.IIR（无限脉冲响应）模拟滤波器设计-72. IIR（无限脉冲响应）数字滤波器设计-83. FIR（有限脉冲响应）数字滤波器设计-114.利用FFT进行频谱分析-14四、课设心得-15五、参考文献-17一、课程设计目的1.使学生增进对MATLAB的认识，学会MATLAB的使用，掌握MATLAB的程序设计方法，并加深对数字信号处理理论方面的理解。2.使学生掌握数字信号处理中IIR和FIR滤波器的设计。3.使学生了解和掌握用MATLAB实现IIR和FIR滤波器的设计方法、过程，为以后的设计打下良好基础。二、课程设计原理1. IIR滤波器设计原理：滤波是信号处理的一种最基本而重要的技术。利用滤波从，复杂的信号中提取所需要的信号，抑制不需要的部分。滤波器是具有一定传输特性的信号处理装置。数字滤波器的工作原理：数字滤波器是具有一定传输特性的数字信号处理装置。它的输入和输出均为离散的数字信号，借助数字器件或一定的数值计算方法，对输入信号进行处理，改变输入信号的波形或频谱，达到保留信号中有用成分去 除无用成分的目的。如果加上A/D、D/A转换，则可以用于处理模拟信号。滤波器的分类：滤波器的种类很多，分类方法也不同。按处理的信号划分：模拟滤波器、数字滤波器；按频域特性划分；低通滤波器、高通滤波器、带通滤波器、带阻滤波器；按时域特性划分：FIR滤波器、IIR滤波器。滤波器设计步骤：（1） 按任务要求确定Filter的性能指标；（2） 用IIR或FIR系统函数去逼近这一性能要求；（3） 选择适当的运算结构实现这个系统函数； （4） 软件还是用硬件实现。IIR滤波器设计：由于它的脉冲响应序列h(n)是无限长的，称为无限长冲激响应滤波器。IIR滤波器的设计根据滤波器某些性能指标要求，设计滤波器的分子和分母多项式。设计方法：（1） 模拟滤波器变换法（经典设计法）；（2） 直接设计法（3） 参数模型设计法（4） 最大平滑滤波器设计法IIR设计方法比较：（1）借助模拟filter的设计方法（经典设计法）a. 对设计性能指标中的频率指标进行转换使其满足模拟滤波器原型设计性能指标；b. 估计模拟滤波器最小阶数和边界频率。Matlab提供的函数buttord，cheb1ord，cheb2ord，ellipord；c. 设计模拟低通滤波器原型。Matlab提供的函数buttap，cheb1ap，cheb2ap,ellipap；d. 由模拟低通原型经频率变换获得模拟滤波器。Matlab提供的函数lp2lp，lp2hp，lp2bp，lp2bs；e. 将模拟滤波器离散化获得IIR数字滤波器。Matlab提供的函数bilinear，impinvar。IIR滤波器完全设计函数：在MATLAB信号处理工具箱中提供了IIR滤波器设计的完全工具函数，用户只要调用这些工具函数即可一次性完成设计，而不需像上面分步实现。MATLAB提供的函数有：Butter、cheby1、cheby2、ellip等。在使用这些函数设计数字滤波器时，数字频率采用标准化频率（归一化频率）。归一化频率为：频率的取值范围在01之间，标准化频率1对应的数字频率为对应的模拟频率为采样频率的一半。归一化的处理方法：归一化频率=模拟频率/采样频率一半。巴特沃斯数字滤波器：格式：b,a=butter(n,wn,ftype)。其中，n为滤波器阶数；wn为滤波器截止频率；ftype为滤波器类型：high为高通滤波器，截止频率wn ，stop为带阻滤波器，截止频率wn=w1,w2；缺省时为低通或带通滤波器；b,a分别为滤波器传递函数分子和分母的系数向量。（2） IIR滤波器直接设计法经典设计法只限于几种标准的低通、高通、带通、带阻滤波器，而对于具有任意形状或多频带滤波器的设计是无能为力的。直接设计法采用最小二乘法拟合给定的幅频响应，使设计的滤波 器幅频特性逼近期望的频率特性。MATLAB提供的工具函数。函数调用格式：b,a=yulewalk(n,f,m)说明：n：滤波器阶数；f：给定的频率点向量（标准频率），第一个频点必须为0，最后一个必须为1；m：和频率向量对应的理想幅值响应向量；b,a：所设计的滤波器分子和分母多项式系数向量。2.用窗函数法设计FIR滤波器原理：根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为w，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2采样N点，采用离散傅里叶反变换（IDFT）即可求出。用窗函数wd(n)将hd(n)截断，并进行加权处理，得到如果要求线性相位特性， 则h(n)还必须满足：根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。 验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。3.用FFT对信号作频分析的原理：用FFT对信号作频分析是学习数字信号处理的重要内容，经常需要进行分析的信号是模拟信号的时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是2/N，因此要求2/N小于等于D。可以根据此式选择FFT的变换区间N。误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时，离散谱的包络才能逼近连续谱，因此N要适当选择大一些。周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号的周期，可以尽量选择信号的观察时间长一些。对模拟信号进行频谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。三、课程设计过程与结果1.IIR（无限脉冲响应）模拟滤波器设计：(1)15阶椭圆低通滤波器设计：a. 设计指标：通带截止频率为fp=50Hz，通带最大衰减Rp=3dB，阻带最小衰减Rs=40dB，椭圆阶数N=15。b.设计程序：N=15;wp=50*2*pi;Rp=3;Rs=40;W=0:0.1:100*2*pi;B,A=ellip(N,Rp,Rs,wp,low,s);H=freqs(B,A,W);subplot(2,1,1);plot(W/2/pi,20*log10(abs(H);xlabel(频率/Hz);ylabel(幅值/dB);title(15阶椭圆模拟低通滤波器幅频特性);c.运行结果：（2） 切比雪夫I型带通滤波器设计：a.设计指标：通带上下边界频率为100Hz和200Hz，阻带上下边界频率为80Hz和220Hz，通带最大衰减rp=2dB，阻带最小衰减rs=50dB。b.设计程序：wp=100,200*2*pi;ws=80,220*2*pi;rp=2;rs=50;W=0:1:300*2*pi;n,wn=cheb1ord(wp,ws,rp,rs,s);B,A=cheby1(n,rp,wn,bandpass,s);H=freqs(B,A,W);subplot(2,1,1);plot(W/2/pi,20*log10(abs(H);xlabel(频率/Hz);ylabel(幅值/dB);title(切比雪夫I型模拟带通滤波器幅频特性);c.运行结果：2.IIR（无限脉冲响应）数字滤波器设计：（1） 设计巴特沃斯数字低通，并对录取的声音信号进行处理a. 设计指标：通带截止频率为fp=0.2Hz，阻带截止频率为fp=0.3Hz，通带最大衰减Rp=0.1dB，阻带最小衰减Rs=45dB，采样频率Fs=8000Hz。b. 设计程序：A,Fs,bits=wavread(E:/dp.WAV);N=length(A);t=0:1/Fs:(N-1)/Fs;figure(name,原始音频信号);subplot(211);plot(t,A);xlabel(t/s);ylabel(y/V);title(原始信号时域波形);h1=fft(A,N);w=0:Fs/N:(N/2-1)/N*Fs;subplot(212);plot(w,abs(h1(1:N/2);xlabel(t/s);ylabel(幅值);title(原始信号频谱图);wp=0.2*2*pi;ws=0.3*2*pi;Rp=0.1;As=45;ripple=10(-Rp/20);Attn=10(As/20);T=1/Fs;Omgp=wp*Fs;Omgs=ws*Fs;n,Omgc=buttord(Omgp,Omgs,Rp,As,s)z0,p0,k0=buttap(n);ba1=k0*real(poly(z0);aa1=real(poly(p0);ba,aa=lp2lp(ba1,aa1,Omgc);bd,ad=impinvar(ba,aa,Fs)Y=filter(bd,ad,A);figure(name,滤波器特性及处理后声音信号);h2=fft(Y,N);subplot(211);freqz(bd,ad,256,Fs);title(巴特沃斯数字低通滤波器幅频特性);subplot(212);plot(w,abs(h2(1:N/2);xlabel(t/s);ylabel(幅值);title(滤波后信号频谱图);c.运行结果：n = 18Omgc = 1.1308e+04bd = 0.0000 -0.0000 0.0000 0.0000 0.0004 0.0047 0.0212 0.0441 0.0466 0.0262 0.0079 0.0013 0.0001 0.0000 0.0000 0.0000 0.0000 -0.0000ad = Columns 1 through 18 1.0000 -3.4929 7.5957 -11.8036 14.2908 -13.9860 11.3137 -7.6586 4.3657 -2.0992 0.8495 -0.2876 0.0806 -0.0184 0.0033 -0.0005 0.0000 -0.0000 Column 190.0000(2) 设计切比雪夫I型带通滤波器，并对录取的声音信号进行处理a.设计指标：通带上下边界频率为1000Hz和2500Hz，阻带上下边界频率为800Hz和2800Hz，通带最大衰减rp=2dB，阻带最小衰减rs=50dB，采样频率Fs=8000Hz。b.设计程序：A,Fs,bits=wavread(E:/dp.WAV);N=length(A);t=0:1/Fs:(N-1)/Fs;figure(name,原始音频信号);subplot(211);plot(t,A);xlabel(t/s);ylabel(y/V);title(原始信号时域波形);h1=fft(A,N);w=0:Fs/N:(N/2-1)/N*Fs;subplot(212);plot(w,abs(h1(1:N/2);xlabel(t/s);ylabel(幅值);title(原始信号频谱图);wp=1000,2500*2*pi;ws=800,2800*2*pi;rp=2;rs=50;W=0:1:4000*2*pi;n,wn=cheb1ord(wp,ws,rp,rs,s);B,C=cheby1(n,rp,wn,bandpass,s);bd,ad=impinvar(B,C,Fs)Y=filter(bd,ad,A);figure(name,滤波器特性及处理后声音信号);h2=fft(Y,N);H=freqs(B,C,W);subplot(2,1,1);plot(W/2/pi,20*log10(abs(H);xlabel(频率/Hz);ylabel(幅值/dB);title(切比雪夫I型带通滤波器幅频特性);subplot(212);plot(w,abs(h2(1:N/2);xlabel(t/s);ylabel(幅值);title(滤波后信号频谱图);axis(0,4000,0,8)c.运行结果：bd = Columns 1 through 19 -0.0000 0.0000 0.0000 0.0000 -0.0004 0.0004 0.0021 -0.0063 0.0050 0.0056 -0.0140 0.0093 0.0012 -0.0049 0.0023 -0.0000 -0.0002 0.0000 0.0000 Columns 20 through 21 0.0000 0ad = Columns 1 through 19 1.0000 -4.6370 15.9377 -38.9633 81.1653 -141.2146 218.9300 -297.9495 368.2681 -407.3919 412.3852 -375.5223 312.8779 -233.1286 157.7425 -93.5296 49.4397 -21.7568 8.1911 Columns 20 through 21 -2.1796 0.44173.FIR（有限脉冲响应）数字滤波器设计:(1)利用矩形窗设计带阻滤波器，并对录取的声音信号进行处理a.设计指标：N=61，wc1=0.3，wc2=0.7，采样频率Fs=8000Hz。b.设计程序：A,Fs,bits=wavread(E:/dp.WAV);N=length(A);t=0:1/Fs:(N-1)/Fs;figure(name,原始音频信号);subplot(211);plot(t,A);xlabel(t/s);ylabel(y/V);title(原始信号时域波形);h1=fft(A,N);w=0:Fs/N:(N/2-1)/N*Fs;subplot(212);plot(w,abs(h1(1:N/2);xlabel(t/s);ylabel(幅值);title(原始信号频谱图);M=61;n=0:M-1;%输入设计指标wc1=0.3;wc2=0.7;f=0,wc1,wc1,wc2,wc2,1;%建立理想幅频特性频率向量m=1,1,0,0,1,1;%建立理想幅频特性幅度向量windows=boxcar(M);%使用矩形窗b=fir2(M-1,f,m,windows);%求FIR滤波器的系数Y=filter(b,1,A);figure(name,滤波器特性及处理后声音信号);h2=fft(Y,N);subplot(211);freqz(b,1,1024,Fs);title(boxcar窗FIR数字带阻滤波器频率特性);subplot(212);plot(w,abs(h2(1:N/2);xlabel(t/s);ylabel(幅值);title(滤波后信号频谱图);c.运行结果：(2)利用汉宁窗设计高通滤波器，并对录取的声音信号进行处理a.设计指标：N=41，wc=0.5，采样频率Fs=8000Hz。b.设计程序：A,Fs,bits=wavread(E:/dp.WAV);N=length(A);t=0:1/Fs:(N-1)/Fs;figure(name,原始音频信号);subplot(211);plot(t,A);xlabel(t/s);ylabel(y/V);title(原始信号时域波形);h1=fft(A,N);w=0:Fs/N:(N/2-1)/N*Fs;subplot(212);plot(w,abs(h1(1:N/2);xlabel(t/s);ylabel(幅值);title(原始信号频谱图);wc=0.5;M=41;n=0:M-1;f=0,wc,wc,1;m=0,0,1,1;b=fir2(M-1,f,m)Y=filter(b,1,A);figure(name,滤波器特性及处理后声音信号);h2=fft(Y,N);subplot(211);freqz(b,1,1024,Fs);title(汉明窗FIR数字高通滤波器频率特性);subplot(212);plot(w,abs(h2(1:N/2);xlabel(t/s);ylabel(幅值);title(滤波后信号频谱图);axis(0,4000,0,8)c.运行结果：b = Columns 1 through 19 0.0000 0.0011 -0.0000 -0.0020 -0.0000 0.0039 -0.0000 -0.0073 -0.0000 0.0125 -0.0000 -0.0206 -0.0000 0.0330 -0.0000 -0.0542 -0.0000 0.1002 -0.0000 Columns 20 through 38 -0.3163 0.5000 -0.3163 0.0000 0.1002 0.0000 -0.0542 0.0000 0.0330 0.0000 -0.0206 0.0000 0.0125 0.0000 -0.0073 0.0000 0.0039 0.0000 -0.0020 Columns 39 through 410.0000 0.0011 -0.00004.利用FFT进行频谱分析:a. 设计指标：利用FFT对信号s=10*cos(2*pi*2*t)+10*cos(2*pi* 2.05*t)+10*cos(2*pi*1.9*t)做频谱分析。b. 设计程序：fs=10;N=128;t=0:1/fs:N/fs;s=10*cos(2*pi*2*t)+10*cos(2*pi*2.05*t)+10*cos(2*pi*1.9*t);subplot(221);plot(t,s);xlabel(t/s);ylabel(幅值);title(混合频率信号);h1=fft(s(1:64),64);w1=0:fs/64:(31)/64*fs;subplot(222);plot(w1,abs(h1(1:32);xlabel(f/Hz);ylabel(|H(f)|);title(64点FFT);h2=fft(s(1:64),128);w2=0:fs/128:63/128*fs;subplot(223);plot(w2,abs(h2(1:64);xlabel(f/Hz);ylabel(|H(f)|);title(128点补零FFT);h3=fft(s(1:128),128);w3=0:fs/128:63/128*fs;subplot(224);plot(w3,abs(h3(1:64);xlabel(f/Hz);ylabel(|H(f)|);title(128点FFT);C.运行结果：4、 课设心得在数字信号与处理的学习中，已经使用过MATLAB，对其有了一定的了解和认识。这次课程设计给了我一个充分了解和学习MATLAB的机会，通过这