实验方法及数字信号分析处理-滤波器设计

2014 年春季学期研究生课程考核大作业一

考核科目:实验方法及数字信号分析处理学生所在院（系）:

学生所在学科:

学生姓名:

学号:

学生类别:

考核结果阅卷人

一、题目：

给定信号()sin(25)sin(250)sin(2300)x t t t t πππ=?+?+?。 要求：

1、离散化（t ?），画出{k x }曲线（满足采样定理）

2、用卷积滤波器滤波

1) 低通：保留5Hz 信号（带修正/不带） 2) 带通：保留50Hz 信号

3、上述信号加上白噪声（0.2rand （1，v ））再滤波

二、图像及matlab 程序

1、离散化（t ?），画出{k x }曲线（满足采样定理）

根据采样定理：max

1

2t f ?≤

，max f =300Hz ， 则t ?≤0.00167，取采样周期t ?=0.001，满足采样定理。

在MATLAB 中，用横坐标表示时间（取时间长度为0.4s ），纵坐标表示幅值，原信号如下图所示：

-200

-150-100-50

50100150200

-3-2

-1

1

2

3

原信号离散化图形

时间t(ms)

信号幅值

离散化matlab 程序如下：

clear; clc;

Dt=0.001; t=-0.2:Dt:0.2;

x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*300*t); figure(1);

plot(t*1000,x); grid on;

title('原信号离散化图形'); xlabel('时间t(ms)'); ylabel('信号幅值');

2、用卷积滤波器滤波

a.低通：保留5Hz 信号（带修正/不带）

FIR 滤波器传递函数为：0?()N

i i

i H z f z -==∑其中，f i 滤波因子，N 为滤波长度， 对于低通和带通时其滤波因子求取如下： {f i }结果：

sin(2)

i 0 f i Fi t i

ππ≠=

时； i i 0f =2F t =时 ；取截止频率F=20Hz.

图像中包含五条曲线，分别表示标准信号','原始信号','未修正低通滤波','汉宁修正低

通滤波','汉明修正低通滤波，第一幅图是整体的图像，第二幅图为细节图。

分析对比：

观察第三幅细节图大致可知，卷积低通滤波效果较好，相位基本没有滞后，但是未经过修正的信号因泄露原因幅值较小，且信号曲线不太平滑。汉明修正结果较汉宁修正结果要好，更符合标准曲线。

-0.2

-0.15-0.1-0.0500.050.10.150.2

-3-2

-1

1

2

3

时间t(ms))

信号幅值）

低通滤波器滤波及对比

0.04

0.045

0.05

0.055

0.06

0.8

0.850.90.9511.05

1.1时间t(ms))

信号幅值）

低通滤波器滤波及对比

低通滤波器滤波程序如下：

主程序：

clc;

clear;

Dt=1/1000;%设置采样时间间隔

t=-0.2:Dt:0.2;%设置采样时间长度

x1=sin(2*pi*5*t);

x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*300*t);%输入信号（未加噪声）

%x=x+0.2*rand(1,size(x,2));%在原信号上加上一个均值为0，方差为1，序列长度与原信号等长的白噪声

figure(2)

set( figure(2), 'Color','w');

[y,tl]=Lowpass(70,70,20,1);%原信号通过低通滤波器后返回的信号值y及时间长度，N1=70，N2=70，截止频率F=20

[y1,tl]=Lowpass(70,70,20,2);

[y2,tl]=Lowpass(70,70,20,3);

plot(t,x1,t,x,tl,y,tl,y1,'g--',tl,y2,'r-.');

hold on;

grid on;

xlabel('时间t(ms))','Fontsize',18);

ylabel('信号幅值）','Fontsize',18);

title('低通滤波器滤波及对比','Fontsize',18);

legend('标准信号','原始信号','未修正低通滤波','汉宁修正低通滤波','汉明修正低通滤波');

功能函数：

function [ y,tl ] = Lowpass(N1,N2,F,cor)

%低通卷积滤波器,cor--修正参数设置，1—不修正2—汉宁修正3—汉明修正

Dt=0.001;

t=-0.2:Dt:0.2;

x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*300*t);

for i=-N1:1:N2

if i ~=0

switch cor

case 1

f(i+N1+1)=sin(2*pi*F*i*Dt)/(pi*i);

case 2

f(i+N1+1)=1/2*(1+cos(pi*(i)/N1))*sin(2*pi*F*i*Dt)/(pi*i);

case 3

f(i+N1+1)=(0.54+0.46*cos(pi*(i)/N1))*sin(2*pi*F*i*Dt)/(pi*i);

end

end if i==0

f(i+N1+1)=real(2*F*Dt); end end

for k=(N2+1):1:(size(x,2)-N1-1) y(k-N2)=0; for i=-N1:1:N2

y(k-N2)=f(i+N1+1)*x(k-i)+y(k-N2); end end

tl=t(N2+1):Dt:t(size(t,2)-N1-1); end

