MATLAB下的数字信号处理实现示例10号江城xin要点.docx

1、MATLAB下的数字信号处理实现示例10号江城xin要点MATLAB下的数字信号处理实现示例 附录一 信号、系统和系统响应 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; x=A*exp(-a*n*T).*sin(w0*n*T); %pi是MATLAB定义的，信号乘可采用“.*” close all %清除已经绘制的x(n)图形 subplot(3,1,1);stem(x); %绘制x(n)

2、的图形 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）改变参数为：T=1, =2.0734, a=0.4, A=1n=0:50; %定义序列的长度是50 A=1;

3、%设置信号有关的参数 a=0.4; T=1; %采样率 w0=2.0734; 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(理想采样信号序列); 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

4、=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);

5、stem(magX);title(单位冲击信号的幅度谱); angX=angle(X); %绘制x(n)的相位谱 subplot(3,1,3);stem(angX) ; title (单位冲击信号的相位谱) 4、特定冲击串：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(单位冲击信号序列); k=-25:25; X=x*(exp(-j*pi/12.5).(n*k); magX=abs(X)

6、; %绘制x(n)的幅度谱 subplot(3,1,2);stem(magX);title(单位冲击信号的幅度谱); angX=angle(X); %绘制x(n)的相位谱 subplot(3,1,3);stem(angX) ; title (单位冲击信号的相位谱) #在MATLAB中。提供了卷积函数conv，即y=conv(x,h)，调用十分方便。例如： 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);st

7、em(hb);title(系统hbn); m=1:50; %定义序列的长度是50 A=444.128; %设置信号有关的参数 a=50*sqrt(2.0)*pi; T=0.001; %采样率 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(输入信号xn); y=conv(x,hb); subplot(3,1,3);stem(y);title(输出信号yn); 6、卷积定律验证 k=-25:25; X=x*(exp(-j*pi/12.5)

8、.(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(

9、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(输出信号的幅度谱); angY=angle(Y); %绘制y(n)的相位谱 subplot(3,2,6);stem(angY) ; title (输出信号的相位谱) %以下将验证的结果显示 XHb=X.*Hb; figure,Subplot(2,1,1);stem(abs(XHb);title(x(n)的幅度谱与hb(n)幅度谱相乘); Sub

10、plot(2,1,2);stem(abs(Y);title(y(n)的幅度谱);axis(0,60,0,8000) %连续信号的傅里叶变化，在matlab中实现命令为fft( )，求连续信号的福谱图则用abs(fft(x)。n=0:15; %定义序列的长度是15 p=8;q=2; x=exp(-1*(n-p).2/q); 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); sub

11、plot(3,1,3); stem(abs(fft(x) 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) for i=0:3 x(i)=i+1;x(i+4)=8-(i+4); end for i=8:15 x(i)=0; end close all; subplot(2,1,1); stem(x); subplot(2,1,2); stem(abs(fft(x,16) % fft(

12、x,16)求关于x的16个点的傅立叶变换附录三 窗函数法设计FIR滤波器 一、在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）布拉克曼窗（Bla

13、ckman Window） 调用格式：w=blackman(n) ，根据长度n产生一个布拉克曼窗w。 （6）恺撒窗（Kaiser Window） 调用格式：w=kaiser(n,beta) ，根据长度n和影响窗函数旁瓣的参数产生一个恺撒窗w。 二、基于窗函数的FIR滤波器设计 利用MATLAB提供的函数firl 来实现 调用格式：firl (n,Wn,ftype,Window)，n为阶数、Wn是截止频率，取值范围为0 Wn 1.0（如果输入是形如W1 W2的矢量时，本函数将设计带通滤波器，其通带为W1W2）、ftype是滤波器的类型（低通-省略该参数、高通-ftype=high、带阻-ftyp

14、e=stop）、Window是窗函数。 例 设计一个长度为8的线性相位FIR滤波器。 Window=boxcar(8); %产生长度为8的矩形窗b=fir1(7,0.4,Window); freqz(b,1) % freqz（）计算滤波器的频率响应，此处freqz(b,1)等同于freqz(b),只是看作带分母且系数为1Window=blackman(8); %加布拉克曼窗b=fir1(7,0.4,Window); freqz(b,1) 例 设计线性相位带通滤波器，其长度N=15，上下边带截止频率分别为W1= 0.3，w2=0.5 Window=blackman(16); b=fir1(15,

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; tr_width=ws-wp; M=ceil(6.6*pi/tr_width)+1; N=0:1:M-1; wc=(ws+wp)/2; hd=ideal_lp(wc,M); w_ham=(boxcar(M); h=hd.*w_ham; db,mag,pha,grd,w=freqz_m(h,1); delta_w=2*pi/1000; Rp=-(min(db(1:1:wp/delta_w+1); As=-

16、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(hdn); subplot(2,2,2);stem(w_ham);title(汉明窗); axis(0 M-1 0 1.1);ylabel(wn); subplot(2,2,3);stem(h);title(实际冲击响应); axis(0 M-1 0.1 0.3);ylabel(hn); subplot(2,2,4);plot(w/pi,db); title(衰减幅度); ax

17、is(0 1 -100 10);ylabel(Decibles); 附录四 IIR滤波器的实现 MATLAB中滤波器的分析和实现 1、freqs函数：模拟滤波器的频率响应 例 系统传递函数为 的模拟滤波器，在MATLAB中可以用以下程序来实现： a=1 0.4 1; b=0.2 0.3 1; w=logspace(-1,1); %产生从10-1到101之间的50个等间距点，即50个频率点freqs(b,a,w) %根据输入的参数绘制幅度谱和相位谱 2、freqz函数：数字滤波器的频率响应 程序来实现： a=1 0.4 1; b=0.2 0.3 1; %根据输入的参数绘制幅度谱和相位谱，得到0到

18、之间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)，根据阶

19、数n和截止频率Wn计算ButterWorth滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。 相同参数条件下的模拟滤波器则调用格式为：b,a=butter(n,Wn,s) 例 采样频率为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); 4、Chebyshev模拟和数字滤波器 （1）cheb1ord函数：Chebyshev型滤波

20、器阶数计算。 调用格式：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=

21、cheb1ord(0.2,0.3,1,25); b,a=cheby1(n,1,Wn); freqz(b,a,512,1); （1）脉冲响应不变法设计数字ButterWorth滤波器 调用格式：bz,az=impinvar(b,a,Fs)，再给定模拟滤波器参数b，a和取样频率Fs的前提下，计算数字滤波器的参数。两者的冲激响应不变，即模拟滤波器的冲激响应按Fs取样后等同于数字滤波器的冲激响应。 （2）利用双线性变换法设计数字ButterWorth滤波器 调用格式：bz,az=bilinearb,a,Fs，根据给定的分子b、分母系数a和取样频率Fs，根据双线性变换将模拟滤波器变换成离散滤波器，具有分

22、子系数向量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)