原创 设计FIR数字滤波器

FIR滤波器

王葵军

2010/8/4

设计一个数字FIR低通滤波器

针对上述指标，我们选择Hamming窗，拥有较小的过渡带，最小阻带衰减50dB.

wp=0.2*pi;

ws=0.3*pi;

wc=(wp+ws)/2;%截止频率

width=ws-wp;%过渡带带宽

N=ceil(6.6*pi/width)+1;%滤波器阶数

disp('the filter length:N=');

disp(N);

alpha=(N-1)/2;%求中间样本

n=[0:1N-1)];

m=n-alpha+eps;

hd=sin(wc*m)./(pi*m);%理想的冲击响应

w_hamming=hamming(N)';%HAMMING窗

h=hd.*w_hamming;%实际脉冲响应

[H,w]=freqz(h,[1],1000,'whole');%频域响应

W=w(1:501);

magH=abs(H(1:501));%幅度

angH=angle(H(1:501));%相位

dB=20*log10(magH/max(magH));

deltaw=2*pi/1000;

Rp=-(min(dB(1:1:wp/deltaw+1)))%通带波纹求最小值 也就是 最大值

disp('Rp=');

disp(Rp);

As=-round(max(dB(ws/deltaw+1:1:501)));%阻带波纹 求最大值 也就是最小值

disp('As=');

disp(As);

subplot(2,3,1)

stem(n,hd);

grid on

axis([0 66 -0.06 0.26]);

title('理想的冲击响应')

subplot(2,3,2)

stem(n,w_hamming);

grid on

axis([0 66 0 1.1]);

title('Hamming window')

subplot(2,3,3)

stem(n,h);

grid on

title('加窗的冲击响应')

axis([0 66 -0.05 0.26]);

subplot(2,3,4)

plot(W/pi,magH);

grid on

title('幅度谱')

axis([0 1 0 1.1]);

subplot(2,3,5)

plot(W/pi,angH);

grid on

title('相位谱')

axis([0 1 -3.2 3.1]);

subplot(2,3,6)

plot(W/pi,dB)

grid on

title('幅度谱dB')

axis([0 1 -106 0.2]);

按照FIR设计步骤进行设计

使用FDATOOL进行设计

设计一个数字FIR高通滤波器

在设计高通滤波器，等于全通滤波器-低通滤波器

针对上述指标，我们选择Hamming窗，拥有较小的过渡带，最小阻带衰减50dB.

close all

clear all

clc

ws=0.2*pi;%阻带频率

wp=0.3*pi;%通带起始频率

wc=(wp+ws)/2;%求的是截止频率

width=abs(ws-wp);

N=ceil(6.6*pi/width)+1;

disp('the filter length:N=');

disp(N);

alpha=(N-1)/2;

n=[0:1N-1)];

m=n-alpha+eps;

hd=sin(pi*m)./(pi*m)-sin(wc*m)./(pi*m);

w_hamming=hamming(N)';

h=hd.*w_hamming;

[H,w]=freqz(h,[1],1000,'whole');

W=w(1:501);

magH=abs(H(1:501));%幅度

angH=angle(H(1:501));%相位

dB=20*log10(magH/max(magH));

deltaw=2*pi/1000;

As=-round(max(dB(1:1:ws/deltaw+1)));

disp('As=');

disp(As);

Rp=-(min(dB(wp/deltaw+1:1:501)));

disp('Rp=');

disp(Rp);

subplot(2,3,1)

stem(n,hd);

grid on

axis([0 66 -0.3 0.76]);

title('理想的冲击响应')

subplot(2,3,2)

stem(n,w_hamming);

grid on

axis([0 66 0 1.1]);

title('Hamming window')

subplot(2,3,3)

stem(n,h);

grid on

title('加窗的冲击响应')

axis([0 66 -0.3 0.76]);

subplot(2,3,4)

plot(W/pi,magH);

grid on

title('幅度谱')

axis([0 1 0 1.1]);

subplot(2,3,5)

plot(W/pi,angH);

grid on

title('相位谱')

axis([0 1 -3.2 3.2]);

subplot(2,3,6)

plot(W/pi,dB)

grid on

title('幅度谱dB')

axis([0 1 -70 0.2]);

按照FIR设计步骤进行设计

使用FDATOOL进行设计

设计一个数字FIR带通滤波器

在设计带通通滤波器，等于截止频率为WC1的低通滤波器-截止频率为WC2的低通滤波器

针对上述指标，我们选择blackman窗，

close all

clear all

clc

wp1=0.2*pi;%

ws1=0.35*pi;%

wp2=0.65*pi;%

ws2=0.8*pi;%

wc1=(wp1+ws1)/2;%求的是截止频率

wc2=(wp2+ws2)/2;%求的是截止频率

width=min(ws2-wp2,ws1-wp1);

N=ceil(11*pi/width)+1;

disp('the filter length:N=');

disp(N);

alpha=(N-1)/2;

n=[0:1N-1)];

m=n-alpha+eps;

hd=sin(wc2*m)./(pi*m)-sin(wc1*m)./(pi*m);

w_blackman=blackman(N)';

h=hd.*w_blackman;

[H,w]=freqz(h,[1],1000,'whole');

W=w(1:501);

magH=abs(H(1:501));%幅度

angH=angle(H(1:501));%相位

dB=20*log10(magH/max(magH));

