数字信号处理实验二

实验二
信号的分析与处理综合实验
一、实验目的
目的：综合运用数字信号处理的理论知识进行信号的采样，重构，频谱分析和滤波器的设计，通过理论推导得出相应结论，再利用Matlab作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。

二、基本要求
1．掌握数字信号处理的基本概念、基本理论和基本方法；
2．学会MA TLAB的使用，掌握MA TLAB的程序设计方法；
3．掌握用MA TLAB设计简单实验验证采样定理的方法；
4．掌握在Windows环境下语音信号采集的方法；
5．学会用MA TLAB对信号进行频谱分析；
6．掌握MA TLAB设计FIR和IIR数字滤波器的方法；
三、实验内容
1．利用简单正弦信号设计实验验证采样定理：
（1）Matlab产生离散信号的方法，作图的方法，以及基本运算操作
（2）对连续正弦信号以不同的采样频率作采样
（3）对采样前后信号进行傅立叶变换，并画频谱图
（4）分析采样前后频谱的有变化，验证采样定理。

掌握画频谱图的方法，深刻理解采样频率，信号频率，采样点数，频率分辨率等概念
2．真实语音信号的采样重构：录制一段自己的语音信号，并对录制的信号进行采样；画出采样前后语音信号的时域波形和频谱图；对降采样后的信号进行插值重构，滤波，恢复原信号。

（1）语音信号的采集
（2）降采样的实现（改变了信号的采样率）
（3）以不同采样率采样后，语音信号的频谱分析
（4）采样前后声音的变化
（5）对降采样后的信号进行插值重构，滤波，恢复原信号
3．带噪声语音信号的频谱分析
（1）设计一频率已知的噪声信号，与实验2中原始语音信号相加，构造带噪声信号
（2）画出原始语音信号和加噪声后信号，以及它们的频谱图
（3）利用频谱图分析噪声信号和原语音信号的不同特性
4．对带噪声语音信号滤波去噪：给定滤波器性能指标，采样窗函数法或双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采样的语音信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；
（1）分析带噪声信号频谱，找出噪声所在的频率段
（2）利用matlab中已有的滤波器滤波
（3）根据语音信号特点，自己设计滤波器滤波
（4）比较各种滤波器性能（至少四种），选择一种合适的滤波器将噪声信号滤除
（5）回放语音信号，比较滤波前后声音的变化
四、实验原理
参考《数字信号处理》教材
《数字信号处理的MATLAB实现》万永革编著
五、主要实验仪器及材料
微型计算机、Matlab。

六、实验步骤
1．设计一简单正弦信号，通过改变采样率观察采样前后的信号变化。

例如：假设有一振幅为1，频率为10Hz，相位为0.3的模拟信号，即sin(2100.3)
π⨯⨯+，用0.01s
t
的采样间隔（采样频率为100Hz）来表示原始信号（注意：实际上模拟信号不能用离散值表示，此处为了在计算机上表示，用采样率非常高的离散信号表示模拟信号）。

分别以5Hz，10Hz（每秒采样10次，即采样间隔为0.1s），20Hz，40Hz，80Hz，200Hz对原始信号进行采样，画出采样前后的信号，并画出其频谱图，对比前后的变化，验证采样定理。

（1）可以用t=0:1/fs:9/f; 取9个周期，通过改变采样率，自动改变采样点数。

（2）也可以通过设置dt1(采样间隔)，已知采样点数n1，t1=n1*dt1，
如图所示，采样率为40Hz时的原始信号，采样过程和采样后的信号时域图和频谱图，可见，当采样率大于原始信号频率的两倍时，采样前后信号频率基本不发生变化，信号不失真。

图2-1 采样的过程
图2-2 原始信号和采样后信号的频谱
2．语音信号的采集。

利用windows下的录音机或其他软件，录制一段自己的话音，时间控制在1秒左右。

然后在MA TLAB 软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。

通过wavread函数的使
用，要求理解采样频率、采样位数等概念。

wavread函数调用格式：
y=wavread(file)，读取file所规定的wav文件，返回采样值放在向量y中。

