matlab信号处理学习总结

常用函数
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。

Ω=2pi*f； 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为，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.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(w*fs/(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 < ω < 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,Window);
freqz(b,1)
3 hamming;
4 hanning;
5 kaiser
窗函数第一旁瓣相对于主瓣衰减
/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长度相同，且长度为偶数
Graph the desired and actual frequency responses of a 17th-order Parks-McClellan bandpass 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,阻带衰减大于 40dB.
rp = 3; % Passband ripple
rs = 40; % Stopband ripple
fs = 2000; % Sampling frequency
f = [500 600]; % Cutoff frequencies
a = [1 0]; % Desired amplitudes
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);
18 y=resample(x,p,q) 数字信号中的重采样。

这时输出信号y的采样频率就
是x的p/q倍，其长度为length(x)*p/q
19 conv 卷积 deconv 反卷积或者求多项式乘法。

xcorr互相关函数 cov协方差 fft2二维FFT fft2二维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);hold on;
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<Ws；
3．带通滤波器：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滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。

在期望通带平滑的情况下，可使用 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=[W1 W2],W1≤ω≤W2。

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)的中心n－2个点等效于 triang(n－2)。

4．汉明窗：w=hamming(n)，产生一个 n点的汉明窗函数。

汉明窗系数为w(k+1)=0.54－0.46cos（） k=0,…,n－1
5．汉宁窗：w=hanning(n)，产生一个 n点的汉宁窗函数。

汉宁窗系数为w(k)=0.5[1－cos( )] k=1,…,n
6．布莱克曼窗：w=Blackman(n),产生一个 n点的布莱克曼窗函数。