deltaw=2*pi/1000;

tempN=[1:1:wp1/deltaw,ws2/deltaw+1:1:501];

As=-round(max(dB(tempN)));

disp('As=');

disp(As);

Rp=-(min(dB(ws1/deltaw+1:1:wp2/deltaw)));

disp('Rp=');

disp(Rp);

subplot(2,3,1)

stem(n,hd);

grid on

axis([0 74 -0.32 0.46]);

title('理想的冲击响应')

subplot(2,3,2)

stem(n,w_blackman);

grid on

axis([0 74 0 1.1]);

title('Blackman window')

subplot(2,3,3)

stem(n,h);

grid on

title('加窗的冲击响应')

axis([0 74 -0.32 0.46]);

subplot(2,3,4)

plot(W/pi,magH);

grid on

title('幅度谱')

axis([0 1 0 1.1]);

subplot(2,3,5)

plot(W/pi,angH);

grid on

title('相位谱')

axis([0 1 -3.2 3.2]);

subplot(2,3,6)

plot(W/pi,dB)

grid on

title('幅度谱dB')

axis([0 1 -100 0.25]);

按照FIR设计步骤进行设计

使用FDATOOL进行设计

设计一个数字FIR带阻滤波器

设计带阻滤波器=低通滤波器1+高通滤波器2=低通滤波器1+(全通滤波器-低通滤波器2)

clc

clear all

close all

wp1=0.2*pi;

ws1=0.3*pi;

wp2=0.7*pi;

ws2=0.8*pi;

wc1=(wp1+ws1)/2;%截止频率

wc2=(wp2+ws2)/2;%截止频率

width=min(ws1-wp1,ws2-wp2);%过渡带带宽

N=ceil(6.6*pi/width)+1;%滤波器阶数

disp('the filter length:N=');

disp(N);

alpha=(N-1)/2;%求中间样本

n=[0:1N-1)];

m=n-alpha+eps;

hd=sin(wc1*m)./(pi*m)+sin(pi*m)./(pi*m)-sin(wc2*m)./(pi*m);%理想的冲击响应

w_hamming=hamming(N)';%HAMMING窗

h=hd.*w_hamming;%实际脉冲响应

[H,w]=freqz(h,[1],1000,'whole');%频域响应

W=w(1:501);

magH=abs(H(1:501));%幅度

angH=angle(H(1:501));%相位

dB=20*log10(magH/max(magH));

deltaw=2*pi/1000;

tempN=[1:1:wp1/deltaw+1,ws2/deltaw:1:501];

Rp=-(min(dB(tempN)));%通带波纹

disp('Rp=');

disp(Rp);

As=-round(max(dB(ws1/deltaw+1:1:wp2/deltaw)));%阻带波纹

disp('As=');

disp(As);

subplot(2,3,1)

stem(n,hd);

grid on

axis([0 66 -0.11 0.51]);

title('理想的冲击响应')

subplot(2,3,2)

stem(n,w_hamming);

grid on

axis([0 66 0 1.1]);

title('Hamming window')

subplot(2,3,3)

stem(n,h);

grid on

title('加窗的冲击响应')

axis([0 66 -0.1 0.51]);

subplot(2,3,4)

plot(W/pi,magH);

grid on

title('幅度谱')

axis([0 1 0 1.1]);

subplot(2,3,5)

plot(W/pi,angH);

grid on

title('相位谱')

axis([0 1 -3.2 3.1]);

subplot(2,3,6)

plot(W/pi,dB)

grid on

title('幅度谱dB')

axis([0 1 -70 0.2]);

按照FIR设计步骤进行设计

使用FDATOOL进行设计

attachment download

pwd:123456

close all;
clear all;
clc;

f1=1900;
f2=5100;
f3=7000;
Fs=20000;
N=1000;
n=0:1/FsN-1)/Fs;
y=sin(2*pi*f1.*n)+cos(2*pi*f2.*n)+sin(2*pi*f3.*n);
plot(n,y)
xlabel('Time /sec');
ylabel('Amplitude volts');
title('the signal in Time Domain before filter');

Y=fft(y,N);
F=(0:N-1)*Fs/N;
figure
plot(F,abs(Y))%1
xlabel('Frequency /Hz');
ylabel('Amplitude volts');
title('the signal in Frequency Domain before filter');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%the low pass filter design Using MATLAB language
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f1=1900;%this signal frequency will be passed in the filter
f2=5100;%this signal frequency will be filtered in the filter
f3=7000;%this signal frequency will be filtered in the filter
Fs=20000;% the sample frequency

w1=2*pi*f1/Fs %normalization in digital frequency domain
%w1/pi=0.19;
%we can define the pass frequency and stop frequency because we know the
%frequency of the signal
%also we can get the cutoff frequency
wp=0.2*pi;
ws=0.25*pi;
wc=(wp+ws)/2;
width=ws-wp;
%get the filter order
order=ceil(6.6*pi/width)+1;
window=hamming(order);
b=fir1(order-1,wc/pi,window);
figure
freqz(b,1,512)
%filter
sr = filter(b,1,y);
figure
plot(sr)
xlabel('Time /sec');
ylabel('Amplitude volts');
title('the signal in Time Domain after filter');figure
plot(F,abs(fft(sr)))
xlabel('Frequency /Hz');
ylabel('Amplitude volts');
title('the signal in Frequency Domain before filter');