[y,fs,nbits]=wavread(file)，采样值放在向量y中，fs表示采样频率（Hz），nbits表示采样位数。

y=wavread(file,N)，读取前N点的采样值放在向量y中。

y=wavread(file,[N1,N2])，读取从N1点到N2点的采样值放在向量y中。

3．对真实语音信号的采样、重构。

（1）降采样：
利用windows下的录音机录制的音频采样率是固定的fs（=22050），可以选择以下函数实现对语音信号的降采样。

y=x((1:N:length(x))); %对原始信号每隔N个点取一位，即采样率变为原来的1/N
y=resample(yn,L,M); %采样率变为原来的L/M倍
y=downsample(yn,N); %%采样率变为原来的1/N倍
改变采样率为原来的1/2倍，1/4倍，1/20倍，1/50倍，1/100倍等，分别画出降采样前后的信号波形和频谱图，分析采样前后信号的变化。

图2-3 采样率为原来的1/2时的信号波形频谱图
图2-4 采样率为原来的1/10时信号波形频谱图
在MA TLAB中，函数sound可以对声音进行回放。

其调用格式：sound(x,fs,bits)；
通过调用此函数，感觉采样前后声音的变化。

（2）重构原信号：
降采样后，信号的采样率和采样点数同时变化。

如采样率变为原来的1/2，即对原始信号每隔一个点采样。

如果要恢复原始信号，即信号长度和采样率须变为原来同样大小。

所以，必须对降采样后信号进行插值重构。

具体过程参见附件1中例子。

对采样后的真实语音信号进行插值重构，滤波，恢复原始信号。

画出插值前后信号的波形以及频谱图，并将重构后信号与原始信号进行比较。

如，对采样率降为原来1/5的降采样后信号插值重构，结果如下图所示。

图2-5 采样率为原来的1/5时的波形频谱图
图2-6 插值后的信号波形频谱图
图2-7 低通滤波器的频率响应图
图2-8 滤波后的波形频谱图
调用sound函数感受插值后的声音，发现会有高频的噪声。

经过低通滤波器之后，高频噪声被滤出。

但是，因为之前的原始信号经过了降采样，所以插值后的效果一定不如原始声音。

也可用用wavwrite函数将经过处理的语音信号保存下来，调用格式为如wavwrite(y,fs,bits,'sound.wav')，其中，y为所要保存的语音信号，fs为其采样率，bits为采样位数，'sound.wav'是保存的语音信号的文件名。

4．语音信号的频谱分析
要求：首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在MA TLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性；从而加深对频谱特性的理解。

[x,fs,bits] =wavread('e:\sound1.wav',[1024 5120]); %读取[1024 5120]段%语音信号数据；如果不指定所读取读取语音信号数据的长度，则读取整段语音信号
X=fft(x，4096); % fftshift(fft(x))
N=length(x);
k=0:N-1;
subplot(221);plot(k,x);title('原始信号波形');
subplot(222);plot((-N/2:N/2-1)/N*fs,abs(X));title('原始信号频谱');
5．对原始语音信号加噪声：
设计频率已知的噪声信号或者用自然噪声信号加在原始语音信号上，构建带噪声信号。

如对原始信号加上频率为5000的正弦波噪声信号，程序如下：
[x,fs,bits] =wavread('sound1.wav'); %读取原始语音信号
X=abs(fftshift(fft(x))); %原始信号频谱
c1=0.01*sin (2*pi*5000*k/fs); %构建频率为5000Hz的正弦波噪声信号
yn=x+c1'; %构建带噪声信号
Yn=abs(fftshift(fft (yn))); %求带噪声信号频谱
subplot(321);plot(x);title('原始信号波形');
subplot(322);plot((-N/2:N/2-1)/length(k)*fs,X);title('原始信号频谱');
subplot(323);plot(yn);title('带噪声信号波形');
subplot(324);plot((-N/2:N/2-1)/length(k)*fs,Yn);title('带噪声信号频谱');
运行程序，结果如图2-9所示。

从频谱图可以看出，原始信号频谱段集中在0~5000Hz之间，主要频率集中在0Hz附近的低频部分。

