matlab信号处理学习总结_第1页
matlab信号处理学习总结_第2页
matlab信号处理学习总结_第3页
matlab信号处理学习总结_第4页
matlab信号处理学习总结_第5页
已阅读5页,还剩26页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、常用函数1图形化信号处理工具,fdatool(滤波器设计),fvtool(图形化滤波器参数查看)sptool(信号处理),fvtool(b,a) , wintool窗函数设计.或者使用工具箱filter design 设计。当使用离散的福利叶变换方法分析频域中的信号时,傅里叶变换时可能引起漏谱,因此 需要采用平滑窗,2数字滤波器和采样频率的关系。如果一个数字滤波器的采样率为FS,那么这个滤波器的分析带宽为Fs/2。也就是说这个滤波器只可以分析0, Fs/2的信号.举个例字:有两个信号,S1频率为20KHz, S2频率为40KHz,要通过数字方法滤除S2。你的滤波器的采样率至少要为Fs=80HK

2、z,否则就分析不到S2 了,更不可能将它滤掉了!(当然根据采样定理,你的采样率F0也必须大于80HK, ,Fs和F0之间没关系不大,可以任取,只要满足上述关系就行。)3两组数据的相关性分析r=corrcoef(x,y)4 expm求矩阵的整体的 exp4离散快速傅里叶fft信号处理中,傅里叶变换的典型用途是将信号分解成幅值分量和频 率分量)。Ft为连续傅里叶变换。反傅里叶 ifft5 ztrans () , Z变换是把离散的数字信号从时域转为频率6 laplace ()拉普拉斯变换是把连续的的信号从时域转为频域7 sound(x)会在音响里产生 x所对应的声音8 norm求范数,det行列式,

