matlab信号处理学习总结.docx
- 文档编号:17307974
- 上传时间:2023-07-24
- 格式:DOCX
- 页数:56
- 大小:41.26KB
matlab信号处理学习总结.docx
《matlab信号处理学习总结.docx》由会员分享,可在线阅读,更多相关《matlab信号处理学习总结.docx(56页珍藏版)》请在冰点文库上搜索。
matlab信号处理学习总结
常用函数
1图形化信号处理工具,fdatool(滤波器设计),fvtool(图形化滤波器参数查看)sptool(信号处理),fvtool(b,a),wintool窗函数设计.或者使用工具箱filterdesign设计。
当使用离散的福利叶变换方法分析频域中的信号时,傅里叶变换时可能引起漏谱,因此
需要采用平滑窗,
2数字滤波器和采样频率的关系。
如果一个数字滤波器的采样率为FS,那么这个滤波器的分析带宽为Fs/2。
也就是说这
个滤波器只可以分析[0,Fs/2]的信号.举个例字:
有两个信号,S1频率为20KHz,S2频率为40KHz,要通过数字方法滤除S2。
你的滤波器的采样率至少要为Fs=80HKz,否则就分析不到S2了,更不可能将它滤掉
了!
(当然根据采样定理,你的采样率F0也必须大于80HK,,Fs和F0之间没关系不大,可
以任取,只要满足上述关系就行。
)
3两组数据的相关性分析r=corrcoef(x,y)
4expm求矩阵的整体的exp
4离散快速傅里叶fft信号处理中,傅里叶变换的典型用途是将信号分解成幅值分量和频
率分量)。
Ft为连续傅里叶变换。
反傅里叶ifft
5ztrans(),Z变换是把离散的数字信号从时域转为频率
6laplace()拉普拉斯变换是把连续的的信号从时域转为频域
7sound(x)会在音响里产生x所对应的声音
8norm求范数,det行列式,rank求秩
9模拟频率,数字频率,模拟角频率关系
模拟频率f:
每秒经历多少个周期,单位Hz,即1/s;
模拟角频率Ω是指每秒经历多少弧度,单位rad/s;
数字频率w:
每个采样点间隔之间的弧度,单位rad。
Ω=2pi*f;w=Ω*T
10RMS求法
Rms=sqrt(sum(P.^2))或者norm(x)/sqrt(length(x)var方差的开方是std标准差,
RMS应该是norm(x)/sqrt(length(x))吧.求矩阵的RMS:
std(A(:
))
11ftshift作用:
将零频点移到频谱的中间
12filtfilt零相位滤波,
采用两次滤波消除系统的非线性相位,
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角频率,
常用于查看滤波器的频率特性
14zplane(z,p)画图零极点分布图
15beta=unwarp(alpha)相位会在穿越+-180发生回绕,可将回绕的
16stepz求数字滤波器的阶跃响应
[h,t]=stepz(b,a,n,fs)
fvtool(b1,a1,b2,a2,...bn,an)fvtool(Hd1,Hd2,...)h=fvtool(...)
15IIR数字滤波器设计方法
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=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.350.65]);
freqz(b,1,512)
16窗函数长度:
窗函数的长度应等于FIR滤波器系数个数,即滤波器阶数n+1。
17加窗函数的FIR滤波器长度的确定
17.1buttord函数求出最佳滤波器阶数和截止频率,然后用fir1函数调用,窗函数长度
为滤波器最佳阶数n+1
17.2用窗函数方法设计FIR滤波器,由滤波器的过渡带的宽度和选择的窗函数决定
这里举一个选用海明窗函数设计低通滤波器的例子。
低通滤波器的设计要求是:
采样频率为100Hz,通带截至频率为3Hz,阻带截止频率为5Hz,
通带内最大衰减不高于0.5dB,阻带最小衰减不小于50dB。
使用海明窗函数。
确定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(w*fs/(2*pi),db);title('幅度响应(单位:
dB)');grid
axis([050-10010]);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'foranorder2*nbandstopdigitalfilterifWnisatwo-element
vector,Wn=[w1w2].Thestopbandisw1<ω 21窗函数 1矩形窗: Window=boxcar(8); b=fir1(7,0.4,Window); freqz(b,1) 2blackman窗: Window=blackman(8); b=fir1(7,0.4,Window); freqz(b,1) 3hamming;4hanning;5kaiser 窗函数第一旁瓣相对于主瓣衰减 /dB 主瓣宽阻带最小衰减/dB 矩形窗–134π/N21 三角窗–258π/N25 汉宁窗–318π/N44 海明窗–418π/N53 布拉克曼 窗 –5712π/N74 凯塞窗可调可调可调 切比雪夫 窗 可调可调可调 15.1基于firpm函数的最佳fir滤波器设计 例子f和a长度相同,且长度为偶数 Graphthedesiredandactualfrequencyresponsesofa17th-orderParks-McClellan bandpassfilter: f=[00.30.40.60.71];a=[001100]; b=firpm(17,f,a); [h,w]=freqz(b,1,512); plot(f,a,w/pi,abs(h)) legend('Ideal','firpmDesign') 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,阻带衰减大于40dB. rp=3;%Passbandripple rs=40;%Stopbandripple fs=2000;%Samplingfrequency f=[500600];%Cutofffrequencies a=[10];%Desiredamplitudes dev=[(10^(rp/20)-1)/(10^(rp/20)+1)10^(-rs/20)]; [n,fo,ao,w]=firpmord(f,a,dev,fs); b=firpm(n,fo,ao,w); freqz(b,1,1024,fs); 18y=resample(x,p,q)数字信号中的重采样。 这时输出信号y的采样频率就 是x的p/q倍,其长度为length(x)*p/q 19conv卷积deconv反卷积或者求多项式乘法。 xcorr互相关函数cov协方差fft2二维FFTfft2二维FFT逆变换 xcorr2,conv2二维卷积 20平滑滤波filter函数 首先要设计好滤波器,然后调用filter.平滑滤波似乎有些过时,butterworth才显得稍微 有些技术含量用法。 filter本身作用是求卷积和conv filter(B,1,X,[],dim);dim缺省为1,是按列滤波的,如果改为2,则是按行滤波。 y=filter(b,a,x),其中b,a为滤波器系数。 计算系统在输入x作用下的零状态响应y[k] 举例: 计算低通滤波器的冲激响应 例题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);holdon; plot(t,x,'r',t,yy,'g') 例2利用filter函数求滑动平均 Matlab有多种计算滑动平均的方法,现介绍基于filter函数的计算方法。 设原始数据为x, 平均窗口设为a(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最好选用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,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,通带 的截止频率为30Hz,在通带内的最大衰减不超过0.1dB,在阻带内的最小衰减不小于40dB。 程序设计如下: wp=30;ws=10;rp=0.1;rs=40;Fs=100; wp=10*2*pi;ws=30*2*pi; [n,wn]=cheb2ord(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,CT,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) 其中Wp和Ws分别是通带和阻带的拐角频率(截止频率),其取值范围为0至1之间。 当 其值为1时代表采样频率的一半。 Rp和Rs分别是通带和阻带区的波纹系数。 不同类型(高通、低通、带通和带阻)滤波器对应的Wp和Ws值遵循以下规则: 1.高通滤波器: Wp和Ws为一元矢量且Wp>Ws; 2.低通滤波器: Wp和Ws为一元矢量且Wp 3.带通滤波器: Wp和Ws为二元矢量且Wp 4.带阻滤波器: Wp和Ws为二元矢量且Wp>Ws,如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=2*Fs1/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: N-1; 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—滤波器系
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 信号 处理 学习 总结