b.带通：保留50Hz 信号（带修正/不带）

FIR 滤波器传递函数为：0?()N

i i

i H z f z -==∑其中，f i 滤波因子，N 为滤波长度， 对于低通和带通时其滤波因子求取如下：

{f i }结果：

02cos(2)sin(2)

i 0 f i F i t Fi t i

πππ≠=时 i i 0f =4t F =时

其中1202F F F +=，12()2F F

F abs -=；取1235Hz 65Hz F F ==；

则

050F = Hz ，15F = Hz 。

图像中包含五条曲线，分别表示标准信号','原始信号','未修正带通滤波','汉宁修正带

通滤波','汉明修正带通滤波，第一幅图是整体的图像，第二幅图为部分放大图，第三幅图为细节图。

分析对比：

观察第三幅细节图大致可知，和低通滤波结果一致，卷积带通滤波效果较好，相位基本没有滞后，但是未经过修正的信号因泄露原因幅值较小，且信号曲线不太平滑。汉明修正结果较汉宁修正结果要好，更符合标准曲线。

-0.2

-0.15-0.1-0.0500.050.10.150.2

-3-2

-1

1

2

3

时间t(ms))

信号幅值）

带通滤波器滤波及对比

-0.06-0.05-0.04-0.03-0.02-0.0100.010.020.03

0.4

0.60.811.21.4

1.6时间t(ms))

信号幅值）

带通滤波器滤波及对比

带通滤波器滤波程序如下： 主程序： clc; clear;

Dt=1/1000;%设置采样时间间隔 t=-0.2:Dt:0.2;%设置采样时间长度 x1=sin(2*pi*50*t)

x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*300*t);%输入信号（未加噪声）

%x=x+0.2*rand(1,size(x,2));%在原信号上加上一个均值为0，方差为1，序列长度与原信号等长的白噪声 figure(1)

set( figure(1), 'Color','w');

[y,tb]=Bandpass(70,70,35,65,1);%原信号通过带通滤波器后返回的信号值y 及时间长度 [y1,tb]=Bandpass(70,70,35,65,2); [y2,tb]=Bandpass(70,70,35,65,3); plot(t,x1,t,x,tb,y,tb,y1,'g--',tb,y2,'r-.'); hold on; grid on;

-0.08-0.078-0.076-0.074-0.072-0.07

0.84

0.860.880.90.920.940.960.98

11.02时间t(ms))

信号幅值）

带通滤波器滤波及对比

xlabel('时间t(ms))','Fontsize',18);

ylabel('信号幅值）','Fontsize',18);

title('带通滤波器滤波及对比','Fontsize',18);

legend('标准信号','原始信号','未修正带通滤波','汉宁修正带通滤波','汉明修正带通滤波');

功能函数：

function [y,tb]=Bandpass(N1,N2,F1,F2,cor)

%带通滤波器功能函数，N1+N2+1滤波长度，F1，F2截止频率上下限，cor,修正参数，1—不修正2—汉宁修正3—汉明修正

Dt=0.001;

t=-0.2:Dt:0.2;

x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*300*t);

F0=(F1+F2)/2;

DF=(F2-F1)/2;

for i=-N1:1:N2

if i~=0

switch cor

case 1

f(i+N1+1)=(2*sin(2*pi*DF*i*Dt)*cos(2*pi*F0*i*Dt)/(pi*i));

case 2

f(i+N1+1)=1/2*(1+cos(pi*(i)/N1))*(2*sin(2*pi*DF*i*Dt)*cos(2*pi*F0*i*Dt)/(pi*i));

case 3

f(i+N1+1)=(0.54+0.46*cos(pi*(i)/N1))*(2*sin(2*pi*DF*i*Dt)*cos(2*pi*F0*i*Dt)/(pi*i));

end

end

if i==0

f(i+N1+1)=4*Dt*DF;

end

end

for k=(N2+1):1:(size(x,2)-N1-1)

y(k-N2)=0;

for i=-N1:1:N2

y(k-N2)=f(i+N1+1)*x(k-i)+y(k-N2);

end

end

tb=t(N2+1):Dt:t(size(t,2)-N1-1);

end

c. 用递推滤波器滤波：保留5Hz 信号。

低通滤波器设计——基于巴特沃斯（Butterworth ）函数，这里采用二阶递推滤波器进行低通滤波，参数计算公式如下：

222

21

tan()()1()(j f t j f t e F t s T z e B s s j

s ππωπωω

ω

ω

ω

-===+=

+

+；；

H(z)=；取

F=20

21f =

；22f =

；2

3f =；

2121g ω-=

2g = 图像中包含五条曲线，分别表示,'原始信号', 标准信号'递推低通滤波，第一幅图是整体的图像，第二幅图为为细节图。本次选用两个截止频率F=10，和F=20进行滤波，滤波结果如下图所示。

分析对比：

观察第三幅细节图大致可知，递推滤波器滤波时存在明显的相位滞后，但是幅值衰减并不明显，当截止频率F=20时，存在明显的噪声干扰，F=10时，噪声干扰较小，说明递推滤波器的频率响应在过了截止频率以后衰减并不是很迅速，因此选用递推滤波器时应该着重考虑其相位的滞后作用和抗干扰性能。

-0.2

-0.15-0.1-0.0500.050.10.150.2

-3-2

-1

1

2

3

时间t(ms))

