matlab的fir高通数字滤波器的设计及分析.doc

基于窗函数的FIR高通数字滤波器设计
摘要
无限长脉冲数字滤波器的设计方法只考虑了幅度特性，没有考虑相位特性，所设的滤波器一般是某种确定的非线性相位特性。有限脉冲响应（FIR）滤波器在保证了幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性。
本课题利用MATLAB软件实现。MATLAB是“矩阵实验室”（MATrix LABoratoy）的缩写，是一种科学计算软件，它使用方便，输入简捷，运算高效，内容丰富，因此利用MATLAB软件，通过一系列较为系统的函数法，根据已知的技术指标，就可以设计出满足要求的滤波器。
关键字:MATLAB；窗函数；FIR带阻数字滤波器；线性相位

目录
1.FIR滤波器简介 3
1.1 FIR的特点 3
2.2线性相位 3
2．主要设计内容 5
3.窗函数 6
3.1常用窗函数 6
3.2窗函数的指标 9
4应用窗函数法设计 FIR 数字滤波器的步骤 10
4.1数字高通滤波器的设计： 10
总结 11
参考文献 12
附录 13

1.FIR滤波器简介
  数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。
1.1 FIR的特点
FIR滤波器的主要优点为：系统总是稳定的，FIR 滤波器的系统函数可以表示为
(2-1)

易知，H(z) 在 Z 平面上有 N-1个零点，z=0 是 N-1 阶极点，因此FIR 系统总是稳定的（极点都在单位圆内）。FIR 滤波器的优点之二：容易实现线性相位。当 FIR 系统的单位冲激响应满足 时，该系统具有线性相位。
（N为奇数） （2-2）
（N为偶数） （2-3）
FIR 滤波器的优点之三：允许设置多通带（或多阻带）滤波器。FIR 滤波器的优点之四：FIR 滤波器可以采用 FFT 方法实现其功能，从而大大提高效率。FIR 滤波器的缺点：由于 FIR 系统只有零点，因此这类系统不像FIR 滤波器不像 IIR 滤波器那样容易取得比较好的通带与阻带衰减特性。要取得较好的衰减特性，一般要求 H(z) 的阶次较高。综合起来看， FIR 滤波器具有IIR 滤波器没有的许多特点，得到了越来越广泛的应用。
FIR滤波器的设计方法主要有三种：a.窗函数设计法；b.频率抽样发；c.最小平法抽样法；这里我主要讨论在MATLAB环境下通过调用信号分析与处理工具箱的几类窗函数来设计滤波器并分析与比较其性能
2.2线性相位
一个单一频率的正弦信号通过一个系统，假设它通过这个系统的时间需要t，则这个信号的输出相位落后原来信号wt的相位。从这边可以看出，一个正弦信号通过一个系统落后的相位等于它的w*t；反过来说，如果一个频率为w的正弦信号通过系统后，它的相位落后delta，则该信号被延迟了delta/w的时间。在实际系统中，一个输入信号可以分解为多个正弦信号的叠加，为了使得输出信号不会产生相位失真，必须要求它所包含的这些正弦信号通过系统的时间是一样的。因此每一个正弦信号的相位分别落后，w1*t，w2*t，w3*t。因此，落后的相位正比于频率w，如果超前，超前相位的大小也是正比于频率w。从系统的频率响应来看，就是要求它的相频特性是一条直线。在FIR滤波器的设计中，为了得到线性相位的性质，通常利用实偶对称序列的相频特性为常数0和实奇对称序列为相频特性为常数90度的特点。因此得到的是对称序列，不是因果序列，是不可实现系统，为了称为物理可实现系统，需要将它向右移动半个周期，这就造成了相移特性随时间的变化，同时也是线性变化。
单位脉冲响应h(n)（为实数）具有偶对称或奇对称性，则FIR数字滤波器具有严格的线性相位特性。

2．主要设计内容
利用窗函数法、频率取样法及优化设计方法设计FIR滤波器，绘制出滤波器的特性图。利用所设计的滤波器对多个频带叠加的正弦信号进行处理，对比滤波前后的信号时域和频域图，验证滤波器的效果。
基本思路：从时域出发设计 h(n)逼近理想 hd(n)。设理想滤波器的单位响应在时域表达为hd(n),则Hd(n) 一般是无限长的，且是非因果的，不能直接作为FIR 滤波器的单位脉冲响应。要想得到一个因果的有限长的滤波器单位抽样响应 h(n)，最直接的方法是先将hd(n)往右平移，再迕行截断，即截取为有限长因果序列：h(n)=hd(n)w(n)，并用合适的窗函数迕行加权作为FIR滤波器的单位脉冲响应。按照线性相位滤波器的要求，线性相位FIR数字低通滤波器的单位抽样响应h(n)必须是偶对称的。对称中心必须等于滤波器的延时常数，即用矩形窗设计的FIR 低通滤波器，所设计滤波器的幅度函数在通带和阻带都呈现出振荡现象，且最大波纹大约为幅度的9%，返个现象称为吉布斯（Gibbs）效应。为了消除吉布斯效应，一般采用其他类型的窗函数。MATLAB 设计 FIR 滤波器有多种方法和对应的函数。窗函数设计法不仅在数字滤波器的设计中占有重要的地位，同时可以用于功率谱的估计，从根本上讲，使用窗函数的目的就是消除由无限序列的截短而引起的Gibbs现象所带来的影响。

