第7章 MATLAB在数字信号处理中的应用,——本章将介绍数字信号处理的相关知识、傅里叶变换、IIR数字滤波器、FIR数字滤波器的基本理论和MATLAB实现。,(一)典型信号及其表示,1. 单位抽样信号(又称Kronecker函数) 在MATLAB中用zeros函数实现，但需要注意的是，MATLAB只能实现有限长的序列，选择产生N点的单位抽样序列： x=zeros(1, N); x(1)=1;,,2.单位阶跃序列 在MATLAB中，可采用ones函数产生N点单位阶跃序列： x=ones(1,N);,,3.正弦序列 其中，f是正弦波的频率， 是抽样周期， 为正弦波初相， 为正弦波幅度。在MATLAB中定义变量后，下列MATLAB语句给出正弦波序列： n=0:N-1; x=A*sin(2*pi*f*n*TS+pha); %TS为抽样周期，pha为正弦波初相,,,,,4. sinc函数 连续sinc函数定义为： MATLAB函数提供了sinc命令来产生N点采样序列： x=-2*pi:4*pi/N:2*pi; y=sinc(x);,,5. 随机序列 MATLAB提供了多种命令产生不同概率密度分布的随机变量序列，其中产生均匀分布与高斯分布随机序列命令： (1)rand：产生均值为0.5，幅度在[0,1]间均匀分布的伪随机数，调用方式为u=rand(N)或u=rand(M,N)。 (2)randn：产生均值为0，方差为1，服从高斯(正态)分布的随机序列，调用方法与rand相同。 除此之外，MATLAB中还包含其它随机序列命令，如chirp、gauspuls、triplus等，读者若用到，可再进一步查找。,(二)离散时间系统,离散时间系统的特点 离散时间系统，可以抽象为一种变换，或者一种映射，即把输入序列 变换为输出序列 ： 其中T表示一种变换。,,,,在数字信号处理中，通常研究的都是线性时不变(LTI)系统(亦可称作线性移不变LSI系统)。 1.线性：指该系统的输入、输出之间满足叠加原理。 2.时不变性：指对给定的输入，系统的输出和输入施加的时间无关。,离散时间系统的表示方法,离散时间系统可以用时域和频域两种方法描述。在时域，LTI系统的输入输出关系可以用差分方程或单位抽样响应 ；在频域还可以采用转移函数、二次分式和零极点增益等形式。 单位抽样响应 频域响应 零极点增益形式,1. 单位抽样响应,单位抽样响应是当输入为单位抽样信号时系统的输出响应，因此利用单位抽样信号来求单位抽样响应。在数字信号处理中，如果一个离散时间系统是用作对输入信号做滤波处理，离散时间系统又可称为数字滤波器。可以利用MATLAB中的filter函数和单位抽样信号来获得离散时间的单位抽样响应。,,例：利用 求出系统的单位抽样响应。,delta=[1 zeros(1,63)]; %单位抽样信号 b=[0.6 -0.3 0.5]; a=[1 -0.7 -0.2]; h=filter(b,a,delta); %单位抽样响应 h1=impz(b,a,64)； %impz函数可直接根据b、a给出单%位抽样响应。 subplot(2,1,1) %比较两个函数给出的单位抽样响应 stem(h); xlabel('n');ylabel('h(n)'); title('由filter函数求得的单位抽样响应') subplot(2,1,2) stem(h1); xlabel('n');ylabel('h(n)'); title('由impz函数求得的单位抽样响应'),,2.频率响应,对于离散系统，频域响应由MATLAB函数freqz得到，函数用法为： [h,f]=freqz (b, a, n, ’whole’, fs)：b, a同例11.1由离散系统决定；fs为抽样频率，可自行定义，若fs=1，频率轴给出归一化频率；n是频率轴的分点数，建议取2的整次幂；whole指定计算的频率范围是从0~fs，缺省时是0~fs/2。,,,例：分区间显示函数值,function y=f(x) if x0 y=0; elseif x1 f=x; elseif x2 f=2-x; else f=0; end,例：利用 求出系统的频域响应。,b=[0.6 -0.3 0.5]; a=[1 -0.7 -0.2]; [h,f]=freqz(b,a,256,1); hr=abs(h); hph=angle(h)*180/pi; subplot(2,1,1) plot(f,hr); xlabel('归一化频率');ylabel('幅度') title('幅频响应') subplot(2,1,2) plot(f,hph) xlabel('归一化频率');ylabel('相位/角度') title('相频响应'),,3.零极点增益形式,零极点增益形式——利用MATLAB中的求根函数由离散系统差分方程得到系统的零极点。 例：根据上例差分方程求系统零极点 b=[0.6 -0.3 0.5]; a=[1 -0.7 -0.2]; zp=roots(b) pp=roots(a),4.离散系统转换,同一个离散系统，可以用多种表示方式，如转移函数、零极点增益、状态空间等等，MATLAB工具箱中提供了实现不同表示方法之间进行转换的函数，如前一章所述。,(三)傅里叶分析,离散时间信号的Z变换和傅里叶变换 离散傅里叶变换 快速傅里叶变换,1.离散时间信号的Z变换和傅里叶变换,Z变换是对离散序列进行的一种数学变换，定义为： 离散系统的系统函数定义为,,,离散时间信号的傅里叶变换 对于非周期序列，可定义其DTFT为 除离散信号的傅里叶变换之外，还包括离散周期信号的傅里叶级数，不再赘述。重点是时域、频域均离散且有限长的离散傅里叶变换(DFT),,2.离散傅里叶变换,离散傅里叶变换由离散傅里叶级数得到，仅在时域和频域取一个周期，公式为： 其中,,,例：给出有限长序列的[0.5 1 1 0.5]的傅里叶变换，并对其频域值求傅里叶反变换，再与原序列比较。,N=4; %序列长度 n=0:N-1; k=n; xn=[0.5 1 1 0.5]; %产生序列 wn=exp(-j*pi*2/N); nk=n'*k; wnk=wn.^nk; xk=xn*wnk; %计算DFT WNKI=wn.^(-nk); XNI=xk*WNKI/N; %计算 IDFT,离散傅里叶变换的性质,线性 正交性* 圆周移位性质 时域圆周卷积 Parseval定理,,,例 将序列[1 2 3 4]循环右移两位。,N=4; %序列长度 n=0:N-1; M=2; xn=[1 2 3 4]; %产生序列 nm=mod((n-M),N) xm=xn(nm+1); subplot(2,1,1),stem(xn); xlabel(‘n’);ylabel(‘幅值’); title('原始序列'); subplot(2,1,2),stem(xm); xlabel(‘n’);ylabel(‘幅值’); title('移位后序列'),利用圆周移位可以实现圆周卷积,function[y]=circonv(x1,x2,N) % 实现 x1 和x2的圆周卷积 % 输出序列为y % 卷积长度为N % 方法:y(n)=sum(x1(m)*x2((n-m)mod N)) % 检查x1的长度 if (length(x1)N|length(x2)N) error('N must be=length(x1)') end x1=[x1 zeros(1,N-length(x1))]; x2=[x2 zeros(1,N-length(x2))]; for n=1:N for m=1:N y(n,m)=x1(m)*x2(mod( (n-m),N)+1) end end y=sum(y');,例：求序列[1 2 0.5 5 3]和序列(0≤n8)的10个点的圆周卷积。,x1=[1 2 0.5 5 3]; n=0:7; x2=0.6.^n; N=10; y=circonv(x1,x2,N); subplot(3,1,1),stem(x1); title('x1') subplot(3,1,2),stem(x2); title('x2'); subplot(3,1,3),stem(y); title('y'),,3.快速傅里叶变换,J.W.Cooley和J.W.Tukey巧妙的利用因子的周期性及对称性，提出了快速傅里叶变换(FFT)算法， MATLAB中提供了内部函数fft来得到序列的快速傅里叶变换，函数格式： y=fft(x,n)：可以计算序列x的n点FFT。当x的长度小于n时，fft函数在x的尾部补零，以构成n点数据；当x的长度大于n时，fft函数对序列x进行截尾。为了提高运算速度，n通常取2的幂次方。 y=ifft(x,n)：计算序列的逆傅里叶变换,例：模拟信号 ，求其幅度谱和相位谱。,fs=100; %抽样频率 N=128; %抽样点数 n=0:N-1; xn=2*sin(3*pi*n/fs)-sin(6*pi*n/fs); xk=fft(xn,N); f=(0:N-1)*fs/N; xkam=abs(xk); xkan=angle(xk); subplot(2,1,1) plot(f,xkam) xlabel('频率(Hz)'),ylabel('幅值') title('幅频响应') axis([0 50 0 150]) subplot(2,1,2) plot(f,xkan) xlabel('频率(Hz)'),ylabel('弧度') title('相频响应') axis([0 50 -2 2]),,(四)数字滤波器,数字滤波器通过对采样数据信号进行数学运算来达到频域滤波的目的，在数字信号处理中发挥着重要的作用。当滤波器的输入、输出都是离散信号，则该滤波器的冲激响应也必然是离散的，这样的滤波器称为数字滤波器。数字滤波器中的数学运算有两种实现方式，一种是频域方法，即利用FFT快速算法对输入信号进行离散傅里叶变换，分析其频谱，然后根据所希望的频率特性进行滤波，再利用IFFT快速算法对滤波后的频域信号进行反变换得到时域信号，这种方法比较灵活、简单，并且计算快速。第二种是时域方法，通过对离散数据作差分方程运算来达到滤波的目的。,1.模拟低通滤波器,MATLAB中给出了几种不同类型的模拟滤波器原型。 低通模拟Butterworth滤波器原型——buttap [Z,P,K]=buttap(n)； Z，P，K分别是n阶Butterworth滤波器的零点、极点和增益。 低通模拟ChebyshevI型滤波器原型——cheb1ap,[Z,P,K]= cheb1ap (n,rp) Z，P，K分别是n阶ChebyshevI型滤波器的零点、极点和增益。滤波器在通带内的最大衰减为rp。ChebyshevI型滤波器的主要特点是在阻带内达到最大平滑。 低通模拟ChebyshevⅡ型滤波器原型——cheb2ap [Z,P,K]= cheb2ap (n,rs) Z，P，K分别是n阶ChebyshevⅡ型滤波器的零点、极点和增益。滤波器在阻带内的最小衰减为rs。ChebyshevⅡ型滤波器的主要特点是在通带内达到最大平滑。,低通模拟椭圆滤波器原型——ellipap [Z,P,K]=ellipap (n,rp,rs) Z，P，K分别是n阶椭圆滤波器的零点、极点和增益。滤波器在通带内的最大衰减为rp，在阻带内的最小衰减为rs。 低通模拟Bessel滤波器——besselap [Z,P,K]= besselap (n) Z，P，K分别是n阶低通模拟Bessel滤波器的零点、极点和增益。,例：分别设计20阶ChebyshevI型低通模拟滤波器，通带内的最大衰减为0.3dB；20阶低通模拟ChebyshevⅡ型滤波器阻带内的最小衰减为45dB，并给出其频率特性图。,[z1,p1,k1]=cheb1ap(20,0.3); [num1,den1]=zp2tf(z1,p1,k1); %将零极点形式转换为系统函数形式 [z2,p2,k2]=cheb2ap(20,45); [num2,den2]=zp2tf(z2,p2,k2); figure(1) %绘图 freqs(num1,den1) figure(2) freqs(num2,den2),对于模拟高通、带通、带阻滤波器，其设计方法为先将要设计的滤波器的技术指标通过某种频率转换关系转换成低通滤波器的技术指标，并依据这些指标设计出低通滤波器的转移函数，然后再依据频率转换关系变成所设计的滤波器的转移函数。 MATLAB的信号处理工具箱提供了从低通滤波器向低通、高通、带通、带阻滤波器转换的函数。,(1)低通向低通的转变——lp2lp [numt,dent] = lp2lp(num,den,Wo) [At,Bt,Ct,Dt] = lp2lp(A,B,C,D,Wo) 这两条语句分别将表示成传递函数形式与状态方程形式的具有任意截止频率的模拟低通滤波器原型转换成截止频率为Wo的低通滤波器。 (2)低通向高通的转变——lp2hp [numt,dent] = lp2hp(num,den,Wo) [At,Bt,Ct,Dt] = lp2hp(A,B,C,D,Wo) 这两条语句分别将表示成传递函数形式与状态方程形式的具有任意截止频率的模拟低通滤波器原型转换成截止频率为Wo的高通滤波器。,(3)低通向带通的转变——lp2bp [numt,dent] = lp2bp(num,den,Wo,Bw) [At,Bt,Ct,Dt] = lp2bp(A,B,C,D,Wo,Bw) 这两条命令分别将表示成传递函数形式与状态方程形式的具有任意截止频率的模拟低通滤波器原型转换成中心频率为Wo、带宽为Bw的带通滤波器。 (4)低通向高通的转变——lp2bs [numt,dent] = lp2bs(num,den,Wo,Bw) [At,Bt,Ct,Dt] = lp2bs(A,B,C,D,Wo,Bw) 这两条命令分别将表示成传递函数形式与状态方程形式的具有任意截止频率的模拟低通滤波器原型转换成中心频率为Wo、带宽为Bw的带阻滤波器。,例：采用两种方法设计十阶模拟椭圆低通滤波器，通带的最大衰减为3dB，阻带内的最大衰减为40dB，截止频率为8πrad。,[z,p,k]=ellipap(10,3,40); %根据要求设计零极点形式低通模拟椭圆滤波器 [A,B,C,D]=zp2ss(z,p,k); %零极点形式向状态方程形式转换 [at,bt,ct,dt]=lp2lp(A,B,C,D,8*pi); %给定截止频率的低通滤波器 [num2,den2]=ss2tf(at,bt,ct,dt) figure freqs(num2,den2) %%第二种方式 [z2,p2,k2]=ellipap(10,3,40); [num,den]=zp2tf(z,p,k) %零极点形式向传递函数形式转换 [num1,den1]=lp2lp(num,den,8*pi) figure freqs(num1,den1),为了设计数字滤波器，需要将模拟滤波器离散化，两个离散化基本方法为冲激响应不变法和双线性变换法， MATLAB中也有两种相应的函数。 1.冲激响应不变法 冲激响应不变法的基本原理是从滤波器的冲激响应g(t) 出发，对其以周期T进行采样得到离散序列g(nT) 作为数字滤波器的冲激响应。相应的MATLAB函数为impinvar，具体格式为 [bz,az] = impinvar(b,a,fs) 将具有[b,a]模拟滤波器传递函数模型转换成采样频率为fs的数字滤波器的传递函数[bz,az]。如果在函数中没有指定采样频率fs，默认为1Hz。,2.双线性变换法 为了克服冲激响应不变法产生的频率混叠现象，模拟向数字的转变通常采用双线性变换法。相应的MATLAB函数为bilinear，具体格式为 (1)[zd,pd,kd] = bilinear(z,p,k,fs) (2)[numd,dend] = bilinear(num,den,fs) (3)[Ad,Bd,Cd,Dd] = bilinear(A,B,C,D,fs) 以上三个命令分别将零极点形式、传递函数形式和状态方程形式的模拟滤波器转换成相应形式的数字滤波器，其中fs是采样频率。,2.无限长单位脉冲响应滤波器 --IIR,MATLAB设计IIR数字滤波器可分为以下几步来实现： (1) 按一定规则将数字滤波器的技术指标转换成模拟低通滤波器的技术指标； (2) 根据转换后的技术指标确定滤波器的最小阶数和截止频率； (3) 利用最小阶数产生模拟低通滤波器原型； (4) 利用截止频率把低通滤波器原型转换成模拟低通、高通、带通或带阻滤波器； (5) 利用冲激响应不变法或双线性变换法把模拟滤波器转换成数字滤波器。,表： 滤波器阶数选择函数,其中，N为滤波器阶数，Wc为截止频率，wp、ws是通带和阻带的截止频率，Wc、wp、ws均为被采样频率归一化后的频率，范围是[0,1]。,例：设计数字信号处理系统，采样频率fs=2000Hz，希望设计成Butterworth型高通数字滤波器，通带中允许的最大衰减为0.4dB，阻带内的最小衰减为40dB，通带上限临界频率为400Hz，阻带下限频率为300Hz,wp=400*2*pi;ws=300*2*pi;rp=0.4;rs=40;fs=2000; % 转换成模拟滤波器 [N,Wc]=buttord(wp,ws,rp,rs,'s'); %计算阶数和截止频率 [Z,P,K]=buttap(N); [A,B,C,D]=zp2ss(Z,P,K); [AT,BT,CT,DT]=lp2hp(A,B,C,D,Wc); [num1,den1]=ss2tf(AT,BT,CT,DT); [num2,den2]=bilinear(num1,den1,fs); %双线性法得到数字滤波器 % [num2,den2]=impinvar(num1,den1,fs) %冲激响应不变法 [H,W]=freqz(num2,den2); F=W*fs/2/pi; semilogy(F,abs(H));grid xlabel('频率(Hz)');ylabel('幅值');,MATLAB还提供了几个直接设计IIR数字滤波器的函数，这些函数也要与数字滤波器的阶数选择配合使用。 如设计Butterworth数字滤波器的butter函数： (1)[b,a] = butter(n,Wn)：设计n阶截止频率为Wn的Butterworth低通数字滤波器的传递函数模型系数[b,a]，系数长度为n+1。Wn为归一化截止频率，其最大值为采样频率的一半。当Wn为两元素向量Wn=[w1,w2]时，函数返回2N阶的带通数字滤波器，通带为w1ww2。 (2)[b,a] = butter(n,Wn,’high’)：设计高通滤波器系数b，a。 (3)[b,a] = butter(n,Wn,’stop’)：设计带阻滤波器系数b，a，频率Wn=[w1,w2]，阻带为w1ww2。 (4)[A,B,C,D] = butter(n,Wn)和[z,p,k] = butter(n,Wn)：分别返回所设计数字滤波器的状态模型系数[A,B,C,D]和零极点增益[z,p,k]。 其它函数可参考书中表格,例：一个数字滤波器的抽样频率为fs=2000Hz，试设计一个为此系统使用的带通滤波器，希望采用椭圆型滤波器。要求(1)通带范围为300Hz到400Hz，在通带边缘频率处的衰减不大于0.5dB，(2)在200Hz以下和500Hz以上衰减不小于15dB。,fs=2000; ws=[300,400]/(fs/2);wp=[200,500]/(fs/2); rp=0.5;rs=15; %滤波器参数转换 [N,wc]=ellipord(wp,ws,rp,rs); %阶数与截止频率 [num,den]=ellip(N,rp,rs,wc); %直接设计椭圆滤波器 [H,W]=freqz(num,den); F=W*fs/2/pi; plot(F,20*log10(abs(H)));grid xlabel('频率(Hz)');ylabel('幅值(dB)'); axis([0 1000 -50 0]),3.有限长单位脉冲响应滤波器,有限长单位脉冲响应滤波器FIR的系统函数为： FIR系统只有零点，因此FIR系统不像IIR系统那样易取得较好的通带和阻带衰减特性。 但FIR系统具有自己突出的优点：稳定且具有线性相位。,,有限长单位脉冲响应FIR滤波器的设计方法建立在对理想滤波器特性作某种近似的基础上，主要包括窗函数法和频率采样法。 (1).窗函数法 窗函数法设计的原理是，用一个截短后的单位抽样响应h(n)来逼近理想的非因果的单位抽样响应hd(n)。即 为窗函数。MATLAB信号处理工具箱提供了Boxcar(矩形)、Bartlet(巴特利特)、Hanning(汉宁)等窗函数，,,矩形窗： w=boxcar(M)：返回M点矩形窗序列。其频率响应的逼近程度好坏完全取决于窗函数的频率特性。 除了提供窗函数命令外，MATLAB还提供了专用窗函数法设计FIR滤波器的命令fir1。 b= fir1(n,Wn,'ftype')：用来设计具有线性相位的n阶FIR数字滤波器，返回的向量b为滤波器的系数。截止频率Wn必须在0到1之间，数值1对应于fs/2。滤波器的归一化增益在Wn处为-6dB。当Wn=[W1 W2]时，fir1返回一个N阶的带通滤波器，通带为[W1 W2]。,其中，’ftype’是可选参数，可以选择以下参数中的一种‘high’、‘low’、‘stop’、‘dc-1’、‘dc-0’，分别表示高通、低通、带阻、第一个频带为通带、第一个频带为阻带的意思。若Wn=[W1,W2, W3,W4,…,Wn]，为多通带滤波器，其频带为0 WW1，W1 W W2，…，Wn W 1。 (2)b = fir1(n,Wn,window)：window表示用指定的window窗设计FIR滤波器。默认情况下，fir1使用hamming窗，还可以指定使用boxcar、bartlett、blackman、kaiser等窗函数。,例： 用矩形窗设计线性相位FIR低通滤波器，通带截止频率为0.2，滤波器单位脉冲响应的长度为21。,M=21;wc=0.2*pi; r=(M-1)/2; nr=-r:r; hdn=sin(wc*nr)/pi./nr; if rem(M,2)~=0,hdn(r+1)=wc/pi;end w=boxcar(M); %矩形窗 hn=hdn.*w'; %时域冲激响应 hw=fft(hn,512);w1=2*[0:511]/512; %频谱特性 subplot(2,1,1),stem(1:M,hn,'o'); xlabel('n'),ylabel('h(n)'); subplot(2,1,2),plot(w1,20*log10(abs(hw))); xlabel('w/pi');ylabel('幅值(dB)'); axis([0 2 -90 10]);,例： 设计一个多通带滤波器，归一化的通带是:[0,0.2]，[0.4,0.6]，[0.8,1]。注意高频端为通带，故滤波器的阶数应该为偶数，取N=30。,wn=[0.2,0.4,0.6,0.8]; %设置多通带 N=30; b=fir1(N,wn,'dc-1'); %令第一个频带为通带 freqz(b) figure(2) stem(b,'.'); xlabel('n'),ylabel('h(n)');,(2)频率采样法 频域采样法的基本原理是对理想的滤波器的频率响应在频域进行采样，以此来确定FIR滤波器的H(k)。在MATLAB信号处理箱中，为频率采样法设计FIR滤波器提供了专用函数命令fir2。该函数所得FIR数字滤波器系数为实数，具有线性相位，且满足对称性。,fir2命令的格式如下： 1)b= fir2(n,f,m)：返回n阶FIR数字滤波器的系数b。向量f和m分别指定滤波器采样点的频率和幅值，f为归一化频率。所期望的滤波器的频率响应可用plot(f,m)绘制。 2)b = fir2(n,f,m,window)：用指定的窗函数来设计FIR数字滤波器，窗函数包括Boxar、Hann、Bartlett、Blackman、Kaiser及Chebwin等。默认情况下，使用hamming窗。,例： 使用fir2命令设计上例中的滤波器。,f=0:0.01:1 m(1:21)=1;m(22:41)=0;m(42:61)=1;m(62:81)=0;m(82:101)=1; %设置多通带 plot(f,m,'k:');hold on N=30; b=fir2(N,f,m); %令第一个频带为通带 [h,f1]=freqz(b); f1=f1./pi; plot(f1,abs(h)); legend('理想滤波器','设计滤波器'); xlabel('归一化频率'),ylabel('幅值');,谢谢！,

展开阅读全文