Matlab与信号处理信号处理.ppt

Matlab与信号处理信号处理4、基于MATLAB的数字信号处理4.1 4.1 信号信号产生生4.2 4.2 离散系离散系统的的MATLABMATLAB实现4.3 4.3 变换的的MATLABMATLAB实现4.4 4.4 数字数字滤波器波器设计4.5 4.5 功率功率谱估估计24.1 信号产生l信号的表示信号的表示x(n)={x(n)}={…,x(-2),x(-1),x(0),x(1),x(2) …}Matlab中，信号都是用矩阵（向量）描述例如：x=[4 3 7 –9 1]l典型离散信号表示典型离散信号表示Ø单位抽样序列 x=zeros(1,N); x(1)=1; δ(n) x(k)=1; δ(n-k)Ø单位阶跃序列: x=ones(1,N);3Ø波形的波形的产生生首先产生时间：如t=0:0.01:2*pi;产生信号：如 y=sin(2*t+pi/4);square sawtooth sinc Ø噪声的噪声的产生生rand , randn4 周期波形产生函数名功能Sawtooth(t,width)产生锯齿波或三角波. Width确定最大值位置Square(t,duty)产生方波。

Duty指定正半周期的比例Sinc(t)产生sinc或sin(pi*t)/(pi*t)Diric(t,n)产生Dirichlet函数5sawtooth(t,0.7)square(t)sinc(t)diric(t,5)6l信号的常规操作信号加x=x1+x2信号必须等长对齐信号乘x=x1.*x2信号必须等长对齐改变比例x=a*x1移位x=x1(n-k)7折叠y=x(-n)取样和y=sum(x(n1:n2))取样积y=prod(x(n1:n2))卷积y=conv(x,h)[y,ny]=conv_m(x,nx,h,nh)相关r= xcorr(x,y)8l产生数据生数据Ø直接直接键盘输入入Ø利用函数利用函数产生，如生，如sin,cos,sawtooh等等Ø用用Matlab从从.Mat文件文件导入入Ø通通过文件文件读写写产生数据，如生数据，如fopen,fread等等Ø用一个用一个.M文件文件专门产生数据生数据94.2 离散系统的MATLAB实现离散系离散系统的表示方法：的表示方法： LSI系系统的的时域表示域表示 常系数线性差分方程 单位抽样响应h(n) LSI系系统的的频域表示域表示 频率响应 （DFT） 转移函数 （z变换） 零极点增益 二次分式 离散系离散系统的内部描述的内部描述 10l离散系离散系统的的输入入输出的差分方程：出的差分方程：l单位抽位抽样响响应 ?filter函数y=filter(b,a,x)impz函数y=impz(b,a)1112l频率响应 ？ [h,f]=freqz(b,a,n,fs)1314l 零极点增益？15l离散系统的表示方法转换函数 tf2zp sos2tf tf2ss sos2zp zp2tf sos2ss zp2sos ss2tf zp2ss ss2zp ss2sos16 4.3 变换的MATLAB实现1、DFT17计算 DFT并画图 1819202、FFT 可利用内部函数fft进行计算，速度快。

y=fft(x,n) x的n点fft x=ifft(y,n)213、z变换[R,P,K]=residuez(b,a);224、DCTly=dct(x,n)lx=idct(y,n)23rp245、Hilbert变换l 解析信号 l DFT方法求解l 函数y=Hilbert(x)25（（1）） 数字数字滤波器的基本概念波器的基本概念（（2）基于）基于MATLAB的的IIR滤波器波器设计（（3）） 基于基于MATLAB的的FIR滤波器波器设计 基于基于MATLAB的数字的数字滤波器波器设计26（1） 数字数字滤波器的基本概念波器的基本概念l数字数字滤波器波器 对数字信号进行滤波的线性时不变系统l数字数字滤波器的分波器的分类 IIR数字滤波器：设计经典选频滤波器，性能好， 阶数低 FIR数字滤波器：设计线性相位经典滤波器，适用范围广、 阶数高 27l数字数字滤波器的波器的设计步步骤1） 按照实际任务的要求，确定滤波器的技术指标。

