中国科学技术大学电子工程与信息科学系多媒体通信实验室(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