MATLAB窗函数法实现FIR的高通-带通和低通滤波器的程序.doc

MATLAB课程设计报告

学 院：地球物理与石油资源学院
班 级： 测井（基）11001
姓 名： 大牛啊啊啊
学 号：
班内编号：
指导教师： 陈义群
完成日期： 2013年6月3日

一、 题目 FIR滤波器的窗函数设计法及性能比较
1. FIR滤波器简介
数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。与IIR滤波器相比，FIR滤波器的主要特点为：
a. 线性相位；b.非递归运算。
2. FIR滤波器的设计
FIR滤波器的设计方法主要有三种：a.窗函数设计法；b.频率抽样发；c.最小平法抽样法；
这里我主要讨论在MATLAB环境下通过调用信号分析与处理工具箱的几类窗函数来设计滤波器并分析与比较其性能。窗函数法设计FIR滤波器的一般步骤如下：
a. 根据实际问题确定要设计的滤波器类型；
b. 根据给定的技术指标，确定期望滤波器的理想频率特性；
c. 求期望滤波器的单位脉冲响应；
d. 求数字滤波器的单位脉冲响应；
e. 应用。
常用的窗函数有

4. 常用窗函数的参数

5. FIR滤波器的MATLAB实现方式
在MATLAB信号分析与处理工具箱中提供了大量FIR窗函数的设计函数，本次用到主要有以下几种：
hanning(N) hanning窗函数的调用
hamming(N) hamming窗函数的调用
blackman(N) blackman窗函数的调用
kaiser(n+1,beta) kaiser窗函数的调用
kaiserord 计算kaiser窗函数的相关参数
freqz 求取频率响应
filter 对信号进行滤波的函数
6. 实验具体步骤
本次实验分别通过调用hanning ,hamming ,Blackman,kaiser窗函数，给以相同的技术参数，来设计低通，带通，高通滤波器，用上述窗函数的选择标准来比较各种窗函数的优劣，并给以一个简谐波进行滤波处理，比较滤波前后的效果。达到综合比较的效果。
二、源代码
1.利用hanning hamming blackman kaiser窗,设计一个低通FIR
function lowpassfilter
clc;
clear all;

Fs=100;%采样频率
fp=20;%通带截止频率
fs=30;%阻带起始频率
wp=2*pi*fp/Fs;%将模拟通带截止频率转换为数字滤波器频率
ws=2*pi*fs/Fs;%将模拟阻带起始频率转换为数字滤波器频率
wn=(wp+ws)/2/pi;%标准化的截止频率响应
Bt=ws-wp;
N0=ceil(6.2*pi/Bt);%滤波器长度
N=N0+mod(N0+1,2);
window1=hanning(N);%使用hanning窗函数
window2=hamming(N);%使用hamming窗函数
window3=blackman(N);%使用blackman窗函数
[n,Wn,beta,ftype]=kaiserord([20 25],[1 0],[0.01 0.01],100);
window4=kaiser(n+1,beta);%使用kaiser窗函数

%设计加窗函数fir1
b1=fir1(N-1,wn,window1);
b2=fir1(N-1,wn,window2);
b3=fir1(N-1,wn,window3);
b4=fir1(n,Wn/pi,window4 ,'noscale');
%求取频率响应
[H1,W1]=freqzb1,1,512,2);
[H2,W2]=freqz(b2,1,512,2);
[H3,W3]=freqz(b3,1,512,2);
[H4,W4]=freqz(b4,1,512,2);