信号幅值）

递推低通滤波器滤波及对比F=20Hz

时间t(ms))

信号幅值）

递推低通滤波器滤波及对比F=20Hz

-0.2

-0.15-0.1-0.0500.050.10.150.2

-3-2

-1

1

2

3

时间t(ms))

信号幅值）

递推低通滤波器滤波及对比F=10Hz

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

0.11

0.2

0.40.60.811.2

1.4时间t(ms))

信号幅值）

递推低通滤波器滤波及对比F=10Hz

递推滤波器的matlab程序如下：

主程序：

clear;

clc;

Dt=0.001;%设置采样时间间隔

t=-0.2:Dt:0.2;%设置采样时间长度

x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*300*t);%输入信号（未加噪声）

x1=sin(2*pi*5*t)

%x=x+0.2*rand(1,size(x,2));%在原信号上加上一个均值为0，方差为1，序列长度与原信号等长的白噪声

figure(1)

set( figure(1), 'Color','w');

[y,ti]=IIRLowpass(10);%截止频率是10Hz

plot(t,x,t,x1,ti,y);

legend('原始信号','标准信号','递推低通滤波');

hold on;

grid on;

xlabel('时间t(ms))','Fontsize',18);

ylabel('信号幅值）','Fontsize',18);

title('递推低通滤波器滤波及对比F=10Hz','Fontsize',18);

功能函数：

function [y,ti]=IIRLowpass(F)

Dt=0.001;%设置采样时间间隔

t=-0.2:Dt:0.2;%设置采样时间长度

x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*300*t);%输入信号（未加噪声）

w=tan(pi*Dt*F);

f(1)=w*w/(1+sqrt(2)*w+w*w);

f(2)=2*w*w/(1+sqrt(2)*w+w*w);

f(3)=w*w/(1+sqrt(2)*w+w*w);

g(1)=-2*(1-w*w)/(1+sqrt(2)*w+w*w);

g(2)=(1-sqrt(2)*w+w*w)/(1+sqrt(2)*w+w*w);

y(1)=f(1)*x(1);

y(2)=f(1)*x(2)+f(2)*x(1)-g(1)*y(1);

ti=t(1):Dt:t(size(t,2));

for k=3:1:size(x,2)

y(k)=f(1)*x(k)+f(2)*x(k-1)+f(3)*x(k-2)-g(1)*y(k-1)-g(2)*y(k-2);

end

end

3.上述信号加上白噪声（0.2rand （1，v ））再滤波

原始信号为x=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*300*t);加入随机白噪声后信号变为x=x+0.2*rand(1,size(x,2));即在原信号上加上一个均值为0，方差为1，序列长度与原信号等长的随机白噪声，FIR 与IIR 的滤波程序同之前相同，故不再列出，其滤波图线如下所示。

如上图可知，加入随机噪声后其信号波动更为明显。

-0.2

-0.15-0.1-0.0500.050.10.150.2

-3-2

-1

1

2

3

4时间t(ms))

信号幅值）

加入白噪声后的信号

时间t(ms))

信号幅值）

加入白噪声后的信号

对于有随机噪声的信号，低通滤波器中经过修正和未修正的滤波都有幅值衰减，但相位滞后不明显。

-0.2

-0.15-0.1-0.0500.050.10.150.2

-3-2

-1

时间t(ms))

信号幅值）

低通滤波器滤波及对比

-0.07-0.065

-0.06-0.055

-0.05-0.045-0.04-0.035-0.03-0.025

时间t(ms))

信号幅值）

低通滤波器滤波及对比

对于有随机噪声的信号，带通滤波器中未修正的滤波都有幅值衰减，经过休整的信号衰减程度较小，所有滤波相位滞后不明显。

时间t(ms))

信号幅值）

带通滤波器滤波及对比

时间t(ms))

信号幅值）

带通滤波器滤波及对比

可以看出递推低通滤波器有明显的相位滞后作用，且F=10的滤波效果较F=20的滤波效果好的多，说明其递推滤波器的频率响应在过了截止频率以后衰减并不是很迅速，因此选用递推滤波器时应该着重考虑其相位的滞后作用和抗干扰性能。

-0.2

-0.15-0.1-0.0500.050.10.150.2

-3-2

-1

1

2

3

4

时间t(ms))

信号幅值）

递推低通滤波器滤波及对比F=10Hz

-0.2

-0.15-0.1-0.0500.050.10.150.2

时间t(ms))

信号幅值）

递推低通滤波器滤波及对比F=20Hz