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。
Ω=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