figure(1);
subplot(2,2,1),plot(W1,20*log10(abs(H1)));%绘制频率响应图形
axis([0,1,-100,100]);
title('低通hanning窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

subplot(2,2,2),plot(W2,20*log10(abs(H2)));%绘制频率响应图形
axis([0,1,-100,100]);
title('低通hamming窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

subplot(2,2,3),plot(W3,20*log10(abs(H3)));%绘制频率响应图形
axis([0,1,-100,100]);
title('低通blackman窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

subplot(2,2,4),plot(W4,20*log10(abs(H4)));%绘制频率响应图形
axis([0,1,-100,100]);
title('低通kaiser窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

T=1/Fs;
L=100;%信号长度
t=(0:L-1)*T;%定义时间范围和步长
y=sin(2*pi*5*t)+5*sin(2*pi*15*t)+8*sin(2*pi*40*t);%滤波前的图形
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;%将时域信号变换到频域
f = Fs/2*linspace(0,1,NFFT/2+1);%频域采样
figure(2);
plot(f,2*abs(Y(1:NFFT/2+1)));xlabel('frequency/Hz');ylabel('Amuplitude') ;%滤波前频谱
title('滤波前的频谱');
%滤波后频谱
%采用hanning窗滤波器
yy1=filter(b1,1,y);%调用滤波函数
YY1=fft(yy1,NFFT)/L;%进行傅里叶变换，下同。
f1=Fs/2*linspace(0,1,NFFT/2+1);
figure(3);
subplot(2,2,1),plot(f1,2*abs(YY1(1:NFFT/21))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('hanning窗的滤波效果');
%采用hammning窗滤波器
yy2=filter(b2,1,y);
YY2=fft(yy2,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,2),plot(f1,2*abs(YY2(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('hamming窗的滤波效果');
%采用blackman窗滤波器
yy3=filter(b3,1,y);
YY3=fft(yy3,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,3), plot(f1,2*abs(YY3(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('blackman窗的滤波效果');
%采用kaiser窗滤波器
yy4=filter(b4,1,y);
YY4=fft(yy4,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,4),plot(f1,2*abs(YY4(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
xlabel('frequency/Hz');ylabel('Amuplitude');
title('kaiser窗函数滤波效果');

%滤波前后的信号的时域对比
figure(4);
plot(y);xlabel('时间/s');ylabel('振幅');title('滤波前振幅特性');
figure(5);
subplot(2,2,1),plot(yy1);xlabel('时间/s');ylabel('振幅');title('hanning窗函数滤波振幅特性');
subplot(2,2,2),plot(yy2);xlabel('时间/s');ylabel('振幅');title('hamming窗函数滤波振幅特性');
subplot(2,2,3),plot(yy3);xlabel('时间/s');ylabel('振幅');title('blackman窗函数滤波振幅特性');
subplot(2,2,4),plot(yy4);xlabel('时间/s');ylabel('振幅');title('kaiser窗函数滤波振幅特性');
%滤波前后的信号的相位对比
figure(6);
plot(angle(Y));xlabel('时间/s');ylabel('相位');title('滤波前的相位特性');
figure(7);
subplot(2,2,1),plot(angle(YY1));xlabel('时间/s');ylabel('相位');title('hanning窗函数滤波相位特性');
subplot(2,2,2),plot(angle(YY2));xlabel('时间/s');ylabel('相位');title('hamming窗函数滤波相位特性');
subplot(2,2,3),plot(angle(YY3));xlabel('时间/s');ylabel('相位');title('blackman窗函数滤波相位特性');
subplot(2,2,4),plot(angle(YY4));xlabel('时间/s');ylabel('相位');title('kaiser窗函数滤波相位特性');
2.设计一个hanning hamming blackman kaiser窗函数bandpass_FIR
%设计一个hanning hamming blackman kaiser窗函数bandpass_FIR
function bandpassfilter
Fs=100;%采样频率
fp1=15;%通带下限截止频率
fp2=20;%通带上限截止频率
fs1=10;
fs2=25;
wp1=2*pi*fp1/Fs;%将通带下限截止频率转换为数字滤波器频率
wp2=2*pi*fp2/Fs;%将通带上限截止频率转换为数字滤波器频率
ws1=2*pi*fs1/Fs;%将通带下限截止频率转换为数字滤波器频率
ws2=2*pi*fs2/Fs;%将通带上限截止频率转换为数字滤波器频率

Bt=wp1-ws1;
N0=ceil(6.2*pi/Bt);
N=N0+mod(N0+1,2);
wn=[(wp1+ws1)/2/pi,(wp2+ws2)/2/pi];
window1=hanning(N);%使用hanning窗函数
window2=hamming(N);%使用hamming窗函数
window3=blackman(N);%使用blackman窗函数
%设过渡带宽度为5Hz
[n,Wn,beta,ftype]=kaiserord([10 15 20 25],[0 1 0],[0.01 0.01 0.01],100);%求阶数n以及参数beta
window4=kaiser(n+1,beta);%使用kaiser窗函数
%设计加窗函数fir1
b1=fir1(N-1,wn,window1);
b2=fir1(N-1,wn,window2);
b3=fir1(N-1,wnwindow3);
b4=fir1(n,Wn,window4,'noscale');
%求取频率响应
[H1,W1]=freqz(b1,1,512,2);
[H2,W2]=freqz(b2,1,512,2);
[H3,W3]=freqz(b3,1,512,2);
[H4,W4]=freqz(b4,1,512,2);

figure(1);
subplot(2,2,1),plot(W1,20*log10(abs(H1)));%绘制频率响应图形
axis([0,1,-100,100]);
title('带通hanning窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

subplot(2,2,2),plot(W2,20*log10(abs(H2)));%绘制频率响应图形
axis([0,1,-100,100]);
title('带通hamming窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

subplot(2,2,3),plot(W3,20*log10(abs(H3)));%绘制频率响应图形
axis([0,1,-100,100]);
title('带通blackman窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值')

subplot(2,2,4),plot(W4,20*log10(abs(H4)));%绘制频率响应图形
axis([0,1,-100,100]);
title('带通kaiser窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

T=1/Fs;
L=100;%信号长度
t=(0:L-1)*T;%定义时间范围和步长
y=sin(2*pi*5*t)+5*sin(2*pi*15*t)+8*sin(2*pi*40*t);%滤波前的图形
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;%将时域信号变换到频域
f = Fs/2*linspace(0,1,NFFT/2+1);%频域采样
figure(2);
plot(f,2*abs(Y(1:NFFT/2+1)));xlabel('frequency/Hz');ylabel('Amuplitude') ;%滤波前频谱
title('滤波前的频谱');

%滤波后频谱
%采用hanning窗滤波器
yy1=filter(b1,1,y);%调用滤波函数
YY1=fft(yy1,NFFT)/L;%进行傅里叶变换，下同。
f1=Fs/2*linspace(0,1,NFFT/2+1);
figure(3);
subplot(2,2,1),plot(f1,2*abs(YY1(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('hanning窗的滤波效果');
%采用hammning窗滤波器
yy2=filter(b2,1,y);
YY2=fft(yy2,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,2),plot(f1,2*abs(YY2(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('hamming窗的滤波效果');
%采用blackman窗滤波器
yy3=filter(b3,1,y);
YY3=fft(yy3,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,3), plot(f1,2*abs(YY3(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('blackman窗的滤波效果');
%采用kaiser窗滤波器
yy4=filter(b4,1,y);
YY4=fft(yy4,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,4),plot(f1,2*abs(YY4(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
xlabel('frequency/Hz');ylabel('Amuplitude');
title('kaiser窗函数滤波效果');

%滤波前后的信号的时域对比
figure(4);
plot(y);xlabel('时间/s');ylabel('振幅');title('滤波前振幅特性');
figure(5);
subplot(2,2,1),plot(yy1);xlabel('时间/s');ylabel('振幅');title('hanning窗函数滤波振幅特性');
subplot(2,2,2),plot(yy2);xlabel('时间/s');ylabel('振幅');title('hamming窗函数滤波振幅特性');
subplot(2,2,3),plot(yy3);xlabel('时间/s');ylabel('振幅');title('blackman窗函数滤波振幅特性');
subplot(2,2,4),plot(yy4);xlabel('时间/s');ylabel('振幅');title('kaiser窗函数滤波振幅特性');
%滤波前后的信号的相位对比
figure(6);
plot(angle(Y));xlabel('时间/s');ylabel('相位');title('滤波前的相位特性');
figure(7);
subplot(2,2,1),plot(angle(YY1));xlabel('时间/s');ylabel('相位');title('hanning窗函数滤波相位特性');
subplot(2,2,2),plot(angle(YY2));xlabel('时间/s');ylabel('相位');title('hamming窗函数滤波相位特性');
subplot(2,2,3),plot(angle(YY3));xlabel('时间/s');ylabel('相位');title('blackman窗函数滤波相位特性');
subplot(2,2,4),plot(angle(YY4));xlabel('时间/s');ylabel('相位');title('kaiser窗函数滤波相位特性');

3.分别设计hanning hamming blackman kaiser窗函数highpass_FIR
function highpassfilter
clc;
clear all;
Fs=100;%采样频率
fs=35;%高通阻带模拟截止频率
fp=40;%高通通带模拟起始频率
ws=2*pi*fs/Fs;
wp=2*pi*fp/Fs;
wn=(wp+ws)/2/pi;
Bt=wp-ws;
N0=ceil(55*pi/Bt);
N=N0+mod(N0+1,2);
%调用窗函数
window1=hanning(N);
window2=hamming(N);
window3=blackman(N);
[n,Wn,beta,ftype]=kaiserord([35,40],[0 1],[0.01 0.01],100);
window4=kaiser(n1,beta);
%设计加窗函数fir1
b1=fir1(N-1,wn,'high',window1);
b2=fir1(N-1,wn,'high',window2);
b3=fir1(N-1,wn,'high',window3);
b4=fir1(n,Wn,'high',window4 ,'noscale');
%求取频率响应
[H1,W1]=freqz(b1,1,512,2);
[H2,W2]=freqz(b2,1,512,2);
[H3,W3]=freqz(b3,1,512,2);
[H4,W4]=freqz(b4,1,512,2);

figure(1);
subplot(2,2,1),plot(W1,20*log10(abs(H1)));%绘制频率响应图形
axis([0,1,-100,100]);
title('高通hanning窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

subplot(2,2,2),plot(W2,20*log10(abs(H2)));%绘制频率响应图形
axis([0,1,-100,100]);
title('高通hamming窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

subplot(2,2,3),plot(W3,20*log10(abs(H3)));%绘制频率响应图形
axis([0,1,-100,100]);
title('高通blackman窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

subplot(2,2,4),plot(W4,20*log10(abs(H4)));%绘制频率响应图形
axis([0,1,-100,100]);
title(' 高通kaiser窗的频率响应图形');
xlabel('频率（Hz)');
ylabel('幅值');

T=1/Fs;
L=100;%信号长度
t=(0:L-1)*T定义时间范围和步长
y=sin(2*pi*5*t)+5*sin(2*pi*15*t)+8*sin(2*pi*40*t);%滤波前的图形
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;%将时域信号变换到频域
f = Fs/2*linspace(0,1,NFFT/2+1);%频域采样
figure(2);
plot(f,2*abs(Y(1:NFFT/2+1)));xlabel('frequency/Hz');ylabel('Amuplitude') ;%滤波前频谱
title('滤波前的频谱');

%滤波后频谱
%采用hanning窗滤波器
yy1=filter(b1,1,y);%调用滤波函数
YY1=fft(yy1,NFFT)/L;%进行傅里叶变换，下同。
f1=Fs/2*linspace(0,1,NFFT/2+1);
figure(3);
subplot(2,2,1),plot(f1,2*abs(YY1(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('hanning窗的滤波效果');
%采用hammning窗滤波器
yy2=filter(b2,1,y);
YY2=fft(yy2,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,2),plot(f1,2*abs(YY2(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('hamming窗的滤波效果');
%采用blackman窗滤波器
yy3=filter(b3,1,y);
YY3=fft(yy3,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,3), plot(f1,2*abs(YY3(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
title('blackman窗的滤波效果');
%采用kaiser窗滤波器
yy4=filter(b4,1,y);
YY4=fft(yy4,NFFT)/L;
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(2,2,4),plot(f1,2*abs(YY4(1:NFFT/2+1))) ;xlabel('frequency/Hz');ylabel('Amuplitude');
xlabel('frequency/Hz');ylabel('Amuplitude');
title('kaiser窗函数滤波效果');

%滤波前后的信号的时域对比
figure(4);
plot(y);xlabel('时间/s');ylabel('振幅');title('滤波前振幅特性');
figure(5);
subplot(2,2,1),plot(yy1);xlabel('时间/s');ylabel('振幅');title('hanning窗函数滤波振幅特性');
subplot(2,2,2),plot(yy2);xlabel('时间/s');ylabel('振幅');title('hamming窗函数滤波振幅特性');
subplot(2,2,3),plot(yy3);xlabel('时间/s');ylabel('振幅');title('blackman窗函数滤波振幅特性');
subplot(2,2,4),plot(yy4);xlabel('时间/s');ylabel('振幅');title('kaiser窗函数滤波振幅特性');
%滤波前后的信号的相位对比
figure(6);
plot(angle(Y));xlabel('时间/s');ylabel('相位');title('滤波前的相位特性');
figure(7);
subplot(2,2,1),plot(angle(YY1));xlabel('时间/s');ylabel('相位');title('hanning窗函数滤波相位特性');
subplot(2,2,2),plot(angle(YY2));xlabel('时间/s');ylabel('相位');title('hamming窗函数滤波相位特性');
subplot(2,2,3),plot(angle(YY3));xlabel('时间/s');ylabel('相位');title('blackman窗函数滤波相位特性');
subplot(2,2,4),plot(angle(YY4));xlabel('时间/s');ylabel('相位');title('kaiser窗函数滤波相位特性');
三、运行结果
1.给定的简谐信号：

图一：输入简谐信号滤波前的频谱

图二：输入简谐信号滤波前的振幅

图三：输入简谐信号滤波前的相位

2.低通滤波器的设计
低通滤波器的技术指标：
采样频率 Fs=100Hz；通带截止频率 fp=20Hz； 阻带起始频率fs=30Hz
Hanning Hamming Blackman Kaiser采用相同的技术指标。以下即是四个窗函数的频响图及对简谐信号滤波后的效果图。
图四 不同低通窗函数低通滤波器的归一化频响图

图五 不同窗函数低通滤波器对信号的滤波后频率域效果

图六 不同窗函数低通滤波器对信号的滤波后时间域效果

图七 不同窗函数低通滤波器对信号的滤波后相位变化

由以上滤波后频率，相位，振幅变化以观察到：hanning窗 hamming窗及blackman窗的滤波效果基本相当，但三者相比：hamming窗的过渡带衰减最快，blackman窗旁瓣幅度最小。而kaiser窗只有5Hz信号，15Hz信号被截断，与设计要求有出入。

3. 带通滤波器的设计

采样频率 Fs=100Hz；阻带截止频率1：fs1=10Hz；通带起始频率1：fp1=15；通带截止频率2：fp2=20；阻带截止频率2：fs2=25Hz；
Hanning Hamming Blackman Kaiser采用相同的技术指标。以下即是四个窗函数的频响图及对简谐信号滤波后的效果图。

图八 不同带通窗函数低通滤波器的归一化频响图

图九 不同窗函数带通滤波器对信号的滤波后频率域效果
图十 不同窗函数带通滤波器对信号的滤波后时间域振幅效果

图十一 不同窗函数带通滤波器对信号的滤波后相位变化

由滤波后的相位，振幅及频率谱可看出：四个窗函数所设计的带通滤波器的滤波效果基本相当，但hanning窗对振幅谱有失真。对比图八可以看出：Blackman窗函数过渡带较窄，旁瓣幅度较小。因此，同等情况下，blackman窗函数设计带通滤波器效果较好。

4. 高通滤波器的设计

高通滤波器的技术指标：
采样频率 Fs=100Hz；阻带截止频率 fs=35Hz；通带起始频率fp=40Hz
Hanning Hamming Blackman Kaiser采用相同的技术指标。以下即是四个窗函数的频响图及对简谐信号滤波后的效果图。

图十二 不同窗函数高通滤波器的归一化频响图

图十三 不同窗函数高通滤波器对信号的滤波后频率域效果

图十四 不同窗函数高通滤波器对信号的滤波后时间域振幅效果

图十五 不同窗函数高通滤波器对信号的滤波后相位变化

由上面滤波后振幅，相位及频率谱可看出：hanning Blackman窗函数的振幅及频率特性有失真。hamming窗的过渡带较 Kaiser窗函数要陡，旁瓣幅度要小。但从频率域与振幅特性上看，kaiser效果较好。因此，设计高通滤波器kaiser窗滤波效果比较好。

四、总结体会
信号的分析与处理对我们学习地球物理的学生来说，是一个必须要掌握的知识，它对提高重，磁，电，地震信号信噪比，分辨率极其有用，同时也能用于对测井信号的分析与处理。这也是我选这个题目的原因。结合已学的《信号分析与处理》课程知识调用MATLAB信号处理的相关函数，通过窗函数设计法实现FIR滤波器的设计。一方面使我更加深刻的理解了窗函数法设计FIR滤波器的精髓，另一方面也让我认识到MATLAB功能的强大性，特别是它简洁的语句，突出的运算能力，丰富的函数库以及十分强大的绘图功能。不过，在做这个题目的过程中，我也遇到的一些困难，比如：kaiserord函数的调用方式，FFT傅立叶变换的实现，以及一些语法等细节方面的问题。通过查阅相关资料，调用MATLAB中的HELP功能，咨询老师等，最后，问题基本得到了解决。总的来时，程序仍存在一些运算速度，简洁性及稳定性方面的问题。这也是我以后需要进一步学习和加强的地方。