3、rank求秩9模拟频率,数字频率,模拟角频率关系模拟频率f :每秒经历多少个周期,单位Hz,即1/s ;模拟角频率是指每秒经历多少弧度,单位 rad/s ;数字频率w:每个采样点间隔之间的弧度,单位 rad。Q =2pi*f ; w = Q *T10 RMS求法Rms = sqrt(sum(P.A2) 或者 norm(x)/sqrt(length(x) var方差的开方是 std 标准差,RMW该是 norm(x)/sqrt(length(x) 吧.求矩阵的 RMS std(A(:) 11 ftshift 作用:将零频点移到频谱的中间 12 filtfilt零相位滤波,采用两次滤波消除系统的非

4、线性相位,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

5、(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为,0-1之间。Ws为阻带角频率,0-1之间。Rp为通带波纹,或者通带衰减,Rs为阻带衰减。若果给出的是模拟频率fp1通带截止频率,fp2阻带截止频率,则Wp=fp1*2/fs,Ws

6、=fp2*2/fs.如果给出的是实际数字频率比如0.3*pi,如果给出的是y=filter(b,a,x);或者采用零相位滤波y=filfilt(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函数调用,窗函数长度

7、为滤波器最佳阶数 n+117.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/Delta

8、W=6.6pi/0.04pi=165所以滤波器的阶数至少是165。在该帖子中是用理想低通滤波器的方法来计算的,这里用fir1函数来计算,相应的程序有fs=100; % 采样频率wp = 3*pi/50; ws = 5*pi/50;deltaw= ws - wp; %过渡带宽 A w的计算N = ceil(6.6*pi/ deltaw) + 1; %按海明窗计算所需的滤波器阶数N0wdham = (hamming(N+1)' %海明窗计算Wn=(3+5)/100; % 计算截止频率b=fir1(N,Wn,wdham);H,w=freqz(b,1);db=20*log10(abs(H);%

9、画频响曲线plot(w*fs/(2*pi),db);title('幅度响应(单位:dB ) ');gridaxis(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

10、,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')'f

11、type '对应'high'是高通滤波器的归一化截止频率'low'低通滤波器的归一化截止频率'stop' for an order 2*n bandstop digital filter if Wn is a two-elementvector, Wn = w1 w2. The stopband is w1 <w < w2.21窗函数1 矩形窗:Window=boxcar(8);b=fir1(7,0.4,Window);freqz(b,1)2 blackman 窗:Window=blackman(8);b=fir1(7,0.4

12、,Window);freqz(b,1)3 hamming;4 hanning;5 kaiser窗函数第一旁瓣相对于主瓣衰减/dB主瓣宽阻带最小衰减/dB矩形窗-13 4兀/N21三角窗-25 8兀/N25汉宁窗-31 8兀/N44海明窗-41 8兀/N53布拉克曼窗-57 12 兀 /N 74凯塞窗可调可调可调切比雪夫窗可调可调可调15.1 基于firpm 函数的最佳fir 滤波器设计例子f和a长度相同,且长度为偶数Graph the desired and actual frequency responses of a 17th-order Parks-McClellan bandpass

13、filter:f = 0 0.3 0.4 0.6 0.7 1; a = 0 0 1 1 0 0;b = firpm(17,f,a);h,w = freqz(b,1,512);plot(f,a,w/pi,abs(h)legend('Ideal','firpm Design')15.2 最佳FIR滤波器阶数估计n,fo,ao,w = firpmord(f,a,dev)n,fo,ao,w = firpmord(f,a,dev,fs)例子2设计一个最低阶低通滤波器通带截止频率500Hz,阻带截止频率6000Hz,采样率2000Hz,通带波纹小于3dB,阻带衰减大于40d

14、B.rp = 3; % Passband ripplers = 40; % Stopband ripplefs = 2000; % Sampling frequencyf = 500 600; % Cutoff frequenciesa = 1 0; % Desired amplitudesdev = (10A(rp/20)-1)/(10A(rp/20)+1) 10A(-rs/20);n,fo,ao,w = firpmord(f,a,dev,fs);b = firpm(n,fo,ao,w);freqz(b,1,1024,fs);18 y=resample(x,p,q)数字信号中的重采样。这时输

15、出信号y的采样频率就是x的p/q倍,其长度为length(x)*p/q19 conv 卷积deconv 反卷积或者求多项式乘法。xcorr互相关函数 cov协方差fft2 二维FFT fft2 二维FFT逆变换 xcorr2 , conv2 二维卷积 20平滑滤波filter 函数首先要设计好滤波器,然后调用filter .平滑滤波似乎有些过时,butterworth才显得稍微有些技术含量用法。filter本身作用是求卷积和convfilter(B,1,X,dim) ; dim缺省为1,是按列滤波的,如果改为 2,则是按行滤波。y = filter(b, a, x),其中b,a为滤波器系数。计

16、算系统在输入x作用下的零状态响应yk举例:计算低通滤波器的冲激响应 例题1点平均滤波f1=3;f2=40;fs=100;t=0:1/fs:1;x=sin(2*pi*t*f1)+.25*sin(2*pi*t*f2);b=ones(1,10)/10; y=filter(b,1,x);求冲激响应stem(y);yy=filtfilt(b,1,x);plot(t,x);hold on;plot(t,x,'r',t,yy,'g')例2利用filter 函数求滑动平均Matlab有多种计算滑动平均的方法,现介绍基于filter函数的计算方法。设原始数据为x,平均窗口设为a

17、(a为正整数),那么无权重滑动平均后的数据y为:windowSize = a;y=filter(ones(1,windowSize)/windowSize,1,x);上述命令实际上计算的是:y(1)=(1/a)*x(1);y(2)=(1/a)*x(2)+(1/a)*x(1);y(a)=(1/a)*x(a)+(1/a)*x(a-1)+.+(1/a)*x(1);y(i)=(1/a)*x(i)+(1/a)*x(i-1)+.+(1/a)*x(i-a+1);4. frezq数字滤波器的频率响应H, W=freqz(B, A, N) 当N是一个整数时返回N点的频率向量 H和N点的幅频响应向量 W N最好选

18、用2的整数次哥,这样便于使用FFT进行快速算法。H为滤波器的复数放大倍数,w为频率向量,如果只想获得放大倍数的幅值,可以用plot(w,abs(h).如果是滤波器设计 plot(w/pi,abs(h)滤波器放大倍数,低频时为1,高频时为0,即低通滤波器N个频率点均匀地分布在单位圆的上半圆上。如果 N没有确定则却缺省为512个点。freqz(B, A, N)将直接绘制频率响应图,而不返回任何值。H, W=freqz(B, A, N, 'whole')运用分布在整个单位圆上的N个点。H=freqz(B, A, W)返回指定在 W向量中频率范围的频率响应,其中W是以弧度为单位在0,

19、pi 范围内。H,F=freqz(B, A, N, Fs) , H,F=freqz(B, A, N, Fs)这两个函数给出了采样频率Fs,则返回频率向量F,它们的单位都是 Hz。invfreqz()是其逆函数,它运用最小二乘法从已知的频率响应中求出传递函数模型。例子2滤波器的特性分析二3.数字滤波器的冲击响应:impz()H, T=impz(B, A)返回滤波器的冲击响应列向量H和时间即采样间隔列向量T。滤波器是用传递函数模型来限定的。impz(B, A)将直接绘制滤波器的冲击响应图。例:设计一个高通 Chebyshev2型滤波器,它的具体要求是阻带的截止频率为10Hz,通带的截止频率为30H

20、z,在通带内的最大衰减不超过0.1dB,在阻带内的最小衰减不小于40dB。程序设计如下:wp=30; ws=10; rp=0.1; rs=40; Fs=100;wp=10*2*pi; ws=30*2*pi;n, wn=cheb20rd(wp/(Fs/2), ws/(Fs/2), rp, rs, 's');z, p, k=cheb2ap(n, rs);wn=wn*(Fs/2)*2*pi;A, B, C, D=zp2ss(z,p,k);AT, BT, CT, DT=lp2hp(A, B, C, D, wn);AT1, BT1, CT1, DT1=bilinear(AT, BT, C

21、T, DT, Fs);num, den尸ss2tf(AT1, BT1, CT1, DT1);H, W=freqz(num, den); figure;plot(W*Fs/(2*pi), abs(H); grid;xlabel('频率/Hz');ylabel('幅值'); impz(num, den);一、巴特沃斯IIR滤波器的设计1首先根据滤波器设计要求用buttord 函数求出最小滤波器阶数和截止频率n,Wn= buttord(Wp,Ws,Rp,Rs)其中Wpm Ws分别是通带和阻带的拐角频率(截止频率),其取值范围为 0至1之间。当其值为1时代表采样频率的一

22、半。Rp和Rs分别是通带和阻带区的波纹系数。不同类型(高通、低通、带通和带阻)滤波器对应的Wp和Ws值遵循以下规则:1 .高通滤波器:Wp和Ws为一元矢量且Wp>Ws2 .低通滤波器:Wp和Ws为一元矢量且 Wp<Ws3 .带通滤波器:Wp和 Ws 为二元矢量且 Wp<Ws 如 Wp=0.2,0.7,Ws=0.1,0.8;4 .带阻滤波器:Wp和 Ws 为二元矢量且 Wp>Ws 如 Wp=0.1,0.8,Ws=0.2,0.7。2求出滤波器系数Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器,其特性为使通带内的幅度响应最大限度地平坦,但同时损失截止频率

23、处的下降斜度。在期望通带平滑的情况 下,可使用butter 函数。butter函数的用法为:b,a=butter(n,Wn,ftype) 。其中n代表滤波器阶数,Wnf弋表滤波器的截止频率二、契比雪夫I型IIR滤波器的设计在期望通带下降斜率大的场合,应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用chebyl函数设计出契比雪夫I型IIR滤波器。chebyl函数可设计低通、高通、带通和带阻契比雪夫 I型滤IIR波器,其通带内为等波纹, 阻带内为单调。契比雪夫 I型的下降斜度比II型大,但其代价是通带内波纹较大。cheby1函数的用法为:b,a=cheby1(n,Rp,Wn,/ftype/

24、)在使用cheby1函数设计IIR滤波器之前,可使用 cheblord 函数求出滤波器阶数 n和截止 频率Wn cheblord函数可在给定滤波器性能的情况下,选择契比雪夫 I型滤波器的最小阶 和截止频率Wncheblord函数的用法为: 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=20

25、dBFs=22050; Fp1=3400; Fs1=5000; Rp=3; Rs=20;流计指标wp1=2*Fp1/Fs ; ws1=2*Fs1/Fs ; %t归一化频率%确定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,

26、 Rs) %N, Wn = cheblord ( (Wp, Ws, Rp, Rs) %N, Wn = cheb2ord (Wp, Ws, Rp, Rs) %N, Wn = ellipord (Wp, Ws, Rp, Rs) %B,A = butter(N,Wn,'type') %B,A = chebyl (N,R,Wn, 'type') %B,A = cheby2(N,R,Wn, 'type') %B,A = ellip(N,Rp,Rs,Wn, 'type') 例子2设计'type'确定 butterworth 的

27、 N 和 Wn确定 Chebyshev滤波器的 N 和 Wn确定 Chebyshev2滤波器的 N和 Wn确定椭圆(Ellipse)滤波器的N和Wn型巴特沃斯(Butterworth) 滤波器filter.设计type'型切比雪夫I滤波器 filter. 设计type'型切比雪夫n滤波器 filter. 设计type' 型椭圆filter.现了对频率为20和200Hz单频叠加cos信号的低通滤波,使输出仅含有20Hz分量clear;fs=1200; %采样频率 N=1200; % N/fs 秒数据 n=0:N-1; t=n/fs; % 时间 fL=20;fH=200;s

28、L=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。); 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,通带纹波小于

29、 1db%古算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wna,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经过低通滤波器以后的新函数

30、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滤波器设计实验 F

31、IR滤波器的设计任务是选择有限长度的h(n)。使传输函数 H()满足技术要求。FIR滤波器的设计方法有多种,如窗函数法、频率采样法及其它各种优化设计方法,本实验介绍窗 函数法的FIR滤波器设计。窗函数法是使用矩形窗、三角窗、巴特利特窗、汉明窗、汉宁窗和布莱克曼窗等设计出标 准响应的高通、低通、带通和带阻FIR滤波器。一、firl函数的使用在MATLAB下设计标准响应 FIR滤波器可使用firl 函数。firl函数以经典方法实现加窗线 性相位FIR滤波器设计,它可以设计出标准的低通、带通、高通和带阻滤波器。firl函数的用法为:b=firl(n,Wn,/ftype/,Window) 各个参数的含

32、义如下:b一滤波器系数。对于一个n阶的FIR滤波器,其n+1个滤波器系数可表示为:b亿尸b(1)+b(2)z 1 +b(n+1)z n。n一滤波器阶数。Wn-截止频率,0wWnc 1, Wn=1对应于采样频率的一半。当设计带通和带阻滤波器时,Wn=W1W2,W1< co < W2 ftype 当指定ftype 时,可设计高通和带阻滤波器。Ftype=high时,设计高通 FIR滤波器;ftype=stop 时设计带阻FIR滤波器。低通和带通 FIR滤波器无需输入ftype 参数。Window-窗函数。窗函数的长度应等于FIR滤波器系数个数,即阶数 n+1。二、窗函数的使用在MATL

33、AB下,这些窗函数分别为: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) 的中心n- 2个点等效于triang(n 2)。4 .汉明窗:w=hamming(n),产生一个n点的汉明窗函数。汉明窗系数为 w

34、(k+1)=0.54 0.46cos ( ) k=0,n15 .汉宁窗:w=hanning(n),产生一个 n点的汉宁窗函数。汉宁窗系数为 w(k)=0.51 cos( ) k=1,n6 .布莱克曼窗:w=Blackman(n),产生一个n点的布莱克曼窗函数。布莱克曼窗系数为 w(k)=0.42 0.5cos(2 兀)+0.8cos(4兀)k=1,n和等长度的汉明窗和汉宁窗相比,布莱克曼窗的主瓣稍宽,旁瓣稍低。7 .凯泽窗:w=Kaiser(n,beta),产生一个n点的凯泽窗数,其中 beta为影响窗函数旁瓣 的3参数,其最小的旁瓣抑制”和3的关系为:0.1102( a 0.87) a &g

35、t;503 = 0.5842( a 21)0.4+0.07886( a 21) 21 < a < 50 0 a <21增加3可使主瓣变宽,旁瓣的幅度降低。8 .契比雪夫窗:w=chebwin(n,r)产生一个n点的契比雪夫窗函数。其傅里叶变换后的旁瓣 波纹低于主瓣r个db数。例子1.FIR低通滤波器截止频率为 200Hz,采样频率为1000Hz,输入信号 x(t)=sin(100 兀t)+sin(500 兀t)滤波,求滤波器的输出 .n=1000;fc=200;fs=1000;w=2*fc/fs;t=(0:1000)/fs;window=boxcar(n+1);b=fir1(

36、n,w,window);b=freqz(b,1,512,fs);s = sin(100*pi*t)+sin(500*pi*t);%混 和正弦波信号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

37、;n2=1:1:20;n3=1:1:50;n4=1:1:100;W h(d)hd1=sin(wc*(n-5.5)./(pi*(n-5.5);hd2=sin(wc*(n-10.5)./(pi*(n-10.5);hd3=sin(wc*(n-24.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);

38、mag4,rad=freqz(h4);%画幅频特性曲线subplot(4,1,1);plot(rad,20*log10(abs(mag1);title('N=10');xlabel('频率/rad');ylabel('幅度/dB');subplot(4,1,2);plot(rad,20*log10(abs(mag2);title('N=20');xlabel('频率/rad');ylabel('幅度/dB');subplot(4,1,3);plot(rad,20*log10(abs(mag3);ti

39、tle('N=50');xlabel('频率/rad');ylabel('幅度/dB');subplot(4,1,4);plot(rad,20*log10(abs(mag4);title('N=100');xlabel('频率/rad');ylabel('幅度/dB');例子21.含有5Hz、15Hz和30Hz的 混和正弦波信号,设计一个FIR带通滤波器,参数要求:采样频率fs=100Hz,通带下限截止频率fc1=10Hz,通带上限截止频率fc2=20Hz,过渡带宽6Hz,通阻带波动0.01 ,采用

40、凯塞窗设计。n,Wn,beta,ftype=kaiserord(7 13 17 23,0 1 0,0.01 0.01 0.01,100);%得出滤波器的阶数n=38,fc1=10;fc2=20;fs=100;w1=2*fc1/fs; w2=2*fc2/fs;% window=kaiser(n+1,beta);% b=fir1(n,w1 w2,window);%beta=3.4将模拟滤波器的技术指标转换为数字滤波器的技术指标使用kaiser窗函数使用标准频率响应的加窗设计函数fir1freqz(b,1,512);%数字滤波器频率响应t = (0:100)/fs;混和正弦波信号subplot(2,

41、1,1);plot(t,s);title('sf = filter(b,1,s);%对信号subplot(2,1,2);plot(t,sf);title('2.设计一低通滤波器,通带为 纹允差1%阻带衰减40dBofs=4000;混合信号);s进行滤波滤波后信号);xlabel(' 时间');ylabel(' 幅度);0至1kHz,阻带从1500Hz到4000Hz,通带允差5%阻带波s = sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30);%dev=0.05 0.01;fedge=1000 1500;N=200;n,

42、Wn,beta,ftype=kaiserord(fedge,1 0,dev,fs);window=kaiser(n+1,beta);b=fir1(n,1000 1500/fs,window);% 使 用标准频率响应的加窗设计函数fir1freqz(b,1,512);%数字滤波器频率响应4 FDATool设计法FDATool(Filter Design & Analysis Tool)是MATLAB信号处理工具箱专用的滤波器设计分析工具,操作简单、灵活,可以采用多种方法设计FIR和IIR 滤波器。在MATLAB命令窗口输入FDATool后回车就会弹出 FDATool界面。4.1 带通滤波

43、器设计已知滤波器的阶数 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还可以看到滤波器的相频响应、组延迟、脉冲响应、阶跃响应、零极点

44、配置等。设计完成后将结果保存为kaiser15.fda 文件。4.2 Simulink 仿真在Simulink 环境下,将滤波器文件kaiser15.fda 导入Digital Filter Design 模块,输入信号为s(t尸sin(10兀t)+sin(30兀t)+sin(60兀t),生成的仿真图和滤波效果如图 2所 示。(1) Simulink仿真图 (2)滤波前后的离散波形图2 Simulink 仿真图和滤波效果图5 SPTool设计法SPTool是MATLAB信号处理工具箱中自带的交互式图形用户界面工具,它包含了信号 处理工具箱中的大部分函数,可以方便快捷地完成对信号、滤波器及频谱的

45、分析、设计和 浏览。在本例中按以下步骤完成滤波器的设计和滤波: 创建并导入信号源。在MATLAB命令窗口输入命令:Fs=100; t = (0:100)/Fs;s = sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30);此时,变量Fs、t、s将显示在 workspace列表中。在命令窗口键入 Sptool ,将弹出Sptool 主界面,如图3所示;点击菜单File/Import 将信号s导入并取名为s。(2)单击Filters 列表下的New,按照参数要求设计出滤波器filt1 ,具体步骤类似于 3.2.1 。(3)将滤波器filt1 使用到s信号序列。分别

46、在 Signals、Filters 、Spectra列表中 选择s、filt1 、mtlbse ,单击Filters列表下的 Apply 按钮,在弹出的 Apply Filter 对话框中将输出彳t号命名为sin15hz。(4)进行频谱分析。在 Signals 中选择s,单击Spectra 下的Create按钮,在弹出的Spectra Viewer界面中选择 Method为FFT, Nfft=512 ,单击 Apply 按钮生成 s的频谱spect1 。同样的步骤可以生成信号 sin15hz的频谱spect2。分别选中信号s、sin15hz、spect1、spect2 ,单击各自列表下方的Vi

47、ew按钮,即可观察他们的波形,如图 4所示。图3 SPTool主界面 图4滤波前后的时域波形和频域特性由图4可以看出,带通滤波器filt1 使输入信号s中频率为15hz的正弦波信号通过,而将频率为5hz和30hz的正弦波信号大大衰减。matlab语音信号处理fs=22050; %语音信号采样频率为22050x1=wavread('Windows Critical Stop.wav'); %读取语音信号的数据,赋给变量x1sound(x1,22050); % 播放语音信号y1=fft(x1,1024); % 对信号做 1024 点 FFT 变换f=fs*(0:511)/1024;

48、figure(1)plot(x1) %做原始语音信号的时域图形title('原始语音信号);xlabel('time n');ylabel('fuzhi n');figure(2)freqz(xl) %绘制原始语音信号的频率响应图title('频率响应图')figure(3)subplot(2,1,1);plot(abs(y1(1:512) %做原始语音信号的FFT频谱图title(' 原始语音信号FFT频谱')subplot(2,1,2);plot(f,abs(y1(1:512);title('原始语音信号频谱)

49、xlabel('Hz');ylabel('fuzhi');程序2:fs=22050; %语音信号采样频率为22050x1=wavread('Windows Critical Stop.wav'); % 读取语音信号的数据,赋给变量x1t=0:1/22050:(size(x1)-1)/22050;y1=fft(x1,1024); %对信号做 1024 点 FFT 变换f=fs*(0:511)/1024;x2=randn(1,length(x1); %产生一和x长度一致的随机信号sound(x2,22050);figureplot(x2) %做原始语

50、音信号的时域图形title(' 高斯随机噪声);xlabel('time n');ylabel('fuzhi n');randn('state',0);m=randn(size(x1);x2=0.1*m+x1;sound(x2,22050);%播放加噪声后的语音信号y2=fft(x2,1024);figure(2) plot(t,x2) title('加噪后的语音信号);xlabel('time n');ylabel('fuzhi n'); figure(3) subplot(2,1,1);plot

51、(f,abs(y2(1:512);title('原始语音信号频谱);xlabel('Hz');ylabel('fuzhi');subplot(2,1,2);plot(f,abs(y2(1:512);title('加噪后的语音信号频谱');xlabel('Hz');ylabel('fuzhi');根据以上代码,你可以修改下面有错误的代码程序3:双线性变换法设计Butterworth 滤波器fs=22050;课程设计 2shuzi.wav');t=0:1/22050:(size(x1)-1)/22050

52、;Au=0.03;d=Au*cos(2*pi*5000*t)'x2=x1+d;wp=0.25*pi;ws=0.3*pi;Rp=1;Rs=15;Fs=22050;Ts=1/Fs;wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标ws1=2/Ts*tan(ws/2);N,Wn=buttord(wp1,ws1,Rp,Rs,'s'); % 选择滤波器的最小阶数 Z,P,K=buttap(N); % 创建 butterworth模拟滤波器Bap,Aap=zp2tf(Z,P,K);b,a=lp21P(Bap,Aap,Wn);bz,az=bi1inear(b,a,Fs

53、); %用双线性变换法实现模拟滤波器到数字滤波器的转换 H,W=freqz(bz,az); %绘制频率响应曲线figure(1) plot(W*Fs/(2*pi),abs(H) grid xlabel('频率/ Hz') ylabel(' 频率响应幅度') title('Butterworth') f1=fi1ter(bz,az,x2);figure(2)subplot(2,1,1)plot(t,x2) %画出滤波前的时域图title(' 滤波前的时域波形');subplot(2,1,2)plot(t,f1); %画出滤波后的时域

54、图title(' 滤波后的时域波形');sound(f1,22050); %播放滤波后的信号F0=fft(f1,1024);f=fs*(0:511)/1024;figure(3)y2=fft(x2,1024);subplot(2,1,1);plot(f,abs(y2(1:512); % 画出滤波前的频谱图 title('滤波前的频谱')xlabel('Hz'); ylabel('fuzhi');subplot(2,1,2)F1=plot(f,abs(F0(1:512); %画出滤波后的频谱图title(' 滤波后的频谱)

55、xlabel('Hz'); ylabel('fuzhi');程序4:窗函数法设计滤波器: fs=22050;课程设计 2shuzi.wav');t=0:1/22050:(size(x1)-1)/22050;Au=0.03;d=Au*cos(2*pi*5000*t)'x2=x1+d;wp=0.25*pi;ws=0.3*pi;wdelta=ws-wp;N=ceil(6.6*pi/wdelta);wn=(0.2+0.3)*pi/2;b=fir1(N,wn/pi,hamming(N+1);figure 现整 %选择窗函数,并归一化截止频率 freqz(b

56、,1,512) f2=filter(bz,az,x2)figure(2)subplot(2,1,1)plot(t,x2)title('滤波前的时域波形);subplot(2,1,2)plot(t,f2);title('滤波后的时域波形');sound(f2,22050); %播放滤波后的语音信号F0=fft(f2,1024);f=fs*(0:511)/1024;figure(3)y2=fft(x2,1024);subplot(2,1,1);plot(f,abs(y2(1:512);title(' 滤波前的频谱')xlabel('Hz');

57、ylabel('fuzhi');subplot(2,1,2)F2=plot(f,abs(F0(1:512);title(' 滤波后的频谱') xlabel('Hz');ylabel('fuzhi');1:使用wavrecord 函数2:使用awgn等函数3:使用fdatool可以设计滤波器5: wavplay6: guidematlab滤波matlab滤波的方法包括平滑滤波,fir , irr,小波滤波,小波包滤波,自适应lms, rls滤波,最佳fir 滤波,fadtool ,或者使用filter design工具箱。filter 适合于平滑滤波,对于无特殊要求的滤波。y = filter(b, a, x) ,其中b,a为滤波器系数。计算系统在输入 x作用下的零状态响应 yk。可以结合filter 函数和fir滤波或 者irr滤波。b,a可以使fir或者irr产生的滤波器系数。fir滤波线性相位滤波器

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论