2）用一个因果、稳定的离散线性时不变系统的系统函数 去逼近这一网络性能根据不同的要求可选择IIR或者 FIR系统函数去逼近3）利用有限精度算法实现系统函数28lIIR数字滤波器一般有两种设计方法： 1）模拟原型法 先设计一个合适的模拟滤波器，然后再变 换为满足要求的数字滤波器 2）直接法 直接在离散时域和离散频域进行数字滤波 器的设计 29lFIR一般有三种设计方法： 窗函数法 频率取样法 切比雪夫逼近法 30（2） 基于基于MATLAB的的IIR滤波器波器设计lButterworth数字数字滤波器波器设计l椭圆数字数字滤波器波器设计lChebyshev-I数字数字滤波器波器设计lChebyshev-II数字数字滤波器波器设计31lButterworth数字数字滤波器波器设计 MATLAB语句： [b,a]=butter(n,wn) (1) [b,a]=butter(n,wn,’ftype’) (2) [z,p,k]=butter(…) (3) [a,b,c,d]=butter(…) (4) 语句(1)~ (4) 都可以设计一个阶数为n，截止频率为wn的数字低通滤波器。

对于语句(2) ，若参数ftype取‘hi-gh’或‘stop’，则可以来设计相应的高通或带阻滤波器 返回值a和b为系统函数的分子和分母的系数32系统函数为： 返回值z,p,k分别是滤波器的零点、极点和增益返回值A、B、C、D构造滤波器的状态方程状态方程为： 其中u是输入信号，x是状态变量，y是输出信号33例例 ： 设计一个阶数为9阶、截止频率为300Hz的高通 Butterworth滤波器，采样频率为1000并画出它的频率特性程序： n=9; wn=300/500; [b,a]=buffer(n,wn,’high’); freqz(b,a,128,1000);34Buttordworth高通滤波器的频率响应35 n=236例例 设计一个10阶的带通Butterworth滤波器，它的通带 范围是100~200Hz,采样频率为1000并画出它的频率特性和冲击响应 程序：n=10; wn=[100 200]/500; [b,a]=butter(n,wn); freqz(b,a,128,1000); [y,t]=impz(b,a,101); stem(t,y);37 Buttordworth带通滤波器的频率响应38Buttordworth滤波器的时域冲击响应39l椭圆数字数字滤波器波器设计 MATLAB语法： [b,a]=ellip(n,Rp,Rs,wn) (5) [b,a]=ellip(n, Rp,Rs, wn,’ftype’) (6) [z,p,k]=ellip(…) (7) [a,b,c,d]=ellip(…) (8) 语句(5)~ (8)都可以设计一个阶数为n，截止频率为wn,通带波纹最大衰减为Rp、阻带波纹最小衰减为Rs的数字低通滤波器。

对于语句(6)，若参数ftype取‘high’或‘stop’，则可以来设计相应的高通或带阻滤波器 返回值a和b为系统函数的分子和分母的系数返回 40值z,p,k分别是滤波器的零点、极点和增益返回值A、B、C、D可以用来构造滤波器的状态方程例例 采样频率为1000Hz，请设计一个阶数为6阶、截止频 率为300Hz的低通ellip数字滤波器其中滤波器在通 带的波纹为3dB，阻带的波纹为50dB 程序： [b,a]=ellip(6,3,50,300/500); freqz(b,a,512,1000); title('n=6 低通椭圆滤波器');41椭圆低通滤波器的频率响应42例例 设计一个10阶的带通滤波器，通带范围为100~200Hz， 信号的采样频率为1000Hz，滤波器在通带的波纹为 dB，阻带的波纹为20dB设计出这个滤波器并画出它 的冲击响应。

程序： [b,a]=ellip(10,0.5,20,[100 200]/500); [y,t]=impz(b,a,101); freqz(b,a,128,1000); stem(t,y); title('n=10 带通椭圆滤波器');43椭圆带通滤波器的频率响应44椭圆滤波器的冲击响应45lChebyshev-I数字数字滤波器波器设计 MATLAB语法： [b,a]=cheby1(n,Rp,wn) (9) [b,a]= cheby1(n,Rp,wn,’ftype’) (10) [z,p,k]= cheby1(…) (11) [a,b,c,d]= cheby1(…) (12) 语句(9)~ (12)都可以设计一个阶数为n，截止频率为wn,通带波纹最大衰减为Rp的数字低通滤波器。