3.窗函数
加窗处理使得得滤波器的频率响应与理想滤波器的频率响应之间产生差异，表现为过渡带和波动的出现。我们希望所设计的滤波器尽量逼近理想滤波器，就要设法减少波动的幅度，同时使过渡带变窄。
在设计 FIR 数字滤波器时，窗函数的频谱应该满足：
1主瓣宽度尽可能的窄，以使过渡带尽量陡峭；
2最大旁瓣相对于主瓣尽可能的小，使能量尽可能集中于主瓣内，这样能够使得波动减小。
3.1常用窗函数
1汉宁（Hanning）窗
汉宁窗函数的时域形式表示为：
(3-1)

频域形式为

（3-2）
汉宁窗函数的最大旁瓣值比主瓣值低31dB，但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍，为8π/N。
hanning函数：生成汉宁窗
调用方式：
(1) w = hanning(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。
注意：此函数不返回是零点的窗函数的首尾两个元素。
(2) w = hanning(n,'symmetric')：与上面相类似。
(3) w = hanning(n,'periodic')：此函数返回包括为零点的窗函数的首尾两个元素。

图3-1汉宁窗及其频谱特性
注释：汉宁窗又被称为升余弦窗，汉宁窗可以看做三个矩形时间窗的频谱之和，汉宁窗主瓣加宽并降低，旁瓣显著减小，分辨率下降
2汉明（Hamming）窗:
函数的时域形式可以表示为
（3-3）
频域形式为： （3-4）
其中，为矩形窗函数的幅度频率特性函数。
海明窗函数的最大旁瓣值比主瓣值低41dB，但它和汉宁窗函数的主瓣宽度是一样大的。
Hamming函数：生成海明窗
调用方式
(1) w = hamming(n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。
(2) w = hamming(n,sflag)：参数sflag用来控制窗函数首尾的两个元素值；其取值为symmetric或periodic；默认值为symmetric。

图3-2 汉明窗及其频谱特性
注释：汉明窗和汉宁窗都是余弦窗，只是加权系数不同，汉明窗旁瓣更小
3布莱克曼窗函数的时域形式可以表示为
(3-5)
它的频域特性为
（3-6）
其中，为矩形窗函数的幅度频率特性函数。
布莱克曼窗函数的最大旁瓣值比主瓣值低57dB，但是主瓣宽度是矩形窗函数的主瓣宽度的3倍，为12π/N。
Blackman函数：生成海明窗
调用方式
(1) w = blackman (n)：输入参数n是窗函数的长度；输出参数w是由窗函数的值组成的n阶向量。
(2) w = blackman (n,sflag)：参数sflag用来控制窗函数首尾的两个元素值；其取值为symmetric或periodic；默认值为symmetric。

图3-3 布莱克曼窗及其频谱特性
注释：布莱克曼窗最大旁瓣值比主瓣值滴57db，但主瓣宽度是矩形窗的3倍

3.2窗函数的指标
表3-1如图所示对于任意一个窗函数，求出它的频域值，并求出在主瓣边笫1个零点的位置；2，求出主瓣在-3dB处的位置；3，把笫1个零点的位置-(-3dB)处的位置，就是过渡带的精确带宽。矩形窗最简单，但其 -21dB 的阻带最小衰减在实际应用中远远不够。另外，矩形窗还会造成很强的吉布斯效应。三角窗的阻带衰减性能与矩形窗相比有所改善，但代价是过度带加宽。

4应用窗函数法设计 FIR 数字滤波器的步骤
4.1数字高通滤波器的设计：
令高通滤波器的频率响应为

（4-1）
则

（4-2）
从上述结果可以看出，一个高通滤波器相当于用一个全通滤波器减去一个低通滤波器。

总结
设计带通滤波器时首先要计算出过渡带，然后查表得到不同窗函数所需要的阶数，不同的窗函数所设计的滤波器的形状各有差异，尤其在主瓣宽度、旁瓣的形状以及主瓣与旁瓣的高度差上有比较明显得差别，实际应用中应根据实际情况，折衷处理，兼顾各项指标。
为了这次课程设计，自己自学了数字信号处理领域中窗函数的有关知识。实际中遇到的离散时间信号总是有限长的，因此不可避免地要遇到数据截断问题。而在信号处理中，对离散序列的数据截断是通过序列与窗函数相乘来实现的。而且，有关滤波器的设计、功率谱估计等基本概念也要用到窗函数。本次课程设计对经常用到的下面6窗函数：矩形窗函数、三角窗函数、汉宁窗函数、哈明窗函函数、布莱克曼窗函数、凯塞窗函窗，先是做了基本概念上的阐释，然后对其MATLAB实现函数做出了说明，最后又结合具体的实例，对这些窗函数的频域特性等进行了介绍。
通过这次学习，我不但掌握了FIR数字滤波器窗函数的基本知识及其实际应用的技巧了，还提高了自己的编程和写报告的能力，收获颇多。

参考文献
[1]《数字信号处理》（第三版），丁玉美，高西全.西安电子科技大学出版社，2000.
[2]《MATLAB及在电子信息课程中的应用》，陈怀堔，吴大正，高西全.电子工业出版社，2006.
[3]《MATLAB 7.0从入门到精通》，求是科技.人民邮电出版社，2006.
[4]《数字信号处理（第三版）》学习指导，高西全，丁玉美.西安科技大学出版社，2001.

附录
FIR高通数字滤波器程序：
clear;clc;close all;
[S,Fs,Bit] = wavread('C:\Users\Administrator\Desktop\FIR和IIR_\FIR\SHE.wav'); %读取音频信号
sound(S,Fs);

fnoise1 = 10000; % 加入频率为fnoise的正弦噪声信号
fnoise2 = 50; % 加入频率为fnoise的正弦噪声信号
N = length(S);
T= N/Fs;
t=T/N:T/N:T;
for j=1:N
Noise1(j)=sin(2*pi*fnoise1*t(j));
Noise2(j)=sin(2*pi*fnoise2*t(j));
end
% 对于低通滤波器：
%S1 = S+0.1*Noise1';
%fp = 3000;
%fs = 4000;

%对于高通滤波器：
S1 = S+Noise2';
fp = 3200;
fs = 3000;

% 对于带通滤波器：
% S1 = S+0.1*Noise1'+Noise2';
% fp1 = 500;
% fs1 = 100;
% fp2 = 3000;
% fs2 = 4000;
sound(S1,Fs);

ap=1;as=100;
Fs=44100; %Fs抽样频率
wp=(fp*2*pi)/Fs; %关于π的归一化通带截止频率
ws=(fs*2*pi)/Fs; %关于π的归一化阻带截止频率
DB=wp-ws; %过渡带宽
N0=ceil(6.2*pi/DB); %计算所需h(n)长度N0，ceil(x)取大于等于x的最小整数
N=N0+mod(N0+1,2); %确保h(n)长度N是奇数
wc=(wp+ws)/2; %计算理想高通滤波器通带截止频率（关于π归一化）
h1=fir1(N-1,wc/pi,'high',hamming(N));%调用fir1计算高通FIRDF的h(n)

[H,f]=freqz(h1,1,N,Fs);%求滤波器幅度响应，设置最大幅度为1
figure(1);
plot(f,abs(H));
figure(2);
freqz(h1); %滤波器幅度响应(db)和相位响应
title('hamming滤波器频率响应（幅频上/相频下）');

%h2 = fir1(N-1,wc/pi,hanning(N));

subplot(2,1,1);stem(h1); %滤波器单位冲击响应序列及幅频响应
F=abs(fft(h1));
xlabel = 0:Fs/(N-1):Fs/2;
subplot(2,1,2);plot(xlabel,F(1:(N+1)/2));

X1 = abs(fft(S));
Z1 = (0:150000)/150000;
subplot(3,1,1);plot(Z1,X1(1:150001));

X2 = abs(fft(S1));
Z2 = (0:150000)/150000;
subplot(3,1,2);plot(Z2,X2(1:150001));

S2=conv(S1,h1);
sound(S2,Fs);
X3=abs(fft(S2));
Z3 = (0:150000)/150000;
subplot(3,1,3);plot(Z3,X3(1:150001));

上图为滤波器单位冲击响应序列及幅频响应

注释：FIR滤波器的相位在通带内是线性的。在通频带窗函数都几乎不衰减。hanning窗可以更好的滤除阻带内的信号。

注释：第一张为没有加入噪声前的频谱图，第二张为噪声频谱图，第三章为滤波后频谱图从上面滤波前后语音信号的幅频图可以看出，原始的语音信号大部分500Hz-1000Hz以内。而高通滤波器的通带的频率为3200Hz左右。因此，它可以滤除语音信号低频段的信号。从图三中可以看出，滤波后的语音信号低频段信号幅度几乎为0，其信号大部分集中在高频段。通过回放这段语音信号可以明显听出语音信号滤波前后的区别，滤波器滤除了原始语音信号的低频成分，而保留了高频成分。

注释：此图为高通滤波器的幅度响应

18