信号处理工具箱.ppt

第7章信号处理箱,7.1信号及其表示7.2信号的基本运算7.3线性时不变系统7.4傅里叶(Fourier)变换7.5IIR数字滤波器的设计方法7.6FIR数字滤波器设计7.7信号处理的图形用户界面工具,7.1信号及其表示,7.1.1工具箱中的信号产生函数,7.1.2离散时间信号的表示,在MATLAB中，离散时间信号x(n)的表示：需用一个向量x表示序列幅值，用另一个等长的定位时间变量n，才能完整地表示一个序列。,例7-10绘制离散时间信号的棒状图。其中x(-1)=-1,x(0)=1,x(1)=2,x(2)=1,x(3)=0,x(4)=-1。MATLAB源程序为：n=-3:5;%定位时间变量x=0,0,-1,1,2,1,-1,0,0;stem(n,x);grid;%绘制棒状图line(-3,5,0,0);%画x轴线xlabel(n);ylabel(xn)运行结果如图7.10所示。,图7.10离散时间信号图形,7.1.3几种常用离散时间信号的表示,1单位脉冲序列,直接实现：x=zeros(1,N);x(1,n0)=1;,2单位阶跃序列,直接实现：n=ns:nf;x=(n-n0)=0;,3实指数序列,直接实现：n=ns:nf;x=a.n;,4复指数序列,直接实现：n=ns:nf;x=exp(sigema+jw)*n);,5正(余)弦序列,直接实现：n=ns:nf;x=cos(w*n+sita);,7.2信号的基本运算,7.2.1信号的相加与相乘y(n)=x1(n)+x2(n)y(n)=x1(n)x2(n)MATLAB实现：y=x1+x2;y=x1.*x2,7.2.2序列移位与周期延拓运算,序列移位：y(n)=x(n-m)。MATLAB实现：y=x;ny=nx-m序列周期延拓：y(n)=x(n)M，MATLAB实现：ny=nxs:nxf；y=x(mod(ny,M)+1),7.2.3序列翻褶与序列累加运算,序列翻褶：y(n)=x(-n)。MATLAB可实现:y=fliplr(x),序列累加的数学描述为：,MATLAB实现：y=cumsum(x),7.2.4两序列的卷积运算,两序列卷积运算：,MATLAB实现：y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有限。,7.2.5两序列的相关运算,两序列相关运算：,。MATLAB实现：y=xcorr(x1,x2)。,7.2.6信号的能量和功率,1.信号能量,数字定义：,MATLAB实现:E=sum(x.*conj(x);或E=sum(abs(x).2);,数字定义：,2.信号功率,MATLAB实现:P=sum(x.*conj(x)/N;或E=sum(abs(x).2)/N;,7.3.1二次分式模型,连续系统：,离散系统：,7.3线性时不变系统,在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间模型、传递函数模型和零极点增益模型。其中sos表示二次分式，g为比例系数，sos为L6的矩阵，即,线性系统模型的变换函数,1.求连续LTI系统的零输入响应函数initial()格式：y,t,x=initial(a,b,c,d,x0)功能：计算出连续时间LTI系统由于初始状态x0所引起的零输入响应y。其中x为状态记录，t为仿真所用的采样时间向量。,2.求离散系统的零输入响应函数dinitial()格式：y,x,n=dinitial(a,b,c,d,x0)功能：计算离散时间LTI系统由初始状态x0所引起的零输入响应y和状态响应响应x，取样点数由函数自动选取。n为仿真所用的点数。,7.3.2线性时不变系统的响应,例7.16二阶系统,当初始状态,时，求系统的零输入响应。,MATLAB程序如下：a=-0.55,-0.78;0.78,0;b=1;0;c=5.966.45;d=0;x0=1;0;t0=0:0.1:20;initial(a,b,c,d,x0,t0);title(LTI系统的零输入响应曲线),运行结果如图7.15所示。,3滤波函数filter,MATLAB信号处理工具箱中提供了滤波函数filter。格式：y=filter(B,A,x)功能：对向量x中的数据进行滤波处理,例7.18设系统差分方程为,求该系统对信号,的响应。,MATLAB程序如下：B=1;A=1,-0.8;n=0:31;x=0.8.n;y=filter(B,A,x);subplot(2,1,1);stem(x)subplot(2,1,2);stem(y),程序运行结果如图7.17所示。,时间离散、连续频率序列傅里叶变换,离散时间、离散频率离散傅里叶级数,离散时间、离散频率离散傅里叶变换(DFT),正变换：,逆变换：,正变换：,逆变换：,正变换：,逆变换：,7.4傅里叶变换,1一维快速正傅里叶变换函数fft格式：X=fft(x,N)功能：采用FFT算法计算序列向量x的N点DFT变换，当N缺省时，fft函数自动按x的长度计算DFT。当N为2整数次幂时，fft按基-2算法计算，否则用混合算法。2一维快速逆傅里叶变换函数ifft格式：x=ifft(X,N)功能：采用FFT算法计算序列向量X的N点IDFT变换。,例7-19用快速傅里叶变换FFT计算下面两个序列的卷积。,并测试直接卷积和快速卷积的时间。,图7.18快速卷积框图,MATLAB程序(部分)：%线性卷积xn=sin(0.4*1:15);%对序列x(n)赋值,M=15hn=0.9.(1:20);%对序列h(n)赋值,N=20yn=conv(xn,hn);%直接调用函数conv计算卷积%园周卷积L=pow2(nextpow2(M+N-1);Xk=fft(xn,L);Hk=fft(hn,L);Yk=Xk.*Hk;yn=ifft(Yk,L);,图7.19x(n),h(n)及其线性卷积波形,7.5IIR数字滤波器的设计方法,1.数字滤波器的频率响应函数,幅度响应：,相位响应：,理想低通、高通、带通、带阻数字滤波器幅度特性,2.滤波器的技术指标幅度响应指标、相位响应指标,数字低通滤波器的幅度特性,通带要求：,阻带要求：,通带最大衰减：,阻带最小衰减：,7.5.1冲激响应不变法,2.MATLAB信号处理工箱中的专用函数impinvar()：格式：BZ,AZ=impinvar(B,A,Fs)功能：把具有B,A模拟滤波器传递函数模型转换成采样频率为Fs(Hz)的数字滤波器的传递函数模型BZ,AZ。采样频率Fs的默认值为Fs=1。,1.冲激响应不变法设计IIR数字滤波器的基本原理：,例7.20利用函数impinvar将一个低通模拟滤波器转变成数字滤波器，采样频率为10Hz。,MATLAB程序如下：b,a=butter(4,0.3,s);bz,az=impinvar(b,a,10);sys=tf(b,a);impulse(sys);绘制模拟滤波器的冲激响应holdonimpz(10*bz,az,10);绘制数字滤波器的冲激响应axis(020-0.040.12),程序运行结果如图7.20所示。,MATLAB信号处理工具箱中的专用双线性变换函数bilinear()格式：numd,dendbilinear(num,den,Fs)功能：把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。,7.5.2双线性变换法,双线性变换利用频率变换关系：,例7.21一个三阶模拟Butterworth低通滤波器的传递函数为,试用双线性变换法求出数字Butterworth低通滤波器的传递函数。,MATLAB程序如下：num=1;模拟滤波器系统函数的分子den=1,sqrt(3),sqrt(2),1;模拟滤波器系统函数的分母num1,den1=bilinear(num,den,1)求数字滤波器的传递函数程序运行结果为num1=0.05330.15990.15990.0533den1=1.0000-1.33820.9193-0.1546,7.5.3IIR数字滤波器的频率变换设计法,图7.21IIR数字滤波器MATLAB设计步骤流程图,1MATLAB的典型设计,利用在MATLAB设计IIR数字滤波器可分以下几步来实现(1)按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标；(2)根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频率Wc；(3)利用最小阶数N产生模拟低通滤波原型；(4)利用截止频率Wc把模拟低通滤波器原型转换成模拟低通、高通、带通或带阻滤波器；(5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。,例7.24设计一个在通带内的最大衰减为3dB，在阻带内的最小衰减为40dB的4阶低通模拟椭圆滤波器原型。MATLAB程序如下：n=4;rp=3;rs=40;z,p,k=ellipap(n,rp,rs);b,a=zp2tf(z,p,k);w=logspace(-1,1);freqs(b,a),运行结果如图7.24所示。,2MATLAB的直接设计,图7.25IIR数字滤波器MATLAB设计步骤流程图,例7-25试设计一个带阻IIR数字滤波器，其具体的要求是：通带的截止频率：wp1650Hz、wp2850Hz；阻带的截止频率：ws1700Hz、ws2800Hz；通带内的最大衰减为rp0.1dB；阻带内的最小衰减为rs50dB；采样频率为Fs2000Hz。MATLAB源程序设计如下：wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs=2000;wp=wp1,wp2/(Fs/2);ws=ws1,ws2/(Fs/2);%利用Nyquist频率频率归一化N,wc=ellipord(wp,ws,rp,rs,z);%求滤波器阶数num,den=ellip(N,rp,rs,wc,stop);%求滤波器传递函数H,W=freqz(num,den);%绘出频率响应曲线plot(W*Fs/(2*pi),abs(H);grid;xlabel(频率/Hz);ylabel(幅值)该程序运行后的幅频响应曲线如图7.26所示。,图7.26椭圆带阻滤波器的幅频响应,7.6FIR数字滤波器设计,格式：w=boxcar(M)功能：返回M点矩形窗序列。MATLAB信号处理工具箱中的窗函数法设计FIR数字滤波器的专用命令fir1()。格式：Bfir1(N,wc)功能：设计一个具有线性相位的N阶(N点)的低通FIR数字滤波器，返回的向量B为滤波器的系数(单位冲激响应序列)，其长度为N+1。,7.6.1窗函数设计法窗函数设计的基本原理：h(n)=w(n)hd(n)w(n)为窗函数，hd(n)理想数字滤波器的单位冲激响应。在MATLAB信号处理工具箱中为用户提供了Boxcar(矩形)、Bartlet(巴特利特)、Hanning(汉宁)等窗函数，这些窗函数的调用格式相同。,FIR数字滤波器的单位冲激响应h(n)满足偶(奇)对称h(n)=h(N-n-1)或h(n)=-h(N-n-1)FIR数字滤波器具有线性相位：,或,例7.27用矩形窗、三角窗、汉宁窗、哈明窗分别设计低通数字滤波器。信号的采样频率为1000Hz，数字滤波器的截止频率为100Hz，滤波器的阶数为80。,MATLAB程序如下：passrad=0.2*pi;w1=boxcar(81);w2=triang(81);w3=hanning(81);w4=hamming(81);n=1:1:81;hd=sin(passrad*(n-41)./(pi*(n-41);hd(41)=passrad/pi;h1=hd.*rot90(w1);h2=hd.*rot90(w2);h3=hd.*rot90(w3);h4=hd.*rot90(w4);MAG1,RAD=freqz(h1);,MAG2,RAD=freqz(h2);MAG3,RAD=freqz(h3);MAG4,RAD=freqz(h4);subplot(2,2,1);plot(RAD,20*log10(abs(MAG1);gridonsubplot(2,2,2);plot(RAD,20*log10(abs(MAG2);gridonsubplot(2,2,3);plot(RAD,20*log10(abs(MAG3);gridonsubplot(2,2,4);plot(RAD,20*log10(abs(MAG4);gridon,程序运行结果如图7.28所示。,7.6.2频率抽样法,1.频率抽样法的基本原理对所期望的滤波器的频率响应Hd(ejw)进行等间隔采样获得H(k)，利用h(n)=IDFTH(k)求得FIR的单位冲激响应。2.MATLAB信号处理工具箱中的频率抽样法专用函数命令fir2()格式：Bfir2(N,F,A)功能：设计一个N阶的FIR数字滤波器，其频率响应由向量F和A指定，滤波器的系数（单位冲激响应）返回在向量B中，长度为N+1。向量F和A分别指定滤波器的采样点的频率及其幅值，F中的频率必须在0.0到1.0之间，1.0对应于采样频率的一半。它们必须按递增的顺序从0.0开始到1.0为结束。,例7-33试用频率抽样法设计一个FIR低通滤波器，该滤波器的截止频率为0.5pi，频率抽样点数为33。MATLAB源程序为：N=32;F=0:1/32:1;%设置抽样点的频率，抽样频率必须含0和1。A=ones(1,16),zeros(1,N-15);%设置抽样点相应的幅值B=fir2(N,F,A);fr