加上噪声信号后，在5000Hz处有个幅值非常大的的高频成份，即以上所加的高频正弦波噪声信号。

图2-9 加噪声前后信号的波形频谱图
6．设计数字滤波器和画出频率响应
根据分析所得的原始信号的频谱和噪声信号频谱特点，给出有关滤波器的性能指标。

首先用窗函数法或者最优化法设计高通，低通，带通，带阻滤波器，在MA TLAB中，可以利用函数fir1，firls设计FIR滤波器；
然后在用双线性变换法或脉冲响应法设计上面几种滤波器，在MA TLAB中，可以利用函数butte、cheby1和ellip设计IIR滤波器；
最后，利用MA TLAB中的函数freqz画出各滤波器的频率响应。

具体方法参加附件3种滤波器设计的步骤和实例。

如，根据以上语音信号的特点给出有关IIR滤波器的性能指标：
1）低通滤波器性能指标，fp=4500，fc=6500，Rs=100，Rp=1。

（fp：通带截至频率；fc：阻带截至频率；Rs：通带波纹；Rp：阻带波纹）
2）带阻滤波器性能指标，fp1=4800 Hz，fp2=5200 Hz，fc1=4600 Hz，fc2=5400 Hz，Rs=30dB，Rp=1dB。

（[fp1 fp2]：阻带截至频率；[fc1 fc2]：通带截至频率）
程序如下：
%椭圆带阻滤波器
figure,
wp=[4800 5200]*2/fs;ws=[4600 5400]*2/fs; %通带和阻带边界频率
Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数
[NN,Wn]=ellipord(wp,ws,Rp,Rs); %求取数字滤波器的最小阶数和归一化截至频率
[b,a]=ellip(NN,Rp,Rs,Wn,'stop'); %设计滤波器
freqz(b,a,512,fs); title('椭圆带阻滤波器');
图 2-10 椭圆带阻滤波器的频率响应
%双线性变换法设计的椭圆低通滤波器：
figure,
fp=4500;fc=6500;Rs=100;Rp=1;fs=22050;
wc=2*fc/fs;wp=2*fp/fs;
[n,wn]=ellipord(wp,wc,Ap,As);
[b,a]=ellip(n,Rp,Rs,wn);
freqz(b,a,512,fs);
若设计FIR滤波器，要阻止5000Hz的高频信号通过，即设计如下滤波器特性，% 最优化设计法设计带阻滤波器
figure,
n=100; %滤波器阶数
f=[ 0 4600*2/fs 4800*2/fs 5400*2/fs 5600*2/fs 1 ]; %频率向量
a=[ 1 1 0 0 1 1 ]; %振幅向量
b=firls(n,f,a); %采用 firls 设计滤波器
[h,w1]=freqz(b); %计算其频率响应
plot(w1/pi,abs(h)); %绘制滤波器的幅频响应
title('最优化法设计的带阻滤波器');
图 2-11 最优化法设计的带阻滤波器的幅频响应
7．用滤波器对信号进行滤波
比较各种滤波器的性能，然后用性能好的各滤波器分别对采集的信号进行滤波。

比较滤波前后语音信号的波形及频谱，要求在一个窗口同时画出滤波前后的波形及频谱。

在MA TLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。

xx=filter(b,a,y); %对信号y滤波
XX=fftshift(fft(xx)); %求经过滤波后信号的傅立叶变换
图2-12 对带噪声信号滤波
8、图形用户界面的设计
为使编制的程序操作方便，有兴趣的同学，可以利用Matlab进行图形用户界面的设计。

在所设计的系统界面上可以选择滤波器的类型，输入滤波器的参数，显示滤波器的频率响应，选择信号等。

七、实验报告要求
1.按照实验指导书的格式简述实验目的，基本要求，实验内容以及实验原理。

2.详细写明每一个实验的具体步骤，如何实现，在具体编程中遇到的问题以及如何解决。

3.按照实验要求编程，给出详细程序实现代码，以及程序结果（图，表格等）。

4.在每一步中，分析各种不同的情况所得出的结果，并进行对比，得出结论。

5.在本次实验报告的最后，总结本次实验的主要内容，以及所掌握的内容。