对于语句(10)，若参数ftype取‘high’或‘stop’，则可以来设计相应的高通或带阻滤波器 返回值a和b为系统函数的分子和分母的系数返回46值z,p,k分别是滤波器的零点、极点和增益返回值A、B、C、D可以用来构造滤波器的状态方程例例 对采样频率为1000Hz的采样信号，设计一个阶数为9 阶，截止频率为300Hz的低通cheby1数字滤波器，其 中滤波器在通带的波纹为 程序： [b,a]=cheby1(9,0.5,300/500); freqz(b,a,512,1000); 47 Chebyshev-I数字低通滤波器的频率响应48lChebyshev-II数字数字滤波器波器设计 MATLAB语法： [b,a]=cheby2(n,Rs,wn) (13) [b,a]= cheby2(n,Rs,wn,’ftype’) (14) [z,p,k]= cheby2(…) (15) [a,b,c,d]= cheby2(…) (16) 语句(13)~ (16)都可以设计一个阶数为n，截止频率为wn, 阻带纹波最小衰减为Rs的数字低通滤波器。

对于语句(13)式，若参数ftype取‘high’或‘stop’，则可以来设计相应的高通或带阻滤波器 返回值a和b为系统函数的分子和分母的系数返回49值z,p,k分别是滤波器的零点、极点和增益返回值A、B、C、D可以用来构造滤波器的状态方程例例 对于采样频率为1000Hz的信号，设计一个阶数为9阶，截止频率为300Hz的低通Chebyshev-II数字滤波器，其中滤波器在阻带的波纹为20dB 程序： [b,a]=cheby2(9,20,300/500); freqz(b,a,512,1000); 50 Chebyshev-II数字低通滤波器的频率响应 51（（3 ）基于）基于MATLAB的的FIR滤波器波器设计l窗函数法窗函数法设计FIR数字数字滤波器波器l频率取率取样法法l切比雪夫逼近法切比雪夫逼近法设计FIR滤波器波器52l 窗函数法窗函数法设计FIR数字数字滤波器波器Ø 时域中进行的 Ø 常见的窗函数 w=boxcar(n) w=triang(n) w=hanning(n) w=hamming(n) w=blackman(n) w=chebwin(n,R) w=bartlett(n) w=kaiser(n, ) 53例：用8种窗函数设计数字低通滤波器。

