数字信号处理实验

1一月二月三月产品名称数量金额利润产品名称数量金额利润产品名称数量金额利润合计合计合计四月五月六月产品名称数量金额利润产品名称数量金额利润产品名称数量金额利润合计合计合计实验七一、实验目的加深理解IIR数字滤波器的特性，掌握IIR数字滤波器的设计原理与设计方法，以及IIR数字滤波器的应用。二、实验原理N阶IIR数字滤波器的系统函数为:IIR数字滤波器的设计主要通过成熟的模拟滤波器设计方法来实现：将数字滤波器设计指标转换为模拟滤波器设计指标，设计出相应的模拟滤波器H(s)，再经过脉冲响应不变法或双线性变换法得到所需的IIR数字滤波器H(z)。IIR数字滤波器设计的重要环节是模拟原型低通滤波器的设计，主要包括Butterworth、Chebyshev和椭圆等滤波器MATLAB信号处理工具箱中提供了IIR滤波器设计的函数。IIR滤波器阶数选择buttord-巴特沃斯(Butterworth)滤波器阶数选择。cheb1ord-切比雪夫(Chebyshev)I型滤波器阶数选择。cheb2ord-切比雪夫(Chebyshev)II型滤波器阶数选择。ellipord-椭圆(Elliptic)滤波器阶数选择。IIR滤波器设计butter-巴特沃斯(Butterworth)滤波器设计cheby1-切比雪夫(Chebyshev)I型滤波器设计cheby2-切比雪夫(Chebyshev)II型滤波器设计ellip-椭圆(Elliptic)滤波器设计maxflat-通用的巴特沃斯(Butterworth)低通滤波器设计yulewalk-Yule-Walker滤波器设计(直接数字滤波器设计法)1.Butterworth滤波器设计Butterworth滤波器是通带、阻带都单调衰减的滤波器。(1)调用buttord函数确定巴特沃斯滤波器的阶数，格式为[N,Wc]=buttord(Wp,Ws,Ap,As)输入参数：Ap,As为通带最大衰减和阻带最小衰减，以dB为单位。Wp,Ws为归一化通带截频和阻带截频，0Wp,Ws1。输出参数：N为滤波器的阶数；Wc为截频，0Wc1。(2)调用butter函数设计出巴特沃斯滤波器，格式为[b,a]=butter(N,Wc,options)iiNijjMjzazbzH101)(NNNNMMMMzazazazbzbzbb)1(111)1(11101输入参数：N和Wc是buttord函数返回的参数，含义见上。Options=’low’,’high’,’bandpass’,’stop’,分别对应低通、高通、带通、带阻，默认情况下为低通或带通。输出参数：b和a为设计出的IIR数字滤波器H(s)的分子多项式和分母多项式的系数矩阵。2.ChebyshevII型滤波器设计ChebyshevII型滤波器为阻带纹波控制器：在阻带呈现纹波特性。[N,Wc]=cheb2ord(Wp,Ws,Ap,As)[b,a]=cheby2(N,As,Wc,options)3.椭圆滤波器设计椭圆滤波器在通阻带都呈现纹波特性。[N,Wc]=ellipord(Wp,Ws,Ap,As)[b,a]=ellip(N,Ap,As,Wc,options)三、实验内容1.信号,确定设计指标，实现各种IIR数字滤波器以实现以下信号处理。（1）设计IIR低通滤波器，滤除的成分。（2）设计IIR高通滤波器，滤除的成分。（3）设计IIR带通滤波器，滤除的成分。（4）设计IIR带阻滤波器，滤除的成分。要求利用butterord函数求解滤波器的阶数；利用butter函数设计各IIR数字滤波器；画出滤波器的幅度相应和相位响应；给出IIR数字滤波器的系统函数。(1)clear;fsam=100;t0=1/fsam;t=6;k=0:t0:t;N=512;x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0);f=(-N/2:(N/2-1))/N*2*pi;[N1,Wc]=buttord(1/4,1/2,3,60);[b,a]=butter(N1,Wc,'low');freqz(b,a);figure;axis([0,1,-120,0]);y=filter(b,a,x);Y=fftshift(fft(y,N));stem(f,abs(Y));)32cos()4cos(1][kkkx)32cos(k)4cos(1k)32cos(1k)4cos(k(2)clear;fsam=100;t0=1/fsam;t=6;k=0:t0:t;N=512;x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0);f=(-N/2:(N/2-1))/N*2*pi;[N1,Wc]=buttord(2/3,10/24,3,100);[b,a]=butter(N1,Wc,'high');freqz(b,a);figure;axis([0,1,-120,0]);y=filter(b,a,x);Y=fftshift(fft(y,N));stem(f,abs(Y));（3）clear;fsam=100;t0=1/fsam;t=6;k=0:t0:t;N=512;x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0);f=(-N/2:(N/2-1))/N*2*pi;[N1,Wc]=buttord([1/88/24],[1/322/3],3,60);[b,a]=butter(N1,Wc,'bandpass');freqz(b,a);figure;axis([0,1,-120,0]);y=filter(b,a,x);Y=fftshift(fft(y,N));stem(f,abs(Y));(4)clear;fsam=100;t0=1/fsam;t=6;k=0:t0:t;N=512;x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0);f=(-N/2:(N/2-1))/N*2*pi;[N1,Wc]=buttord([1/322/3],[0.210/27],3,60);[b,a]=butter(N1,Wc,'stop');freqz(b,a);figure;axis([0,1,-120,0]);y=filter(b,a,x);Y=fftshift(fft(y,N));stem(f,abs(Y));2.某带通滤波器的设计指标为（1）试分别利用巴特沃斯、切比雪夫Ⅰ型、切比雪夫Ⅱ型和椭圆模拟滤波器，通过脉冲响应不变法设计该带通数字滤波器，画出其频率特性，比较设计结果。（2）试分别利用巴特沃斯、切比雪夫Ⅰ型、切比雪夫Ⅱ型和椭圆模拟滤波器，通过双线性变换法设计该带通数字滤波器，画出其频率特性，比较设计结果。（3）分析比较以上设计结果，有何结论？(1).脉冲响应不变法clear;omegas=[0.20.72]*pi;omegap=[0.30.6]*pi;Ap=1;As=42;Fs=1;ws=omegas*Fs;wp=omegap*Fs;[N,Wc]=buttord(wp,ws,Ap,As,'s');[b,a]=butter(N,Wc,'bandpass','s');%[N,Wc]=cheb1ord(wp,ws,Ap,As,'s');%[b,a]=cheby1(N,Ap,Wc,'bandpass','s');%[N,Wc]=cheb2ord(wp,ws,Ap,As,'s');%[b,a]=cheby2(N,As,Wc,'bandpass','s');%%[N,Wc]=ellipord(wp,ws,Ap,As,'s');dBAdBAradradradradspspps42,1,)(72.0,)(6.0,)(3.0,)(2.02211%[b,a]=ellip(N,Ap,As,Wc,'bandpass','s');[numd,dend]=impinvar(b,a,Fs);w=linspace(0,pi,512);h=freqz(numd,dend,w);norm=max(abs(h));numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm));gridon右图从上到下分别是buttordCheby1cheby2ellipord型滤波器的频率特性(脉冲响应不变法)(2)双线性变换法clear;omegas=[0.20.72]*pi;omegap=[0.30.6]*pi;Ap=1;As=42;Fs=1;T=1/Fs;wp=2*tan(omegap/2)/T;ws=2*tan(omegas/2)/T;%[N,Wc]=buttord(wp,ws,Ap,As,'s');%[b,a]=butter(N,Wc,'bandpass','s');%[N,Wc]=cheb1ord(wp,ws,Ap,As,'s');%[b,a]=cheby1(N,Ap,Wc,'bandpass','s');%%[N,Wc]=cheb2ord(wp,ws,Ap,As,'s');%[b,a]=cheby2(N,As,Wc,'bandpass','s');[N,Wc]=ellipord(wp,ws,Ap,As,'s');[b,a]=ellip(N,Ap,As,Wc,'bandpass','s');[numd,dend]=bilinear(b,a,Fs);w=linspace(0,pi,512);h=freqz(numd,dend,w);norm=max(abs(h));numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm));gridontitle('ellipord')（1）利用DTMF的产生原理合成数字“4”，合成公式为：，画出其时域波形和频谱图。（2）根据图中的低通和高通滤波器，画出两个滤波器的幅频响应；并利用对信号进行滤波；利用fft函数求得输出的频响，绘出幅频响应曲线。（3）设计图中的第一组的四个带通滤波器，画出滤波器的幅频响应；并利用对图中相应的输入信号进行滤波，在一张图上用不同颜色画出这四条输出曲线。比较其幅值，BP770Hz滤波器的输出最大，因此确定该DTMF的行频为770Hz。(1)clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);stem(k,x);title('时域波形')figure;X=fftshift(fft(x));omega=((-(N-1)/2):((N-1)/2))/N*2*pi;stem(omega,abs(X));title('频谱图'))800012092sin()80007702sin(][4kkkx204:0k),,(xabfiltery),,(xabfiltery][4kx(2)clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);omega=((-(N-1)/2):((N-1)/2))/N*2*pi;%低通wp=2*941/8000;ws=2*1200/8000;Ap=3;As=60;[N,wc]=buttord(wp,ws,Ap,As);[b,a]=butter(N,wc);freqz(b,a);figure%滤波输出y=filter(b,a,x);Y=fftshift(fft(y));stem(omega,abs(Y))clear;k=0:204;N=205;x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k);omega=((-(N-1)/2):((N-1)/2))/N*2*pi;%高通wp=2*1200/8000;ws=2*1000/8000;Ap=3;As=60;[N,wc]=buttord(wp,ws,Ap,As);[b,a]=butter(N