中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) 附录C】 MATLAB下的数字信号处理实现示例 本部分内容是本讲义中数据信号处理实验部分实验项目在MatLab下实现代码。之所以 提供这些代码，是希望通过研究以下代码，能够更快、更好地掌握用MatLab进行数据信号 处理实验的方法：提高实验质量。希望同学们在阅读代码的时候，注意学习方法，在最短的 时间内熟悉MatLab,提高应用能力。示例中有些部分是实验项目中的内容实现，有些是一 些典型例题的实现。研究示例代码，倡导个性化编程是我们的目标，希望同学们能在在进行 实验项目的过程中提高MatLab的应用能力：在学习MatLab编程的同时加强对数字信号处 理有关实验项目的理解。 以下代码段均在MatLab5.3下调试通过，但是由于排版或其他一些原因，可能有部分代 码段不能得到正常结果。您可以在“htp:/202.38.75.33 dsp/matlab/”得到本讲义的修订内容， 同时可以在这个网址获取所有代码。 附录C1信号、系统和系统响应 1、理想采样信号序列 (1)首先产生信号x(n),0<=n<=50 n=0:50: %定义序列的长度是50 A=444.128: %设置信号有关的参数 a=50*sqrt(2.0)*pi; T=0.001; %采样率 w0=50*sqrt(2.0)*pi %w符号在MatLab中不能输入，用w代替 x=A*exp(-a*n*T).*sin(wo*n*T); %pi是MATLAB定义的π，信号乘可采用“*” close all %清除已经绘制的x)图形 subplot(3,1,1);stem(x); %绘制x)的图形 title(理想采样信号序列)： %设置结果图形的标题 (2)绘制信号x(n)的幅度谱和相位谱 k=-25:25; W=(pi/12.5)*k X=x*(exp(-j*pi/12.5).(n'*k): magX=abs(X); %会制x)的幅度谱 subplot(3,l,2);stem(magX),tite(理想采样信号序列的幅度谱')方 angX=angle(X); %6会制x)的相位谱 subplot(3,l,3),stem(angX;title(理想采样信号序列的相位谱) (3)改变参数为：A=1,a=0.4,2。=2.0734,T=1 n=0:50: %定义序列的长度是50 A=1;a=0.4;w0=2.0734;T=1; %设置信号有关的参数和采样率T x=A*exp(-a*n*T).*sin(wo*n*T); %pi是M4TLAB定义的π，信号乘可采用“*” 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.cdu.cn cxh@ustc.ed山cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） 附录 C MATLAB 下的数字信号处理实现示例 本部分内容是本讲义中数据信号处理实验部分实验项目在 MatLab 下实现代码。之所以 提供这些代码，是希望通过研究以下代码，能够更快、更好地掌握用 MatLab 进行数据信号 处理实验的方法；提高实验质量。希望同学们在阅读代码的时候，注意学习方法，在最短的 时间内熟悉 MatLab，提高应用能力。示例中有些部分是实验项目中的内容实现，有些是一 些典型例题的实现。研究示例代码，倡导个性化编程是我们的目标，希望同学们能在在进行 实验项目的过程中提高 MatLab 的应用能力；在学习 MatLab 编程的同时加强对数字信号处 理有关实验项目的理解。 以下代码段均在 MatLab5.3 下调试通过，但是由于排版或其他一些原因，可能有部分代 码段不能得到正常结果。您可以在“http://202.38.75.33/dsp/matlab/”得到本讲义的修订内容， 同时可以在这个网址获取所有代码。 附录 C1 信号、系统和系统响应 1、理想采样信号序列 （1）首先产生信号 x(n),0<=n<=50 n=0:50; %定义序列的长度是 50 A=444.128; %设置信号有关的参数 a=50*sqrt(2.0)*pi; T=0.001; %采样率 w0=50*sqrt(2.0)*pi; %ω符号在 MatLab 中不能输入，用 w 代替 x=A*exp(-a*n*T).*sin(w0*n*T); %pi 是 MATLAB 定义的π，信号乘可采用“.*” close all %清除已经绘制的 x(n)图形 subplot(3,1,1);stem(x); %绘制 x(n)的图形 title(‘理想采样信号序列’); %设置结果图形的标题 （2）绘制信号 x(n)的幅度谱和相位谱 k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5)).^(n’*k); magX=abs(X); %绘制 x(n)的幅度谱 subplot(3,1,2);stem(magX);title(‘理想采样信号序列的幅度谱’); angX=angle(X); %绘制 x(n)的相位谱 subplot(3,1,3);stem(angX) ; title (‘理想采样信号序列的相位谱’) （3）改变参数为： A = 1,α = 0.4,Ω0 = 2.0734,T = 1 n=0:50; %定义序列的长度是 50 A=1; a=0.4; w0=2.0734; T=1; %设置信号有关的参数和采样率 T x=A*exp(-a*n*T).*sin(w0*n*T); %pi 是 MATLAB 定义的π，信号乘可采用“.*” 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) close all %清除已经绘制的x)图形 subplot(3,1,1);stem(x); %绘制x)的图形 title(理想采样信号序列')方 k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5).(n'*k): magX-abs(X); %会制xm)的幅度谱 subplot((3,l,2);stem(magX);tite(理想采样信号序列的幅度谱)方 angX=angle(X); %绘制x)的相位谱 subplot(3,l,3),stem(angX);title(理想采样信号序列的相位谱) 2、单位脉冲序列 在MatLab中，这一函数可以用zeros函数实现： n=1:50: %定义序列的长度是50 x=zeros(1.50): %注意：MATLAB中数组下标从I开始 x(1)=1;close all; subplot(3,l,I):stem(x);title(单位冲击信号序列)： k=-25:25: X=x*(exp(-j*pi/12.5)).(n'*k); magX=abs(X); %会制x)的幅度谱 subplot(3,1,2):stem(magX;title(单位冲击信号的幅度谱)； angX=angle(X); %会制x)的相位谱 subplot(3,l,3),stem(angX;title(单位冲击信号的相位谱) 3、矩形序列 n=1:5;0x=sign(sign(10-n)+1); close all;subplot((3,l,l),stem(x):title(单位冲击信号序列)： k=-25:25,X=x*(exp(-j*pi/25).(n*k, magX=abs(X); %绘制x)的幅度谱 subplot(3,l,2);stem(magX:title(单位冲击信号的幅度谱)； angX=angle(X); %绘制x)的相位谱 ubplot(3,l,3),stem(angX);title(单位冲击信号的相位谱) 4、特定冲击串 x(n)=6(n)+2.56(n-1)+2.56(n-2)+6(n-3) n=1:50; %定义序列的长度是50 x=zeros(1.50): %注意：ATLAB中数组下标从1开始 x(1)=1;x(2=2.5x(3)=2.5,x(4)=1; close all,subplot((3,1,1)stem(x),title(单位冲击信号序列)； 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.edu.cn cxh@ustc.cdu.cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） close all %清除已经绘制的 x(n)图形 subplot(3,1,1);stem(x); %绘制 x(n)的图形 title(‘理想采样信号序列’); k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5)).^(n’*k); magX=abs(X); %绘制 x(n)的幅度谱 subplot(3,1,2);stem(magX);title(‘理想采样信号序列的幅度谱’); angX=angle(X); %绘制 x(n)的相位谱 subplot(3,1,3);stem(angX) ; title (‘理想采样信号序列的相位谱’) 2、单位脉冲序列 在 MatLab 中，这一函数可以用 zeros 函数实现： n=1:50; %定义序列的长度是 50 x=zeros(1,50); %注意：MATLAB 中数组下标从 1 开始 x(1)=1;close all; subplot(3,1,1);stem(x);title(‘单位冲击信号序列’); k=-25:25; X=x*(exp(-j*pi/12.5)).^(n’*k); magX=abs(X); %绘制 x(n)的幅度谱 subplot(3,1,2);stem(magX);title(‘单位冲击信号的幅度谱’); angX=angle(X); %绘制 x(n)的相位谱 subplot(3,1,3);stem(angX) ; title (‘单位冲击信号的相位谱’) 3、矩形序列 n=1:5; 0x=sign(sign(10-n)+1); close all; subplot(3,1,1); stem(x);title(‘单位冲击信号序列’); k=-25:25; X=x*(exp(-j*pi/25)).^(n’*k); magX=abs(X); %绘制 x(n)的幅度谱 subplot(3,1,2);stem(magX);title(‘单位冲击信号的幅度谱’); angX=angle(X); %绘制 x(n)的相位谱 subplot(3,1,3);stem(angX) ; title (‘单位冲击信号的相位谱’) 4、特定冲击串 x(n) = δ (n) + 2.5δ (n −1) + 2.5δ (n − 2) + δ (n − 3) n=1:50; %定义序列的长度是 50 x=zeros(1,50); %注意：MATLAB 中数组下标从 1 开始 x(1)=1; x(2)=2.5; x(3)=2.5; x(4)=1; close all; subplot(3,1,1);stem(x);title(‘单位冲击信号序列’); 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) k=-25:25,X=x*(exp(-j*pi/12.5).n'*k方 magX-abs(X); %会制x)的幅度谱 subplot(3,l,2):stem(magX),title(单位冲击信号的幅度谱)； angX=angle(X); %绘制x)的相位谱 subplot(3,l,3),stem(angX;title(单位冲击信号的相位谱) 5、卷积计算 y(n)=x(n)*h(n)= 芝mMn-m 在MATLAB中。提供了卷积函数conv,即y=conv(x,h),调用十分方便。例如： 系统：h,(nm)=6(m)+2.56(n-1)+2.56(n-2)+6(n-3) 信号：x(t)=Ae sin(2onT),0≤n<50 n=1:50 %定义序列的长度是50 hb=zeros(1,50); %注意：MATLAB中数组下标从1开始 hb(1)）=1;hb(2)=2.5,hb(3)=2.5,hb(4=1; close all;subplot(3,1,1);stem(hb);title(hb[n]'); m=1:50,T=0.001; %定义序列的长度是和采样率 A=444.128,a=50*sqt(2.0)*pi %设置信号有关的参数 w0=50*sqrt(2.0)*pi x=A*exp(-a*m*T).*sin(w0*m*T); %pi是MATLAB定义的π，信号乘可采用《*” subplot((3,1,2),stem(x);t(输入信号x[n]')方 y=conv(x,hb); subplot(3,l,3)stem(y)tite(输出信号y[n')方 6、卷积定律验证 k=-25:25;X=x*(exp(-j*pi/12.5).(n*k); magX=abs(X); %绘制x()的幅度谱 subplot(3,2,1)stem(magX,tite(输入信号的幅度谱)方 angX-angle(X); %绘制x)的相位谱 subplot(3,2,2),stem(angX);title(输入信号的相位谱) Hb=hb*(exp(-j*pi/12.5)).(n'*k); magHb=abs(Hb); %绘制hb)的幅度谱 subplot(3,2,3);stem(magHb),title(系统响应的幅度谱)方： angHb-angle(Hb); %绘制hb)的相位谱 subplot((3,2,4);stem(angHb);title(系统响应的相位谱) n=1:99,k=1:99: Y=y*(exp(-j*pi/12.5)).^(n'*k); mag Y=abs(Y); %绘制y)的幅度谱 subplot(3,2,5);stem(magY):tie(‘输出信号的幅度谱)方 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.cdu.cn cxh@ustc.cdu.cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） k=-25:25; X=x*(exp(-j*pi/12.5)).^(n’*k); magX=abs(X); %绘制 x(n)的幅度谱 subplot(3,1,2);stem(magX);title(‘单位冲击信号的幅度谱’); angX=angle(X); %绘制 x(n)的相位谱 subplot(3,1,3);stem(angX) ; title (‘单位冲击信号的相位谱’) 5、卷积计算 ∑ +∞ =−∞ = ∗ = − m y(n) x(n) h(n) x(m)h(n m) 在 MATLAB 中。提供了卷积函数 conv，即 y=conv(x,h)，调用十分方便。例如： 系统： h (n) = (n) + 2.5 (n −1) + 2.5 (n − 2) + (n − 3) b δ δ δ δ 信号： ( ) sin( ),0 50 = Ω0 ≤ < − x t Ae nT n nT a α n=1:50; %定义序列的长度是 50 hb=zeros(1,50); %注意：MATLAB 中数组下标从 1 开始 hb(1)=1; hb(2)=2.5; hb(3)=2.5; hb(4)=1; close all; subplot(3,1,1);stem(hb);title(‘系统 hb[n]’); m=1:50; T=0.001; %定义序列的长度是和采样率 A=444.128; a=50*sqrt(2.0)*pi; %设置信号有关的参数 w0=50*sqrt(2.0)*pi; x=A*exp(-a*m*T).*sin(w0*m*T); %pi 是 MATLAB 定义的π，信号乘可采用“.*” subplot(3,1,2);stem(x);title(‘输入信号 x[n]’); y=conv(x,hb); subplot(3,1,3);stem(y);title(‘输出信号 y[n]’); 6、卷积定律验证 k=-25:25; X=x*(exp(-j*pi/12.5)).^(n’*k); magX=abs(X); %绘制 x(n)的幅度谱 subplot(3,2,1);stem(magX);title(‘输入信号的幅度谱’); angX=angle(X); %绘制 x(n)的相位谱 subplot(3,2,2);stem(angX) ; title (‘输入信号的相位谱’) Hb=hb*(exp(-j*pi/12.5)).^(n’*k); magHb=abs(Hb); %绘制 hb(n)的幅度谱 subplot(3,2,3);stem(magHb);title(‘系统响应的幅度谱’); angHb=angle(Hb); %绘制 hb(n)的相位谱 subplot(3,2,4);stem(angHb) ; title (‘系统响应的相位谱’) n=1:99; k=1:99; Y=y*(exp(-j*pi/12.5)).^(n’*k); magY=abs(Y); %绘制 y(n)的幅度谱 subplot(3,2,5);stem(magY);title(‘输出信号的幅度谱’); 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) angY=angle(Y); %绘制y的相位谱 subplot(3,2,6);stem(angY);title(输出信号的相位谱) %以下将验证的结果显示 XHb=X.*Hb; Subplot(2,l,l)stem(abs(XHb):title(x(n)的幅度谱与hb(n)幅度谱相乘)， Subplot(2,l,2),stem(abs(Y);tite(y(n)的幅度谱)，axis(0,60,0,8000]) 附录C2用FFT进行信号的频谱分析 1、高斯序列 (n-p)2 xa(n）= e ,0≤n≤15 0,else n=0:15; %定义序列的长度是15 p=8,q=2,X=exp(-1*(n-p).2/gi %利用历函数实现富氏变换 close all;subplot(3,1,1);stem(abs(fft(x))) p=8,q=4,X=exp(-1*(n-p).2/q): %改变信号参数，重新计算 subplot(3,1,2);stem(abs(fft(x))) p=8;q-8,X=exp(-1*(n-p).2/q: subplot(3,1,3);stem(abs(fft(x))) 2、衰减正弦序列 e-msin2gfi,0≤n≤15 x,(n)= 0,else n=0:15 %定义序列的长度是15 a=-0.l;f=0.0625,x=exp(-a*n).*sin(2*pi*f*n; close all;subplot(2,1,1);stem(x); subplot(2,1,2);stem(abs(fft(x))) 3、三角波序列 n+1,0≤n≤3 x.(n)=8-n,4≤n≤7 0,else for i=1:4 %设置信号前4个点的数值 x①)=i; %注意：MATLAB中数组下标从1开始 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.cdu.cn cxh@ustc.ed山cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） angY=angle(Y); %绘制 y(n)的相位谱 subplot(3,2,6);stem(angY) ; title (‘输出信号的相位谱’) %以下将验证的结果显示 XHb=X.*Hb; Subplot(2,1,1);stem(abs(XHb));title(‘x(n)的幅度谱与 hb(n)幅度谱相乘’); Subplot(2,1,2);stem(abs(Y);title(‘y(n)的幅度谱’); axis([0,60,0,8000]) 附录 C2 用 FFT 进行信号的频谱分析 1、高斯序列 = ≤ ≤ − − else x n e n q n p a 0, ( ) ,0 15 2 ( ) n=0:15; %定义序列的长度是 15 p=8; q=2; x=exp(-1*(n-p).^2/q); %利用 fft 函数实现富氏变换 close all; subplot(3,1,1); stem(abs(fft(x))) p=8; q=4; x=exp(-1*(n-p).^2/q); %改变信号参数，重新计算 subplot(3,1,2); stem(abs(fft(x))) p=8; q=8; x=exp(-1*(n-p).^2/q); subplot(3,1,3); stem(abs(fft(x))) 2、衰减正弦序列 ≤ ≤ = − else e fn n x n n b 0, sin 2 ,0 15 ( ) π α n=0:15; %定义序列的长度是 15 a=0.1; f=0.0625; x=exp(-a*n).*sin(2*pi*f*n); close all; subplot(2,1,1); stem(x); subplot(2,1,2); stem(abs(fft(x))) 3、三角波序列 − ≤ ≤ + ≤ ≤ = else n n n n xc n 0, 8 ,4 7 1,0 3 ( ) for i=1:4 %设置信号前 4 个点的数值 x(i)=i; %注意：MATLAB 中数组下标从 1 开始 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) end for i=5:8 %设置信号后4个点的数值 x(1=9-i: end close all;subplot(2,1,1);stem(x); %会制信号图形 subplot(2,1,2);stem(abs(fft(x,16))) %绘制信号的频谱 4、反三角序列 [4-n,0≤n≤3 xa(n)=n-3,4≤n≤7 0,else for i=1:4 %设置信号前4个点的数值 x①=5-i %注意：MATLAB中数组下标从I开始 end for i=5:8 %设置信号后4个点的数值 x0=i-4; end close all;subplot(2,1,1);stem(x); %绘制信号图形 subplot(2,1,2);stem(abs(fft(x,16))) %绘制信号的频谱 附录C3窗函数法设计FIR滤波器 1、在MATLAB中产生窗函数 (I)矩形窗(Rectangle Window) 调用格式：w=boxcar(n),根据长度n产生一个矩形窗w。 (2)三角窗(Triangular Window) 调用格式：w=triang(n),根据长度n产生一个三角窗w。 (3)汉宁窗(Hanning Window) 调用格式：w=hanning(n),根据长度n产生一个汉宁窗w。 (4)海明窗(Hamming Window) 调用格式：w=hamming(n)，根据长度n产生一个海明窗w。 (5)布拉克曼窗(Blackman Window) 调用格式：w=blackman(n)，根据长度n产生一个布拉克曼窗w。 (6)恺撒窗(Kaiser Window) 调用格式：w=kaiser(n,beta),根据长度n和影响窗函数旁瓣的B参数产生一个恺撒窗w。 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.cdu.cn cxh@ustc.ed山cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） end for i=5:8 %设置信号后 4 个点的数值 x(i)=9-i; end close all; subplot(2,1,1); stem(x); %绘制信号图形 subplot(2,1,2); stem(abs(fft(x,16))) %绘制信号的频谱 4、反三角序列 − ≤ ≤ − ≤ ≤ = else n n n n xd n 0, 3,4 7 4 ,0 3 ( ) for i=1:4 %设置信号前 4 个点的数值 x(i)=5-i; %注意：MATLAB 中数组下标从 1 开始 end for i=5:8 %设置信号后 4 个点的数值 x(i)=i-4; end close all; subplot(2,1,1); stem(x); %绘制信号图形 subplot(2,1,2); stem(abs(fft(x,16))) %绘制信号的频谱 附录 C3 窗函数法设计 FIR 滤波器 1、在 MATLAB 中产生窗函数 （1）矩形窗（Rectangle Window） 调用格式：w=boxcar(n)，根据长度 n 产生一个矩形窗 w。 （2）三角窗（Triangular Window） 调用格式：w=triang(n) ，根据长度 n 产生一个三角窗 w。 （3）汉宁窗（Hanning Window） 调用格式：w=hanning(n) ，根据长度 n 产生一个汉宁窗 w。 （4）海明窗（Hamming Window） 调用格式：w=hamming(n) ，根据长度 n 产生一个海明窗 w。 （5）布拉克曼窗（Blackman Window） 调用格式：w=blackman(n) ，根据长度 n 产生一个布拉克曼窗 w。 （6）恺撒窗（Kaiser Window） 调用格式：w=kaiser(n,beta) ，根据长度 n 和影响窗函数旁瓣的β参数产生一个恺撒窗 w。 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) 2、基于窗函数的FIR滤波器设计 利用MATLAB提供的函数fil来实现。 调用格式：firl(n,Wn,'fype',Window),n为阶数、Wn是截止频率（如果输入是形如WI W2]的矢量时，本函数将设计带通滤波器，其通带为W1<o<W2)、fype是滤波器的类型（低 通-省略该参数、高通-ftype=-high、带阻-fype=stop)、Window是窗函数。 [例设计一个长度为8的线性相位FR滤波器。 1,0≤o≤0.4π 其理想幅频特性满足|H.(e)仁 0,else Window=boxcar(8); %根据长度为8的矩形窗indow b=fir1(7,0.4,Window);freqz(b,1) Window=blackman(8)方 %根据长度为8的布拉克曼窗Window b=fir1(7,0.4,Window);freqz(b,1) [例]设计线性相位带通滤波器，设计指标：长度N=15,上下边带截止频率分别为W1= 0.3,w2=0.5π。 Window=blackman(16); b=fir1(15,[0.3 0.5],Window);freqz(b,1) 设计指标为：o0.2亚Rp=0.25dBωa=0.3TAs=50dB的低通数字FIR滤波器。 wp-0.2*pi;ws-0.3*pi;wc-(ws+wp)/2;tr_width-ws-wp; M=ceil(6.6*pi/tr_width)+1;N=[0:1:M-1]; alpha=(M-1)/2; n=0:1:(M-1)1: %n=0,M-1) m=n-alpha eps; hd=sin(wc*m)./(pi*m); w ham=(boxcar(M))': h=hd.*w ham; [H,w]=freqz(h,[1],1000,'whole); H=(H(1:501): w=(w(1:501)片 mag=abs(H); db=20*log10((mag+eps)/max(mag)); pha=angle(H); grd=grpdelay(h,[1],w); delta w=2*pi/1000: Rp=-(min(db(1:1:wp/delta_w+1))); As=-round(max(db(ws/delta_w+1:1:501))); Close all; subplot(2,2,1);stem(hd),title(理想冲击响应) axis([0 M-1 -0.1 0.3]);ylabel('hd[n]); subplot(2,2,2):stem(wham):title(汉明窗)： 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.edu.cn cxh@ustc.cdu.cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） 2、基于窗函数的 FIR 滤波器设计 利用 MATLAB 提供的函数 firl 来实现。 调用格式：firl (n,Wn,’ftype’,Window)，n 为阶数、Wn 是截止频率（如果输入是形如[W1 W2]的矢量时，本函数将设计带通滤波器，其通带为 W1<ω<W2）、ftype 是滤波器的类型（低 通-省略该参数、高通-ftype=high、带阻-ftype=stop）、Window 是窗函数。 [例] 设计一个长度为 8 的线性相位 FIR 滤波器。 其理想幅频特性满足| ≤ ≤ = else H e j d 0, 1,0 0.4 ( ) | ω ω π Window=boxcar(8); %根据长度为 8 的矩形窗 Window b=fir1(7,0.4,Window); freqz(b,1) Window=blackman(8); %根据长度为 8 的布拉克曼窗 Window b=fir1(7,0.4,Window); freqz(b,1) [例] 设计线性相位带通滤波器，设计指标：长度 N=15，上下边带截止频率分别为 W1= 0.3π，w2=0.5π。 Window=blackman(16); b=fir1(15,[0.3 0.5],Window); freqz(b,1) 设计指标为：ωp=0.2π Rp=0.25dB ωa=0.3π As=50dB 的低通数字 FIR 滤波器。 wp=0.2*pi; ws=0.3*pi; wc=(ws+wp)/2; tr_width=ws-wp; M=ceil(6.6*pi/tr_width)+1; N=[0:1:M-1]; alpha=(M-1)/2; n=[0:1:(M-1)]; %n=[0,(M-1)]; m=n-alpha + eps; hd=sin(wc*m)./(pi*m); w_ham=(boxcar(M))'; h=hd.*w_ham; [H,w]=freqz(h,[1],1000,'whole'); H=(H(1:501))'; w=(w(1:501))'; mag=abs(H); db=20*log10((mag+eps)/max(mag)); pha=angle(H); grd=grpdelay(h,[1],w); delta_w=2*pi/1000; Rp=-(min(db(1:1:wp/delta_w+1))); As=-round(max(db(ws/delta_w+1:1:501))); Close all; subplot(2,2,1);stem(hd);title('理想冲击响应') axis([0 M-1 -0.1 0.3]);ylabel('hd[n]'); subplot(2,2,2);stem(w_ham);title('汉明窗'); 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) axis(0 M-1 0 1.1);ylabel(w n); subplot(2,2,3):stem(h):tite(实际冲击响应)； axis([0 M-1 -0.1 0.3]);ylabel('h[n]); subplot(2,2,4)plot(wpi,db),title('衰减幅度)； axis([0 1-100 101);ylabel('Decibles'); 附录C4R滤波器的实现 1、freqs函数：模拟滤波器的频率响应 [例系统传递函数为H)=02s+03s+1 的模拟滤波器，在MATLAB中可以用以下 s2+0.4s+1 程序来实现： a=10.41;b=[0.20.31 %设置分子分母的系数 w=logspace(-l,l)月 %产生从10-1到10'之间地0个等间距点，即50个颜率点 freqs(b,a,w) %根据输入的参数绘制幅度谱和相位谱 2、freqz函数：数字波器的频率响应 [例]系统传递函数为H()= 0.2+0.3z+z2 1+0.4z+22 的模拟滤波器，在MATLAB中可以用以下 程序来实现： a=[10.41]b=[0.20.31; %根据输入的参数绘制幅度谱和相位谱，得到0到π之间128个点处的顿率响应 freqz(b,a,128) 3、ButterWorth模拟和数字滤波器 (1)butterd函数：Butter Worth滤波器阶数的选择。 调用格式：[n,Wn]=butterd(Wp,Ws,Rp,Rs),在给定滤波器性能的情况下（通带临界频率 Wp、阻带临界频率Ws、通带内最大衰减Rp和阻带内最小衰减Rs),计算Butter Worth滤波 器的阶数n和截止频率Wn。相同参数条件下的模拟滤波器则调用格式为： [n,Wn]=butterd(Wp,Ws,Rp,Rs,'s') (2)butter函数：ButterWorth滤波器设计。 调用格式：[b,a=butter(n,Wn),根据阶数n和截止频率Wn计算ButterWorth滤波器分 子分母系数(b为分子系数的矢量形式，a为分母系数的矢量形式)。相同参数条件下的模拟 滤波器则调用格式为：[b,a=butter(n,Wn,'s) [例]采样频率为1Hz,通带临界频率=0.2Hz,通带内衰减小于1dB(ap=1):阻带 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.cdu.cn cxh@ustc.cdu.cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） axis([0 M-1 0 1.1]);ylabel('w[n]'); subplot(2,2,3);stem(h);title('实际冲击响应'); axis([0 M-1 -0.1 0.3]);ylabel('h[n]'); subplot(2,2,4);plot(w/pi,db); title('衰减幅度'); axis([0 1 -100 10]);ylabel('Decibles'); 附录 C4 IIR 滤波器的实现 1、freqs 函数：模拟滤波器的频率响应 [例] 系统传递函数为 0.4 1 0.2 0.3 1 ( ) 2 2 + + + + = s s s s H s 的模拟滤波器，在 MATLAB 中可以用以下 程序来实现： a=[1 0.4 1]; b=[0.2 0.3 1]; %设置分子分母的系数 w=logspace(-1,1); %产生从 10−1 到101 之间地 0 个等间距点，即 50 个频率点 freqs(b,a,w) %根据输入的参数绘制幅度谱和相位谱 2、freqz 函数：数字滤波器的频率响应 [例] 系统传递函数为 2 2 1 0.4 0.2 0.3 ( ) − − + + + + = z z z z H z 的模拟滤波器，在 MATLAB 中可以用以下 程序来实现： a=[1 0.4 1];b=[0.2 0.3 1]; %根据输入的参数绘制幅度谱和相位谱，得到 0 到π之间 128 个点处的频率响应 freqz(b,a,128) 3、ButterWorth 模拟和数字滤波器 （1）butterd 函数：ButterWorth 滤波器阶数的选择。 调用格式：[n,Wn]=butterd(Wp,Ws,Rp,Rs)，在给定滤波器性能的情况下（通带临界频率 Wp、阻带临界频率 Ws、通带内最大衰减 Rp 和阻带内最小衰减 Rs），计算 ButterWorth 滤波 器的阶数 n 和截止频 率 Wn 。 相 同参数条 件下的模拟 滤波器则 调用格式为 ： [n,Wn]=butterd(Wp,Ws,Rp,Rs,’s’) （2）butter 函数：ButterWorth 滤波器设计。 调用格式：[b,a]=butter(n,Wn)，根据阶数 n 和截止频率 Wn 计算 ButterWorth 滤波器分 子分母系数（b 为分子系数的矢量形式，a 为分母系数的矢量形式）。相同参数条件下的模拟 滤波器则调用格式为：[b,a]=butter(n,Wn,’s’) [例] 采样频率为 1Hz，通带临界频率 fp =0.2Hz，通带内衰减小于 1dB（αp=1）；阻带 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) 临界频率f=0.3Hz,阻带内衰减大于25dB(as=25)。设计一个数字滤波器满足以上参数。 [n,Wn]=buttord(0.2,0.3,1,25); [b,a]=butter(n,Wn);fregz(b,a,512,1); 4、Chebyshev模拟和数字滤波器 (I)cheblord函数：Chebyshev I型Ⅱ滤波器阶数计算。 调用格式：[n,Wn]=cheblord(Wp,Ws,Rp,Rs),在给定滤波器性能的情况下（通带临界频 率Wp、阻带临界频率Ws、通带内波纹Rp和阻带内衰减Rs),选择Chebyshev I型滤波器 的最小阶n和截止频率Wn。 (2)chebyl函数：Chebyshev I型滤波器设计。 调用格式：b,a]=butter((n,Rp,Wn),根据阶数n、通带内波纹Rp和截止频率Wn计算 Butter Worth滤波器分子分母系数(b为分子系数的矢量形式，a为分母系数的矢量形式)。 注：ChebyshevⅡ型滤波器所用函数和I型类似，分别是cheb2ord、cheby2。 [例]实现上例中的滤波器 [n,Wn]=cheblord(0.2,0.3,1,25);[b,a]=cheby1(n,1,Wn); freqz(b.a.512,1): 5、滤波器设计 (I)脉冲响应不变法设计数字ButterWorth滤波器 调用格式：[bz,az=impinvar(b,a,Fs),再给定模拟滤波器参数b,a和取样频率Fs的前提 下，计算数字滤波器的参数。两者的冲激响应不变，即模拟滤波器的冲激响应按Fs取样后 等同于数字滤波器的冲激响应。 (2)利用双线性变换法设计数字Butter Worth滤波器 调用格式：[bz,az=bilinear[b,a,Fs],根据给定的分子b、分母系数a和取样频率Fs,根 据双线性变换将模拟滤波器变换成离散滤波器，具有分子系数向量bz和分母系数向量z。 模拟域的butter函数说明与数字域的函数说明相同b,a]=butter(n,Wn,'s)可以得到模拟域的 Butterworth滤波器 [例]采样频率为1Hz,通带临界频率=0.2Hz,通带内衰减小于1dB(ap=1):阻带 临界频率f=0.3Hz,阻带内衰减大于25dB(as=25)。设计一个数字滤波器满足以上参数。 %直接设计数字滤波器 [n,Wn]=buttord(0.2,0.3,1,25)月 [b,a]=butter(n,Wn);freqz(b,a,512,1); %脉冲响应不变法设计数字滤波器 [n,Wn]=buttord(0.2,0.3,1,25,'s) [b,a]=butter(n,Wn,'s');freqs(b,a) bz,az=impinvar(b,a,1);freqz(bz,az,512,1) %双线性变换法设计ButterWorth数字滤波器 [n,Wn=buttord(0.2,0.3,1,25,'s'): [b,a]=butter(n,Wn,'s');freqs(b,a) [bz,az]=bilinear(b,a,1);freqz(bz,az,512,1) 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.edu.cn cxh@ustc.cdu.cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） 临界频率 fs=0.3Hz，阻带内衰减大于 25dB（αs=25）。设计一个数字滤波器满足以上参数。 [n,Wn]=buttord(0.2,0.3,1,25); [b,a]=butter(n,Wn); freqz(b,a,512,1); 4、Chebyshev 模拟和数字滤波器 （1）cheb1ord 函数：ChebyshevⅠ型Ⅱ滤波器阶数计算。 调用格式：[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs)，在给定滤波器性能的情况下（通带临界频 率 Wp、阻带临界频率 Ws、通带内波纹 Rp 和阻带内衰减 Rs），选择 ChebyshevⅠ型滤波器 的最小阶 n 和截止频率 Wn。 （2）cheby1 函数：ChebyshevⅠ型滤波器设计。 调用格式：[b,a]=butter(n,Rp,Wn)，根据阶数 n、通带内波纹 Rp 和截止频率 Wn 计算 ButterWorth 滤波器分子分母系数（b 为分子系数的矢量形式，a 为分母系数的矢量形式）。 注：ChebyshevⅡ型滤波器所用函数和Ⅰ型类似，分别是 cheb2ord、cheby2。 [例] 实现上例中的滤波器 [n,Wn]=cheb1ord(0.2,0.3,1,25); [b,a]=cheby1(n,1,Wn); freqz(b,a,512,1); 5、滤波器设计 （1）脉冲响应不变法设计数字 ButterWorth 滤波器 调用格式：[bz,az]=impinvar(b,a,Fs)，再给定模拟滤波器参数 b，a 和取样频率 Fs 的前提 下，计算数字滤波器的参数。两者的冲激响应不变，即模拟滤波器的冲激响应按 Fs 取样后 等同于数字滤波器的冲激响应。 （2）利用双线性变换法设计数字 ButterWorth 滤波器 调用格式：[bz,az]=bilinear[b,a,Fs]，根据给定的分子 b、分母系数 a 和取样频率 Fs，根 据双线性变换将模拟滤波器变换成离散滤波器，具有分子系数向量 bz 和分母系数向量 az。 模拟域的 butter 函数说明与数字域的函数说明相同[b,a]=butter(n,Wn,’s’)可以得到模拟域的 Butterworth 滤波器 [例] 采样频率为 1Hz，通带临界频率 fp =0.2Hz，通带内衰减小于 1dB（αp=1）；阻带 临界频率 fs=0.3Hz，阻带内衰减大于 25dB（αs=25）。设计一个数字滤波器满足以上参数。 %直接设计数字滤波器 [n,Wn]=buttord(0.2,0.3,1,25); [b,a]=butter(n,Wn); freqz(b,a,512,1); %脉冲响应不变法设计数字滤波器 [n,Wn]=buttord(0.2,0.3,1,25,’s’); [b,a]=butter(n,Wn,’s’); freqs(b,a) [bz,az]=impinvar(b,a,1); freqz(bz,az,512,1) %双线性变换法设计 ButterWorth 数字滤波器 [n,Wn]=buttord(0.2,0.3,1,25,’s’); [b,a]=butter(n,Wn,’s’); freqs(b,a) [bz,az]=bilinear(b,a,1); freqz(bz,az,512,1) 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
中国科学技术大学电子工程与信息科学系多媒体通信实验室(Copyright2000) 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系network@ustc.edu.cn cxh@ustc.cdu.cn
中国科学技术大学电子工程与信息科学系 多媒体通信实验室 （Copyright 2000） 如果您在阅读过程中发现疏漏和错误，请您尽快和编者取得联系 network@ustc.edu.cn cxh@ustc.edu.cn