附件1：信号的插值重构
以对一个简单正弦信号的重构为例，以10倍原信号长度进行插值。

clear;
close all;
f=10;
fs=100;
k=1:100;
x=sin(2*pi*f*k/fs);
n=10;
y=zeros(1,n*length(x));
y(1:n:length(y))=x;
附图2-1 信号插值前后的波形
附图2-2 信号插值前后的频谱图
由上图可看出，插值后信号最大频率同样变为原来的10倍，高频部分进行了拓展（由于在原始信号中插入0的缘故），为原始频谱的周期延拓（类似采样定理）。

故，只需根据以上插值后信号的频率特性，
设计FIR低通滤波器进行滤波，即可得到（滤波器具体的设计方法见附录）。

低频滤波器的通带为0~50Hz，由于1对应最大频率500Hz，故，50Hz对应的转换频率为50/500=0.1，设计通带截止频率为wp为0.09，阻带截止频率为0.11。

N为滤波器阶数，理想状况下，阶数越高，滤波器效果越好。

（工程中涉及到的滤波器阶数问题需另做考虑）
N=100;
B=firls(N,[0 0.09 0.11 1],[1 1 0 0]);
figure,plot((-N/2:N/2)/(N/2),abs(fftshift(fft(B))))，title(‘FIR低通滤波器’);
xx=filter(B,1,y); % 利用低通滤波器滤波
figure,
subplot(211),plot(xx);title('滤波后信号');
subplot(212),plot((-n*50:n*50-1)/(n*100)*n*fs,abs(fftshift(fft(xx))));title('滤
波后信号频谱');
附图2-3 FIR低通滤波器的频率响应
附图2-4 插值后的信号经低通滤波的信号波形和频谱图
附件2：信号频谱
有限长序列的离散傅立叶变换（DFT ）定义为
1
()[()](),01N kn N
n X k DFT x n x n W
k N --===
≤≤-∑
逆变换为：
1
1()[()](),01N kn N
k x n ID FT X k X k W
k N N
--===
≤≤-∑
Matlab 频谱分析相关函数：
（1）一维快速傅立叶变换FFT 函数fft ，调用格式：
y=fft(x) y=fft(x,n)
函数说明：fft 函数用于计算矢量或矩阵的离散傅立叶变换，可通过
1
(1)(1)N kn
N
n X k x n W -=+=
+∑
来实现，式中(2/)
(),j N N N len g th x W e
π-==。

Y=fft(x)利用FFT 算法计算矢量x 的离散傅立叶变换，当x 为矩阵时，y 为矩阵x 的每一列的FFT 。

当x 的长度为2的幂次方时，则fft 采用基-2FFT 算法，否则采用稍慢的混合基算法。

Y=fft(x,n)采用n 点FFT 。

当x 长度小于n 时，fft 函数在x 尾部补零.。

当x 长度大于n 时，函数将序列截断.
（2）一维逆快速傅立叶变换（IFFT ）ifft ，调用格式：
Y=ifft(x) Y=ifft(x,n)
函数说明：ifft 函数用于计算矢量或矩阵的逆傅立叶变换，即
1
1(1)(1)N kn
N
k x n X k W N
--=+=
+∑
式中，(2/)
(),j N N N length x W e
π-==。

（3）fft 函数最通常的应用是计算信号的频谱。

考虑由一个50Hz 和120Hz 正弦信号构成的信号，受零均值随机信号的干扰，数据采样率为1000Hz 。

通过fft 函数来分析其信号频率成分。

相关程序如下：
t=0:0.001:0.6;
x=sin(2*pi*50*t)+sin(2*pi*120*t); y=x+1.5*randn(1,length(t)); Y=fft(y,512);
附件3：滤波器设计
一、IIR 数字滤波器设计
1、经典设计法：主要有两种方法：脉冲响应不变法、双线性变换法
经典设计IIR 滤波器的步骤：
A ．根据给定的性能指标和方法，首先对设计性能指标中的频率指标进行转换，转换后的频率指标作为模拟滤波器原型设计指标；
B ．估计模拟滤波器最小阶数和边界频率，可利用Matlab 工具函数buttord ，cheb1ord 等；如，设计buttord 滤波器的调用格式：
[n,Wn]=buttord(Wp,Ws,Rp,Rs) [n,Wn]=buttord(Wp,Ws,Rp,Rs,’s ’)
buttord 可在给定滤波器性能的情况下，选择模拟或数字滤波器的最小阶数，其中，Wp 和Ws 分别是通带和阻带的截止频率，取值范围为0<Wp(或Ws)<1,其值为1时表示0。

fs 为采样频率，Rp 和Rs 分别是通带和阻带区的波纹系数。

C ．设计模拟低通滤波器原型，可利用buttap ，cheblap 等；
D ．由模拟低通原型经频率变换得到模拟滤波器（低通，高通，带通，带阻），可利用MA TLAB 工具函数lp2lp ，lp2hp ，lp2bs ，lp2bp 等；如低通到低通模拟滤波器变换函数为[bt,at]=lp2lp(b,a,Wo)，lp2lp 函数将截止频率为1rad/s （归一化频率）的模拟低通滤波器原型变换为截止频率为Wo 的低通滤波器。

E ．将模拟滤波器离散化得到IIR 数字滤波器，可利用MA TLAB 工具函数bilinear ，impinvar 等。

例1．用脉冲响应不变法设计Butterworth 低通滤波器，要求通带频率为00.2ωπ≤≤，通带波纹小于1dB ，阻带在0.3πωπ≤≤
内，幅度衰减大于15dB ，采样周期T=0.01s 。

假设一个信号
12()sin 20.5cos 2X t f t f t ππ=+，其中125,30f H z f H z ==。

试将原信号与经过该滤波器的输出信号进
行比较。

注意，该题给出的通带边界频率和阻带边界频率均为数字频率，因此设计时首先要将其转换为模拟频
率。

%Samp1
wp=0.2*pi;ws=0.3*pi;Rp=1;Rs=15; %数字滤波器截止频率、通带波纹和阻带衰减 T=0.01;Nn=128; %采样间隔
Wp=wp/T;Ws=ws/T; %得到模拟滤波器的频率—采用脉冲响应不变法的频率转换形式 [N,Wn]=buttord(Wp,Ws,Rp,Rs,'s'); %计算模拟滤波器的最小阶数 [z,p,k]=buttap(N); %设计低通原型数字滤波器
[Bap,Aap]=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式 [b,a]=lp2lp(Bap,Aap,Wn); %低通滤波器频率转换
[bz,az]=impinvar(b,a,1/T); %脉冲响应不变法设计数字滤波器传递函数 figure(1)
[H,f]=freqz(bz,az,Nn,1/T); %输出幅频响应和相频响应 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; figure(2)
f1=5;f2=30; %输入信号含有的频率
N=100; %数据点数
n=0:N-1;t=n*T; %时间序列
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号 subplot(2,1,1),plot(t,x),title('输入信号') y=filtfilt(bz,az,x); %对信号进行滤波
subplot(2,1,2),plot(t,y),title('输出信号'),xlabel('时间/s')
程序运行结果如下图所示。

该例所要求的通带频率为00.2ωπ≤≤，而该滤波器的采样间隔为0.01s ，采样频率为100Hz ，即2π对应于100Hz ，则0.2π对应于10Hz ，即该滤波器的通带范围为0~10Hz 。

观看附图2-5，其通带范围最大衰减小于1dB ，与该分析一致。

题意要求阻带在0.3πωπ≤≤内，幅度衰减大于15dB 。

其中，0.3π对应于15Hz ，幅度衰减大于15dB 。

完全符合滤波器设计的要求。

附图2-5 例1所设计滤波器的幅频响应和相频响应
附图2-6 例1所设计滤波器的输入输出信号
例2．用双线性变换法设计一个椭圆低通滤波器，其性能指标同例1。

%Samp2
wp=0.2*pi;ws=0.3*pi;Rp=1;Rs=15; %数字滤波器截止频率通带波纹和阻带衰减
Fs=100;Ts=1/Fs;Nn=128; %采样频率
Wp=2/Ts*tan(wp/2.);Ws=2/Ts*tan(ws/2.); %按频率转换公式进行转换——双线性变换法
[N,Wn]=ellipord(Wp,Ws,Rp,Rs,'s'); %计算模拟滤波器的最小阶数
[z,p,k]=ellipap(N,Rp,Rs); %设计模拟原型滤波器
[Bap,Aap]=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式
[b,a]=lp2lp(Bap,Aap,Wn); %低通转换为低通滤波器的频率转换
[bz,az]=bilinear(b,a,Fs); %运用双线性变换法得到数字滤波器传递函数
[H,f]=freqz(bz,az,Nn,Fs); %求出频率特性
subplot(2,1,1),plot(f,20*log10(abs(H)));
xlabel('频率/Hz');ylabel('振幅/dB');grid on;
subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))
xlabel('频率/Hz');ylabel('相位/^o');grid on;
程序运行如下图，在10Hz以前，衰减小于1dB，在15Hz以后，衰减大于15dB，即性能指标完全满足滤波器设计的要求。

