经典matlab信号处理基础知识.docx

经典matlab信号处理基础知识.docx

《经典matlab信号处理基础知识.docx》由会员分享，可在线阅读，更多相关《经典matlab信号处理基础知识.docx（50页珍藏版）》请在冰点文库上搜索。

经典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=.04pi

3，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-elementvector,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

矩形窗

–13

4π/N

21

三角窗

–25

8π/N

25

汉宁窗

–31

8π/N

44

海明窗

–41

8π/N

53

布拉克曼窗

–57

12π/N

74

凯塞窗

可调

可调

可调

切比雪夫窗

可调

可调

可调

15.1基于firpm函数的最佳fir滤波器设计

例子f和a长度相同，且长度为偶数

Graphthedesiredandactualfrequencyresponsesofa17th-orderParks-McClellanbandpassfilter:

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—滤波器系数。

对于一个n阶的FIR滤波器，其n+1个滤波器系数可表示为：

b（z）=b

（1）+b

（2）z－1+…+b（n+1）z－n。

n—滤波器阶数。

Wn—截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。

当设计带通和带阻滤波器时，Wn=[W1W2],W1≤ω≤W2。

ftype—当指定ftype时，可设计高通和带阻滤波器。

Ftype=high时，设计高通FIR滤波器；ftype=stop时设计带阻FIR滤波器。

低通和带