采样频率 1000Hz，截止频率200Hz，滤波器阶数81 程序：5455Ø利用利用fir1函数函数设计滤波器波器（（窗函数法）窗函数法） MATLAB语法： b=fir1(n, wn) b=fir1(n, wn, ’ftype’) b=fir1(n, wn, window) b=fir1(n, wn, ’ftype’, window) 56例例 : 设计一个阶数为48，通带范围为 的带通FIR线性相位滤波器，并分析它的频率特性 程序：b=fir1(48,[0.35,0.65]); freqz(b); 57 利用FIR1设计的带通滤波器的频率响应58Ø利用利用kaiserord函数求函数求凯赛窗函数的参数窗函数的参数 MATLAB语法： [n,wn,beta,ftype]=kaiserord(f,a,dev) [n,wn,beta,ftype]=kaiserord(f,a,dev,Fs）例例 利用凯赛窗函数设计一个低通FIR数字滤波器，通带 范围是0~1000Hz，阻带范围是1500~4000Hz，通带 的波纹最大为，阻带的波纹最大为。

采样频 率为8000Hz59 程序: fsamp=8000; fcuts=[1000 1500]; mags=[1 0]; devs=[0.05 0.01]; [n,wn,beta,ftype]=kaiserord(fcuts,mags, devs,fsamp); b=fir1(n,wn,ftype,kaiser(n+1,beta)); freqz(b); 60 利用凯赛窗设计的低通滤波器的频率响应61Ø利用利用fir2函数函数设计任意响任意响应FIR数字数字滤波器波器 MATLAB语法： b=fir2(n, f, m) b=fir2(n, f, m, window) b=fir2(n, f, m, npt) b=fir2(n, f, m, npt, window) b=fir2(n, f, m, npt, lap) 62例例 设计一个60阶的滤波器，要求设计的滤波器在0~ 的幅度响应为1，在 到 的幅度响应为 ，在 到 的幅度响应为 ，在 到 的 幅度响应为 ，在 到 的幅度响应为 。

并画 出该滤波器和理想滤波器的幅频响应曲线，进行比较 程序：f=[0,0.125,0.125,0.25,0.25,0.5,0.5,0.75,0.75,1]; m=[1,1,0.5,0.5,0.25,0.25,1/6,1/6,0.125,0.125]; b=fir2(60, f, m); [h,w]= freqz(b); plot(f,m,w/pi,abs(h));63 用fir2设计的任意响应的滤波器与理想滤波器比较64l 切比雪夫逼近法设计FIR滤波器Ø 采用最大误差最小准则方法设计滤波器Ø 优化问题Ø 雷米兹（REMEZ)算法 较成功的设计线性相位滤波器的方法65利用利用remez函数函数设计FIR数字数字滤波器波器 MATLAB语法： b=remez(n, f, m) b=remez(n, f, m, w) b=remez(n, f, m, ’ftype’) b=remez(n, f, m, w, ’ftype’) …… 66其它的MATLAB函数 ： firls remezord fircs fircls1 cremez firrcos 67l IIR与FIR数字滤波器的比较Ø在相同的技术指标下，IIR可以用比FIR较少的阶数来满足要求；ØFIR可得到严格的线性相位；ØIIR的设计可利用模拟滤波器的结果ØIIR主要设计规格化的滤波器，而FIR较灵活，可是任意响应的滤波器，适应性较广。

68n一维信号处理：信号滤波信号：s1 = sin(2*pi*t*5); s2=sin(2*pi*t*15); s3=sin(2*pi*t*30);S=s1+s2+s3滤波器：[b,a] = ellip(4,0.1,40,[10 20]*2/Fs);[H,w] = freqz(b,a,512);滤波后的结果sf = filter(b,a,s);695 功率谱估计 概念 ： 平稳随机信号的功率谱密度（PSD）是自相关序列的离散时间傅立叶变换70 功率谱估计的方法71经典功率谱估计方法l直接法（周期图法） 直接由傅立叶变换而得来的72例子：直接法计算功率谱利用FFT直接法：Fs=1000;NFFT=1024;n=0:1/Fs:1;x=sin(2*pi*40*n)+4*sin(2*pi*100*n)+randn(size(n)); %噪声序列X=fft(x,NFFT);Pxx=abs(X).^2/length(n); %求解PSDt=0:round(NFFT/2-1);k=t*Fs/NFFT;P=10*(log10(Pxx(t+1)));plot(k,P)xlabel('Frequency(Hz)');ylabel('PSD(dB/Hz)');7374利用periodogram函数：Fs=1000;NFFT=1024;n=0:1/Fs:1;x=sin(2*pi*40*n)+4*sin(2*pi*100*n)+randn(size(n)); %噪声序列window=boxcar(length(x));[Pxx,f]=periodogram(x,window,NFFT,Fs); %求PSDPlot(f,10*log10(Pxx))75直接法功率谱图76l间接法（自相关法或BT法） 其理论基础是维纳-辛钦定理。

其方法是先由随机信号N个观察值估计出自相关函数 ，然后再求 的傅立叶变换77例子：间接法求功率谱估计Fs=2000;NFFT=1024;n=0:1/Fs:1;x=sin(2*pi*100*n)+4*sin(2*pi*500*n)+randn(size(n));%噪声序列Cx=xcorr(x,'unbiased'); %计算序列的自相关函数Cxk=fft(Cx,NFFT);Pxx=abs(Cxk); %求解PSDt=0:round(NFFT/2-1);k=t*Fs/NFFT;P=10*(log10(Pxx(t+1)));plot(k,P)title('自相关法功率谱估计');xlabel('Frequency(Hz)');ylabel('PSD(dB/Hz)');78间接法功率谱图79。