附图2-7 例2所设计滤波器的幅频响应和相频响应
2、IIR滤波器完全设计函数：
以上介绍了IIR滤波器的基本方法和步骤，并给出一些例子说明如何利用MA TLAB编程分步实现这些步骤。

由这些步骤可知，我们必须多次调用MA TLAB信号处理工具箱中的基本工具函数。

实际上，MA TLAB信号处理工具箱还提供了IIR滤波器设计的完全工具函数，用户只需要调用这些工具函数即可一次性地完成设计，而不需要调用那些基本工具函数分步实现。

数字IIR滤波器的完全设计函数有
[b,a]=butter（n, Wn, ’ftype’)
[b,a]=cheby1（n, Rp, Wn, ’ftype’)
[b,a]= cheby2（n, Rs, Wn, ’ftype’) [b,a]=ellip （n, Rp, Rs, Wn, ’ftype’)
在上面的调用函数中，n 代表滤波器阶数，Wn 代表滤波器的截止频率，取值为0~1，需根据采样频率Fs 来定，如滤波器的截止频率为Fc （Hz ），则
Wn=2*Fc/Fs
这样就转换为0~1的归一化频率。

其中，p ω，s ω等边界频率都要根据此公式进行转换。

’ftype’滤波器的类型为：
可取为’high ’高通，截止频率为Wn ；’stop ’带阻，截止频率为[w1,w2](w1>w2)；’ftype’缺省时为低通或带通滤波器。

