matlab信号处理学习总结.doc

matlab信号处理学习总结 作者： 日期：62 个人收集整理 勿做商业用途常用函数 1 图形化信号处理工具，fdatool（滤波器设计），fvtool（图形化滤波器参数查看）sptool（信号处理)，fvtool(b，a），wintool窗函数设计。或者使用工具箱 filter design设计.当使用离散的福利叶变换方法分析频域中的信号时,傅里叶变换时可能引起漏谱,因此需要采用平滑窗, 2数字滤波器和采样频率的关系。如果一个数字滤波器的采样率为 FS，那么这个滤波器的分析带宽为Fs/2。也就是说这个滤波器只可以分析0，Fs/2的信号。举个例字: 有两个信号,S1频率为20KHz，S2频率为40KHz，要通过数字方法滤除S2。 你的滤波器的采样率至少要为Fs=80HKz，否则就分析不到 S2了，更不可能将它滤掉了！（当然根据采样定理，你的采样率 F0也必须大于80HK，,Fs和 F0之间没关系不大,可以任取，只要满足上述关系就行.） 3 两组数据的相关性分析 r=corrcoef（x，y） 4 expm 求矩阵的整体的 exp 4 离散快速傅里叶 fft信号处理中，傅里叶变换的典型用途是将信号分解成幅值分量和频率分量）.Ft为连续傅里叶变换。反傅里叶 ifft 5 ztrans（),Z变换是把离散的数字信号从时域转为频率 6 laplace（）拉普拉斯变换是把连续的的信号从时域转为频域 7 sound（x)会在音响里产生 x所对应的声音 8 norm求范数，det行列式,rank求秩 9 模拟频率，数字频率，模拟角频率关系 模拟频率f：每秒经历多少个周期，单位Hz，即1/s; 模拟角频率 是指每秒经历多少弧度，单位rad/s； 数字频率w：每个采样点间隔之间的弧度,单位rad。 =2pif; w = T 10 RMS求法 Rms = sqrt（sum(P.2）)或者norm（x）/sqrt(length（x） var方差的开方是std标准差,RMS应该是norm(x）/sqrt(length(x)）吧。 求矩阵的RMS：std（A（：） 11 ftshift 作用:将零频点移到频谱的中间 12 filtfilt零相位滤波， 采用两次滤波消除系统的非线性相位， y = filtfilt(b，a,x)；注意x的长度必须是滤波器阶数的3倍以上，滤波器的阶数由max（length(b)1,length（a）-1)确定. 13 h,t=impz（b，a,n，fs)，计算滤波器的冲激响应 h为n点冲击响应向量 h，x=freqz（b,a,n,fs）计算频响,有fs时，x为频率f,无fs，x为w角频率，常用于查看滤波器的频率特性 14 zplane(z，p) 画图零极点分布图 15 beta=unwarp（alpha) 相位会在穿越+-180发生回绕,可将回绕的 16 stepz 求数字滤波器的阶跃响应 h,t = stepz（b,a，n，fs) fvtool（b1，a1,b2，a2，。bn,an） fvtool(Hd1,Hd2,.。.) h = fvtool（。） 15 IIR数字滤波器设计方法 1 先根据已知带同参数求出最佳滤波器阶数和截止频率 n,Wn = buttord(Wp，Ws,Rp，Rs); n，Wn = buttord（Wp,Ws,Rp,Rs，s） b,a=butter（n,Wn，ftype，s） 其中Wp为，01之间。Ws为阻带角频率，01之间。Rp为通带波纹，或者通带衰减，Rs为阻带衰减. 若果给出的是模拟频率fp1通带截止频率,fp2阻带截止频率,则Wp=fp12/fs, Ws=fp2*2/fs.如果给出的是实际数字频率比如0.3*pi，如果给出的是 y=filter(b,a，x)；或者采用零相位滤波y=filtfilt（b，a，x) 15 传统 FIR滤波器 Ftype为滤波器类型，比如高通，低通，window为窗函数类型. Window窗函数。 例子1 设计一个通带滤波器，带宽为 0。35-0。65 b = fir1(48,0。35 0。65）； freqz(b,1,512) 16窗函数长度:窗函数的长度应等于 FIR滤波器系数个数,即滤波器阶数n+1。 17 加窗函数的 FIR滤波器长度的确定 17.1 buttord函数求出最佳滤波器阶数和截止频率，然后用 fir1函数调用，窗函数长度为滤波器最佳阶数 n+1 17.2 用窗函数方法设计 FIR滤波器，由滤波器的过渡带的宽度和选择的窗函数决定这里举一个选用海明窗函数设计低通滤波器的例子。 低通滤波器的设计要求是：采样频率为100Hz，通带截至频率为 3 Hz,阻带截止频率为 5 Hz,通带内最大衰减不高于0。5 dB，阻带最小衰减不小于50 dB.使用海明窗函数。确定 N的步骤有: 1，从上表可查得海明窗的精确过渡带宽为 6.6pi/N；（在有些书中用近似过渡带来计算，这当然没有错，但阶数增大了,相应也增加计算量。） 2，本低通滤波器的过渡带是： DeltaW=Ws-Wp=（5-3)pi/50=。04pi3，N=6。6pi/DeltaW=6。6pi /0。04pi=165 所以滤波器 的阶数至少是165。在该帖子中是用理想低通滤波器的方法来计算的,这里用fir1函数来计算，相应的程序有 fs=100； % 采样频率 wp = 3*pi/50； ws = 5*pi/50； deltaw= ws wp; 过渡带宽 的计算 N = ceil（6.6*pi/ deltaw） + 1; 按海明窗计算所需的滤波器阶数 N0 wdham = (hamming(N+1）; % 海明窗计算 Wn=(3+5）/100; 计算截止频率 b=fir1(N，Wn，wdham）; H,w=freqz(b，1); db=20*log10(abs（H）; % 画频响曲线 plot(wfs/(2*pi),db）;title（ 幅度响应（单位： dB)）；grid axis(0 50 -100 10）； xlabel（频率（单位：Hz）); ylabel(分贝) set(gca,XTickMode，manual,XTick,0,3,5，50） set(gca,YTickMode，manual，YTick，-50,0) 17 数字滤波器函数 Butter,cheyshev切比雪夫 b，a=cheby1(n，rp,wn，options), b，a=besself(n，wn,options） b，a=ellip（(n,rp,rs,wn，options） n为阶数，wn为截止频率rad/s，rs为阻带起伏.wn在0-1之间，且1对应于采样频率的一半。 b,a=butter(n,Wn,options）, z,p,k = butter(n,Wn，ftype，s） z，p，k = butter(n，Wn，ftype) A,B,C,D = butter(n,Wn，ftype，s） ftype对应 high 是高通滤波器的归一化截止频率 low 低通滤波器的归一化截止频率 stop for an order 2*n bandstop digital filter if Wn is a two-element vector, Wn = w1 w2. The stopband is w1 Ws; 2低通滤波器:Wp和 Ws为一元矢量且WpWs； 3带通滤波器：Wp和 Ws为二元矢量且WpWs，如Wp=0。1，0.8,Ws=0。2,0.7。 2 求出滤波器系数 Butter函数可设计低通、高通、带通和带阻的数字和模拟 IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度.在期望通带平滑的情况下,可使用 butter函数。 butter函数的用法为： b,a=butter（n，Wn，ftype）。其中 n代表滤波器阶数，Wn代表滤波器的截止频率二、契比雪夫 I型 IIR滤波器的设计 在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在 MATLAB下可使用cheby1函数设计出契比雪夫 I型 IIR滤波器。 cheby1函数可设计低通、高通、带通和带阻契比雪夫 I型滤 IIR波器，其通带内为等波纹，阻带内为单调。契比雪夫 I型的下降斜度比 II型大,但其代价是通带内波纹较大。 cheby1函数的用法为: b,a=cheby1(n，Rp,Wn,/ftype/） 在使用 cheby1函数设计 IIR滤波器之前，可使用 cheblord函数求出滤波器阶数 n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫 I型滤波器的最小阶和截止频率Wn。 cheblord函数的用法为： n，Wn=cheblord（Wp,Ws，Rp，Rs） 其中 Wp和 Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为 1时代表采样频率的一半.Rp和 Rs分别是通带和阻带区的波纹系数。 例1 选择设计 IIR的 Butterworth低通滤波器，其 Fs=22050Hz,Fp1=3400Hz, Fs1=5000Hz,Rp=2dB,Rs=20dB Fs=22050;Fp1=3400;Fs1=5000；Rp=3；Rs=20；设计指标 wp1=2*Fp1 /Fs；ws1=2Fs1 /Fs；求归一化频率 确定 butterworth 的最小阶数 N 和频率参数 Wn n，Wn=buttord（wp1，ws1，Rp，Rs); B，A = butter(N，Wn)；%确定传递函数的分子、分母系数 h，f=freqz（b，a，Nn,Fs_value）；%生成频率响应参数 plot（f，20*log(abs（h）) %画幅频响应图 plot(f,angle（h)； 画相频响应图 %N, Wn = buttord（Wp, Ws， Rp, Rs) 确定 butterworth 的 N 和 Wn %N， Wn = cheblord （ (Wp, Ws， Rp, Rs） 确定 Chebyshev滤波器的 N 和 Wn N, Wn = cheb2ord (Wp, Ws, Rp， Rs） 确定 Chebyshev2滤波器的 N 和 Wn %N， Wn = ellipord (Wp， Ws， Rp, Rs）确定椭圆(Ellipse) 滤波器的 N 和 Wn %B，A = butter（N,Wn，type)设计type型巴特沃斯（Butterworth)滤波器 filter. %B，A = cheby1 (N，R,Wn, type）设计type型切比雪夫滤波器 filter。 %B，A = cheby2(N,R,Wn, type)设计type型切比雪夫滤波器 filter. %B,A = ellip(N,Rp，Rs,Wn， type)设计type 型椭圆 filter. 例子 2 实现了对频率为 20和 200Hz单频叠加 cos信号的低通滤波，使输出仅含有 20Hz分量 clear; fs=1200； %采样频率 N=1200; N/fs 秒数据 n=0:N1； t=n/fs; %时间 fL=20; fH=200； sL=cos(2*pi*fL*t)； sH=cos(2*pi*fH*t）； s=sL+sH; s_in_f=fft（s)； i=1：250； plot(i，s_in_f(i); figure（1）； plot（t,s)； title（输入信号）; xlabel(t/s); ylabel(幅度)； 设计低通滤波器: Wp = 50/fs； Ws = 100/fs； 截止频率50Hz,阻带截止频率100Hz，采样频率 200Hz n，Wn = buttord（Wp,Ws,1，50); %阻带衰减大于50db,通带纹波小于 1db %估算得到 Butterworth低通滤波器的最小阶数 N和 3dB截止频率 Wn a,b=butter(n，Wn）; %设计 Butterworth低通滤波器 h，f=freqz(a,b,whole，fs); %求数字低通滤波器的频率响应 f=（0:length（f)1）fs/length(f)； 进行对应的频率转换 figure(2）； plot（f，abs(h)； 绘制 Butterworth低通滤波器的幅频响应图 title(巴氏低通滤波器）; grid; sF=filter（a，b,s）; %叠加函数 s经过低通滤波器以后的新函数 figure（3）; plot（t,sF)； %绘制叠加函数 S经过低通滤波器以后的时域图形 xlabel(t/s）； ylabel（幅度）； SF=fft（sF,N)； 对叠加函数S经过低通滤波器以后的新函数进行256点的基-2快速傅立叶变换 mag=abs（SF）; 求幅值 f=(0：length（SF)-1)*fs/length（SF)； %进行对应的频率转换 figure(4）； plot（f，mag)； 绘制叠加函数 S经过低通滤波器以后的频谱图 title(低通滤波后的频谱图）; 4 窗函数法 FIR滤波器设计实验 FIR滤波器的设计任务是选择有限长度的h(n)。使传输函数H( )满足技术要求。FIR滤波器的设计方法有多种，如窗函数法、频率采样法及其它各种优化设计方法，本实验介绍窗函数法的 FIR滤波器设计。 窗函数法是使用矩形窗、三角窗、巴特利特窗、汉明窗、汉宁窗和布莱克曼窗等设计出标准响应的高通、低通、带通和带阻 FIR滤波器. 一、firl函数的使用 在 MATLAB下设计标准响应 FIR滤波器可使用 firl函数。firl函数以经典方法实现加窗线性相位 FIR滤波器设计，它可以设计出标准的低通、带通、高通和带阻滤波器。firl函数的用法为: b=firl（n,Wn,/ftype/，Window） 各个参数的含义如下： b滤波器系数.对于一个 n阶的 FIR滤波器，其 n+1个滤波器系数可表示为：b(z）=b(1）+b(2）z1+b(n+1）zn。 n滤波器阶数。 Wn截止频率，0Wn1，Wn=1对应于采样频率的一半.当设计带通和带阻滤波器时,Wn=W1 W2,W1W2. ftype当指定 ftype时，可设计高通和带阻滤波器。Ftype=high时，设计高通 FIR滤波器；ftype=stop时设计带阻 FIR滤波器。低通和带通 FIR滤波器无需输入 ftype参数。Window窗函数。窗函数的长度应等于 FIR滤波器系数个数，即阶数n+1。 二、窗函数的使用 在 MATLAB下，这些窗函数分别为： 1矩形窗：w=boxcar(n）,产生一个 n点的矩形窗函数。 2三角窗：w=triang（n)，产生一个 n点的三角窗函数。 当 n为奇数时，三角窗系数为 w（k）= 当 n为偶数时，三角窗系数为 w（k)= 3巴特利特窗:w=Bartlett(n)，产生一个 n点的巴特利特窗函数。 巴特利特窗系数为 w(k)= 巴特利特窗与三角窗非常相似.巴特利特窗在取样点1和n上总以零结束，而三角窗在这些点上并不为零。实际上，当 n为奇数时bartlett（n）的中心n2个点等效于 triang(n2)。 4汉明窗:w=hamming（n），产生一个 n点的汉明窗函数。 汉明窗系数为w(k+1)=0.540。46cos（ ） k=0，n1 5汉宁窗：w=hanning(n)，产生一个 n点的汉宁窗函数。 汉宁窗系数为w（k）=0.51cos( ) k=1,，n 6布莱克曼窗：w=Blackman(n)，产生一个 n点的布莱克曼窗函数。 布莱克曼窗系数为w(k)=0.420.5cos（2 )+0。8cos（4 ) k=1,，n 与等长度的汉明窗和汉宁窗相比，布莱克曼窗的主瓣稍宽，旁瓣稍低。 7凯泽窗:w=Kaiser（n,beta）,产生一个 n点的凯泽窗数，其中 beta为影响窗函数旁瓣的 参数，其最小的旁瓣抑制与的关系为： 0。1102（0。87） 50 = 0。5842(21)0。4+0.07886(21） 2150 0 21 增加 可使主瓣变宽,旁瓣的幅度降低。 8契比雪夫窗：w=chebwin（n,r)产生一个 n点的契比雪夫窗函数。其傅里叶变换后的旁瓣波纹低于主瓣r个db数。例子1。FIR低通滤波器截止频率为200Hz，采样 频率为1000Hz，输入信号 x（t)=sin(100t）+sin（500t)滤波，求滤波器的输出. n=1000; fc=200; fs=1000; w=2fc/fs； t=(0:1000)/fs; window=boxcar（n+1)； b=fir1(n,w，window）; b=freqz(b，1,512,fs); s = sin(100*pit)+sin（500pit）；混 和正弦波信号 subplot（2，1,1）;plot（t，s)； sf = filter(b，1,s)；%对信号 s进行滤波 subplot(2，1，2);plot(t,sf）;3。绘制矩形窗的幅频响应,窗长分别为：N=10;N=20;N=50；N=100 分析四个窗函数窗长增加，幅频响应怎样。 wc=0。2*pi； 设一个截止频率窗函数 w1=boxcar（10)； w2=boxcar(20）； w3=boxcar（50)； w4=boxcar(100); n1=1:1:10； n2=1:1:20； n3=1:1：50； n4=1:1:100;求 h（d） hd1=sin（wc*（n-5。5)。/(pi（n-5.5）; hd2=sin（wc(n10。5)。/(pi*（n-10.5）); hd3=sin(wc*(n24。5)）./(pi（n-24。5）; hd4=sin（wc*（n-49.5）./(pi（n-49.5）； 加窗 h1=hd1。*w1； h2=hd2.*w2； h3=hd3.w3； h4=hd4.w4；求幅频特性曲线 mag1,rad=freqz(h1); mag2,rad=freqz（h2); mag3，rad=freqz(h3）； mag4，rad=freqz（h4);画幅频特性曲线 subplot(4，1，1）;plot(rad，20log10（abs（mag1)）);title（N=10);xlabel（频率 /rad）；ylabel（幅度/dB）; subplot（4，1,2）；plot(rad，20log10(abs（mag2）；title(N=20);xlabel（频率 /rad);ylabel(幅度/dB); subplot（4,1,3);plot(rad，20log10(abs(mag3）)）;title(N=50);xlabel(频率 /rad);ylabel（幅度/dB)； subplot(4,1,4）；plot(rad，20log10(abs(mag4）)）；title（N=100);xlabel（频率 /rad);ylabel(幅度/dB)；例子 2 1.含有5Hz、15Hz和 30Hz的 混和正弦波信号，设计一个 FIR带通滤波器,参数要求:采样频率fs=100Hz，通带下限截止频 率fc1=10Hz,通带上限截止频率fc2=20Hz，过渡带宽6Hz， 通阻带波动0.01,采用凯塞窗设计。 n,Wn,beta，ftype=kaiserord（7 13 17 23，0 1 0,0。01 0.01 0。01，100）; %得出滤波器的阶数n=38，beta=3.4fc1=10； fc2=20； fs=100; w1=2*fc1/fs； w2=2fc2/fs；将模拟滤波器的技术指标转换为数字滤波器的技术指标 window=kaiser(n+1,beta）;%使用 kaiser窗函数 b=fir1（n，w1 w2,window）;%使用标准频率响应的加窗设计函数 fir1 freqz（b，1，512）;数字滤波器频率响应 t = （0：100）/fs； s = sin（2*pi*t*5）+sin(2pit*15)+sin(2pi*t30)；%混和正弦波信号 subplot(2,1，1)；plot(t,s）;title（混合信号); sf = filter（b,1,s);%对信号 s进行滤波 subplot（2,1,2）；plot（t，sf）；title(滤波后信号）；xlabel（时间)；ylabel（幅度);2.设计一低通滤波器，通带为 0至1kHz， 阻带从 1500Hz到4000Hz，通带允差5，阻带波纹允差1, 阻带衰减40dB。 fs=4000; dev=0。05 0.01; fedge=1000 1500； N=200; n,Wn，beta,ftype=kaiserord(fedge，1 0,dev，fs); window=kaiser（n+1,beta); b=fir1(n,1000 1500/fs,window)；使 用标准频率响应的加窗设计函数 fir1 freqz(b,1，512）;%数字滤波器频率响应 4 FDATool设计法 FDATool(Filter Design & Analysis Tool)是 MATLAB信号处理工具箱专用的滤波器设计分析工具，操作简单、灵活，可以采用多种方法设计 FIR和 IIR滤波器.在 MATLAB命令窗口输入 FDATool后回车就会弹出 FDATool界面。 4。1 带通滤波器设计 已知滤波器的阶数n=38，beta=3。4。本例中，首先在 Filter Type中选择Bandpass；在 Design Method选项中选择FIR Window，接着在 Window选项中选取Kaiser,Beta值为3。4；指定 Filter Order项中的 Specify order为38；采样频率Fs=100Hz,截止频率 Fc1=10Hz,Fc2=20Hz。设置完以后点击窗口下方的Design Filter，在窗口上方就会看到所设计滤波器的幅频响应，通过菜单选项 Analysis还可以看到滤波器的相频响应、组延迟、脉冲响应、阶跃响应、零极点配置等.设计完成后将结果保存为 kaiser15。fda文件。 4。2 Simulink仿真 在 Simulink环境下，将滤波器文件 kaiser15.fda导入 Digital Filter Design模块,输入信号为s（t）=sin(10t）+sin(30t)+sin(60t),生成的仿真图和滤波效果如图2所示。 （1）Simulink仿真图 （2）滤波前后的离散波形 图 2 Simulink仿真图和滤波效果图 5 SPTool设计法 SPTool是 MATLAB信号处理工具箱中自带的交互式图形用户界面工具，它包含了信号处理工具箱中的大部分函数,可以方便快捷地完成对信号、滤波器及频谱的分析、设计和浏览.在本例中按以下步骤完成滤波器的设计和滤波： 创建并导入信号源。 在 MATLAB命令窗口输入命令： Fs100；t = (0：100）/Fs； s = sin（2*pit*5）+sin(2*pit*15)+sin（2pit*30）； 此时，变量Fs、t、s将显示在 workspace列表中。在命令窗口键入Sptool，将弹出 Sptool主界面,如图 3所示；点击菜单 File/Import将信号 s导入并取名为s。 （2)单击 Filters列表下的New，按照参数要求设计出滤波器filt1，具体步骤类似于3