布莱克曼窗系数为w(k)=0.42－0.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) 21≤α≤50
0 α<21
增加β可使主瓣变宽，旁瓣的幅度降低。

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(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;
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*(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);
[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)));title('N=50');xlabel('频率
/rad');ylabel('幅度/dB');
subplot(4,1,4);plot(rad,20*log10(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.4
fc1=10;
fc2=20;
fs=100;
w1=2*fc1/fs; w2=2*fc2/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(2*pi*t*15)+sin(2*pi*t*30);%混和正弦波信号
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(10πt)+sin(30πt)+sin(60πt),生成的仿真图和滤波效果如图2所示。

（1）Simulink仿真图（2）滤波前后的离散波形
图 2 Simulink仿真图和滤波效果图
5 SPTool设计法
SPTool是 MATLAB信号处理工具箱中自带的交互式图形用户界面工具，它包含了信号
处理工具箱中的大部分函数，可以方便快捷地完成对信号、滤波器及频谱的分析、设计和浏览。

在本例中按以下步骤完成滤波器的设计和滤波：
创建并导入信号源。

在 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信号序列。

分别在Signals、Filters、Spectra列表中
选择s、filt1、mtlbse，单击 Filters列表下的 Apply按钮，在弹出的 Apply Filter对话框中将输出信号命名为sin15hz。

（4）进行频谱分析。

在 Signals中选择s，单击 Spectra下的 Create按钮，在弹出
的 Spectra Viewer界面中选择 Method为FFT，Nfft=512，单击 Apply按钮生成 s的频谱spect1。

同样的步骤可以生成信号 sin15hz的频谱spect2。

分别选中信号s、sin15hz、spect1、spect2，单击各自列表下方的 View按钮，即可
观察他们的波形，如图 4所示。

图 3 SPTool主界面图4 滤波前后的时域波形和频域特性
由图 4可以看出，带通滤波器 filt1使输入信号 s中频率为 15hz的正弦波信号通过，而将频率为 5hz和 30hz的正弦波信号大大衰减。

matlab语音信号处理
fs=22050; %语音信号采样频率为 22050
x1=wavread('Windows Critical Stop.wav'); %读取语音信号的数据，赋给变量 x1
sound(x1,22050); %播放语音信号
y1=fft(x1,1024); %对信号做 1024点 FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(x1) %做原始语音信号的时域图形
title('原始语音信号');
xlabel('time n');
ylabel('fuzhi n');
figure(2)
freqz(x1) %绘制原始语音信号的频率响应图
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('原始语音信号频谱')
xlabel('Hz');
ylabel('fuzhi');
程序2：
fs=22050; %语音信号采样频率为 22050
x1=wavread('Windows Critical Stop.wav'); %读取语音信号的数据，赋给变量 x1 t=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);
figure(1)
plot(x2) %做原始语音信号的时域图形
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(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;
x1=wavread('h:\课程设计 2\shuzi.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;
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]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换
[H,W]=freqz(bz,az); %绘制频率响应曲线
figure(1)
plot(W*Fs/(2*pi),abs(H))
grid
xlabel('频率／Hz')
ylabel('频率响应幅度')
title('Butterworth')
f1=filter(bz,az,x2);
figure(2)
subplot(2,1,1)
plot(t,x2) %画出滤波前的时域图
title('滤波前的时域波形');
subplot(2,1,2)
plot(t,f1); %画出滤波后的时域图
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('滤波后的频谱')
xlabel('Hz');
ylabel('fuzhi');
程序4:窗函数法设计滤波器：
fs=22050;
x1=wavread('h:\课程设计 2\shuzi.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(1)
%取整
%选择窗函数，并归一化截止频率
freqz(b,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');
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：wavplay
6: guide
matlab 滤波
matlab滤波的方法包括平滑滤波，fir，irr，小波滤波，小波包滤波，自适应lms，rls滤波，最佳 fir滤波，fadtool，或者使用 filter design工具箱。

filter适合于平滑滤波，对于无特殊要求的滤波。

y = filter(b, a, x)，其中b,a为滤波器系数。

计算系统在输入x作用下的零状态响应y[k]。

可以结合filter函数和fir滤波或者
irr滤波。

b,a可以使fir或者irr产生的滤波器系数。

fir滤波线性相位滤波器适用于比较严格的场合，irr对于一般的滤波，对相位无严格
要求。

小波滤波可以实现多频率的，多尺度的信号分解，把信号分为高频和低频，在振动信号处理中，高频信号一般为噪声。

小波包是小波的升级和改进，可以对高频实现更高的分辨率分解。

1 基于 filter平滑滤波（效果没有 fir，iir好），这个滤波器的指标应该是最小
的均方差。

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);
plot(t,x,'r',t,y,'g')。

从图中看出 filter方法可实现平滑滤波，但有相对迟滞。

2 filtfilt 零相位平滑滤波（不推荐，这里 b的长度取值对滤波后信号的幅
值有影响）
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
yy=filtfilt(b,1,x);
plot(t,x);hold on;
plot(t,x,'r',t,yy,'g')
3 IIR滤波器,结合零相位滤波 filtfilt
fp1=3;fs1=10;f2=40;Fs=100;t=0:1/fs:1;Rp=1;As=30;%此处采样率并不是越高越好，一般设为信号最高频率 5-15倍即可。

x=sin(2*pi*fp1*t)+.25*sin(2*pi*f2*t);
[n,wn]=buttord(2*fp1/Fs,2*fs1/Fs,Rp,As);
[b,a]=butter(n,wn);%默认为低通滤波
yf=filter(b,a,x);%用平滑滤波
y=filtfilt(b,a,x) ;
plot(t,x,'r',t,y,'g',t,yf,'b')
4 传统 FIR滤波器
方法一采用窗的精度计算所需滤波器阶数。

fp1=3;fs1=4 ;f2=40;Fs=100;t=0:1/Fs:4;
x=sin(2*pi*t*fp1)+.25*sin(2*pi*t*f2) ;
wp = fp1*pi/50;ws =fs1*pi/50;
deltaw= ws - wp;% 过渡带宽Δω的计算
N = ceil(6.6*pi/ deltaw) + 1; % 按海明窗计算所需的滤波器阶数N0 wdham = (hamming(N+1))';% 海明窗计算
Wn=(3+4)/100;% 计算截止频率
b=fir1(N,Wn,wdham);
y=filter(b,1,x) ;
yt=filtfilt(b,1,x) ;此函数要求 x长度至少为滤波器阶数的 3倍。

方法2 直接设出滤波器阶数，然后凑试
fp1=3;fs1=10;f2=40;Fs=100;t=0:1/Fs:4
x=sin(2*pi*t*fp1)+.25*sin(2*pi*t*f2)
plot(t,x) ;n=50;
b=fir1(n,2*fs1/Fs,hamming(n+1))
y=filter(b,1,x);plot(t,x,t,y)
yz=filtfilt(b,1,x);plot(t,x,t,yz)
5 小波去噪
fp1=3;fs1=10;f2=40;Fs=100;t=0:1/Fs:4
x=sin(2*pi*t*fp1)+.25*sin(2*pi*t*f2)
[c,l]=wavedec(x,3,'db6')；
a3=wrcoef('a',c,l,'db6',3)；
plot(t,a3,'r',t,y,'g')
可以看出效果很好。

6 自适应 lms滤波器 kalman filter rls ,lms,adapltkalman
6.1 自适应 lms滤波器
[y,e,s]=adaptlms(x,d,s)；%x为输入，y为输出，d为期望信号，e为误差信号，s包
含自适应滤波器的结构体。

s=initlms(h0,Mu),ho为滤波器系数的初始值；Mu为步长参数。

6.2 自适应卡尔曼滤波器
Y = adaptkalman(X,D,S)
例子1: 自适应 lms算法：31阶 fir滤波器系统辨识。

（进行 500次迭代）
x = 0.1*randn(1,500); % Input to the filter
b = fir1(31,0.5); % FIR system to be identified
d = filter(b,1,x); % Desired signal
w0 = zeros(1,32); % Initial filter coefficients
mu = 0.8; % LMS step size
S = initlms(w0,mu);
[y,e,S] = adaptlms(x,d,S);
stem([b.',S.coeffs.']);
legend('Actual','Estimated');
title('System Identification of an FIR filter');grid on
例子 2 自适应卡尔曼滤波器
x = 0.1*randn(1,500); % Input to the filter
b = fir1(31,0.5); % FIR system to tbe identified
d = filter(b,1,x); % Desired signal
w0 = zeros(1,32); % Initial filter coefficients
K0 = 0.5*eye(32); % Initial state error correlation matrix
Qm = 2; % Measurement noise covariance
Qp = 0.1*eye(32); % Process noise covariance
S = initkalman(w0,K0,Qm,Qp);
[y,e,S] = adaptkalman(x,d,S);
stem([b.',S.coeffs.']);
legend('Actual','Estimated');
title('System Identification of an FIR filter via Kalman filter');grid on;
小波学习
3 小波变换
是一种信号的时间—尺度(时间—频率)分析方法，它具有多分辨率分析
(Multi-resolutionAnalysis)的特点，而且在时频两域都具有表征信号局部特征的能力，是一种窗口大小固定不变，但其形状可改变，时间窗和频率窗都可以改变的时频局部化分析方法。

即在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，很适合于探测正常信号中夹带的瞬态反常现象并展示其成分，所以被誉为分析信号的显微镜。

在实际的信号处理过程中，尤其是对非平稳信号的处理中，信号在任一时刻附近的频
域特征都很重要。

如柴油机缸盖表面的振动信号就是由撞击或冲击产生的，是一瞬变信号，单从时域或频域上来分析是不够的。

这就促使人们去寻找一种新方法，能将时域和频域结合起来描述观察信号的时频联合特征，构成信号的时频谱。

这就是所谓的时频分析法，亦称为时频局部化方法。

5 小波分析和小波包的区别
小波包分解是在小波基础上发展的，比小波分解更高级，对信号的分解重构更能体
现多分辨率的特征，你应该看看ＭＡＬＬＡＴ算法，看看它们的分解示意图
为了克服小波分解在高频段的频率分辨率较差，而在低频段的时间分辨率较差的缺点，
人们在小波分解的基础上提出了小波包分解。

小波包分解提高了信号的时频分辨率。

是一种更精细的信号分析方法。

由于多分辨率分析只对低频进行分解，对高频部分则保留不动，为了分析振动信号的
高频部分，则用小波包分析，它对低频和高频部分进行再分解.
1 wavemenu，打开图形界面窗口
2 常用小波基函数
常用的
haar,’db1’,’db2’,’bior1.3’,’bior2.4’,’bior2.6’,’coif4’。