在调用完全设计函数之前，需要用到阶数估计函数来估计滤波器的阶数，如：
[n,Wn]=butterd （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]。

其他完全设计函数：Cheby1,cheby2,ellip,besself 阶数估计函数：Cheb1ord,cheb2ord,ellipord
例3．设计一个Buttordworth 高通数字滤波器，通带边界频率为300Hz ，阻带边界频率为200Hz ，通带波纹小于1dB ，阻带衰减大于20dB ，采样频率为1000Hz 。

假设一个信号12()sin 20.5cos 2x t f t f t ππ=+，其中1210,400f H z f H z ==。

试将原信号与通过该滤波器的输出信号进行比较。

%Samp3
Fs=1000; %采样频率
wp=300*2/Fs;ws=200*2/Fs; %根据采样频率将滤波器边界频率进行转换 Rp=1;Rs=20; %通带波纹和阻带衰减
Nn=128; %显示滤波器频率特性的数据长度
[N,Wn]=buttord(wp,ws,Rp,Rs); %求得数字滤波器的最小阶数和截止频率（归一化频率） [b,a]=butter(N,Wn,'high'); %设计Butterworth 高通数字滤波器 figure(1)
[H,f]=freqz(b,a,Nn,Fs); %用Nn 点绘出频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))
xlabel('频率/Hz');ylabel('相位/^o');grid on;
n=0:127;
dt=1/Fs;t=n*dt; %时间序列
f1=10;f2=400; %输入信号频率
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号
figure(2)
subplot(2,1,1),plot(t,x); title('输入信号') %绘制输入信号
y=filter(b,a,x); %对输入信号进行滤波
subplot(2,1,2),plot(t,y),title('输出信号') %绘制输出信号
xlabel('时间/s')
程序输出结果如下图，可看出所设计滤波器在大于300Hz为通带，其衰减均小于1dB；在小于200Hz 为阻带，其衰减大于20dB，完全符合设计要求，但相频特性是非线性的。

当滤波器输入10Hz和400Hz 两种信号后，滤波器滤除了10Hz的信号，使得400Hz的信号通过滤波器，起到滤波的效果。

附图2-8 例3所设计滤波器的幅频响应和相频响应
附图 2-9 例3所设计滤波器的输入和输出信号
二、FIR 数字滤波器
如果所希望的滤波器的理想频率响应函数为()j d H e ω
，则其对应的单位脉冲响应为
1()()2j j n
d d h n H e
e
d π
ω
ωπ
ωπ
-=
⎰
窗函数设计法的基本原理是用有限长单位脉冲响应序列()h n 逼近()d h n 。

由于()d h n 往往是无限长序列，且是非因果的，所以用窗函数()n ω将()d h n 截断，并进行加权处理，得到：
()()()d h n h n n ω=
()h n 就作为实际设计的FIR 数字滤波器的单位脉冲响应序列，其频率响应函数()j H e ω
1
()()N j j n
n H e
h n e
ω
ω--==
∑
式中，N 为所选窗函数()n ω的长度。

这样选定窗函数类型和长度 N 后，求出单位脉冲响应()()()d h n h n n ω=，并求出()j H e ω。

()j H e ω
是
否满足要求，要进行验算。

一般在()h n 尾部加零使长度满足2的整数次幂，以便使用FFT 计算()j H e ω。

如果要观察细节，补零点数增多即可。

如果()j H e ω不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。

如果要求线性相位，则()h n 还必须满足；
=±--
h n h N n
()(1)
FIR滤波器设计的主要方法有窗函数法，最优化设计法，约束最小二乘逼近，升余弦函数法等，在这里主要介绍窗函数法和最优化设计方法。

1．窗函数法：
基于窗函数的FIR滤波器设计函数，常用的窗函数有矩形窗、汉宁窗、哈明窗等。

在Matlab中，将上述方法复合成一个函数，即fir1。

调用格式：
b=firl(n,Wn)
b=firl(n,Wn,’ftype’)
b=firl(n,Wn,Window)
b=firl(n,Wn,’ftype’,Window)
函数说明：以经典方法实现加窗线性相位FIR数字滤波器设计，也可设计出标准的低通，带通，高通和带阻滤波器。

例4．用窗函数设计一个线性相位FIR低通滤波器，并满足性能指标：通带边界的归一化频率wp=0.5，阻带边界的归一化频率ws=0.66，阻带衰减不小于30dB，通带波纹不大于3dB。

假设一个信号，其中f1=3Hz，f2=20Hz。

信号的采样频率为50Hz。

试将原信号与通过滤波器的信号进行比较。

由题意，阻带衰减不小于30dB，根据FIR滤波器各种窗函数的基本参数，选择Hanning窗，因为汉宁窗的第一旁瓣相对主瓣衰减为31dB，满足滤波要求。

在窗函数设计法中，要求设计的频率归一化到0~π之间，Nyquist频率对应于π，因此通带和阻带边界频率为0.5π和0.66π。

程序如下：
%Samp4
wp=0.5*pi;ws=0.66*pi; %滤波器的边界频率
wdelta=ws-wp; %过渡带宽度
N=ceil(8* pi/wdelta);%求解滤波器的最小阶数，根据Hanning 窗主瓣宽
Wn=(0.5+0.66)*pi/2;%截止频率取通带和阻带边界频率的中点
b=fir1(N,Wn/pi,hanning(N+1));%设计FIR滤波器,注意fir1要求输入归一化频率
[H,f]=freqz(b,1,512,50); %采用50Hz的采样频率求出频率响应
subplot(2,1,1),plot(f,20*log10(abs(H)))
xlabel('频率/Hz');ylabel('振幅/dB');grid on;
subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))
xlabel('频率/Hz');ylabel('相位/^o');grid on;
f1=3;f2=20; %检测输入信号含有两种频率成分
dt=0.02; t=0:dt:3; %采样间隔和检测信号的时间序列
x=sin(2*pi*f1* t)+cos(2* pi*f2* t); %检测信号
%y=filter(b,1,x); %可采用此函数给出滤波器的输出
y=fftfilt(b,x); %给出滤波器的输出
figure(2)
subplot(2,1,1), plot(t,x),title('输入信号') %绘输入信号
subplot(2,1,2),plot(t,y) % 绘输出信号
hold on; plot([1 1]*(N-1)/2*dt,ylim, 'r') %绘出延迟到的时刻
xlabel('时间/s'),title('输出信号')
程序运行结果如下图所示。

该例设计通带边界为wp=0.5，阻带边界ws=0.66，对应于50Hz的采样频率通带边界频率为fp=Fs/2*Fnormal=50/2*0.5=12.5Hz，fs=50/2*0.66Hz，其中Fs为采样频率，Fnormal为归一化频率。

由附图2-10可以看出，在小于12.5Hz的频段上，几乎看不到下降，即满足通带波纹不大于3dB 的要求。

在大于16.5Hz的频带上，阻带衰减大于30dB，满足设计要求。

由附图2-11可见，在通带范围内，相位频率为一条直线，表明该滤波器为线性相位。

滤波器输入信号包括3Hz和20Hz的信号，20Hz的信号不能通过滤波器，通过滤波器后只剩下3Hz 的信号。

输出结果证明了这一点。

但要注意，由于FIR滤波器所需的阶数较高，信号延迟（N-1）/2也较大，输出信号前面有一段直线就是延迟造成的。

附图2-10 例4所设计滤波器的幅频响应和相频响应
附图2-11 例4所设计滤波器的输入输出信号
2．最优化设计方法：
Matlab提供了比fir1更为通用的函数firls，其调用函数为：
b = firls(n,f,a)
式中，n为滤波器阶数；f为滤波器期望频率特性归一化频率向量，范围为0~1，为递增向量，允许定义重复频率点；a为滤波器期望频率特性的幅值向量，向量a和f必须为同长度，且为偶数；b为返回的滤波器系数，长度为n+1，且具有偶对称的关系，即b(k)=b(n+2-k)。

若滤波器的阶数为奇数，则在Nyquist频率处（对应于归一化频率1），幅频响应必须为0。

滤波器的阶数为偶数，则无此限制。

例5．用firls设计一个20阶FIR低通滤波器，通带边界频率为0.4，幅值为1，阻带边界频率为0.5，幅值为0。

输入一个采样频率为50Hz，频率为5Hz和15Hz的合成振动，将原信号与通过滤波后的信号进行比较。

%Samp5
clf;n=20; %滤波器阶数
f=[0 0.4 0.5 1]; %频率向量
a=[1 1 0 0]; %振幅向量
b=firls(n,f,a); %采用firls 设计滤波器
[h,w1]=freqz(b); %计算其频率响应
figure(1)
plot(w1/pi,abs(h)); %绘制滤波器的幅频响应
xlabel('归一化频率');ylabel('振幅');
legend('firls','remez'); %给出图例
grid on;
figure(2)
fs=50;t=0:1/fs:2; %采样频率和时间序列
f1=5;f2=15; %输入信号频率
x1=sin(2*pi*f1*t)+8.*cos(2*pi*f2*t);
subplot(2,1,1),plot(t,x1)
title('原始信号')
y1=filter(b,1,x1);
subplot(2,1,2),plot(t,y1);
legend('firls','remez'); %给出图例
title('输出信号')
xlabel('时间/s')
程序运行结果如下所示，输入5Hz和15Hz的两个频率振动信号，5Hz的信号对应的采样频率为50Hz
的归一化频率为0.2 (即：
5
50/2
),15Hz的信号对应的归一化频率为0.6（即：
15
50/2
），对应于该滤波器的
幅频响应，归一化频率为0.2的振动能通过滤波器，而归一化为0.6的振动不能通过该滤波器，这一点从由图例2-13也可以看出。