matlab简单分析数字滤波器FIR

时域FIR滤波器

华中科技大学《数字信号分析理论实践》第六单元 FIR滤波器 学习总结记录

FIR 滤波器脉冲响应函数设计法

• MATLAB 中 FIR 滤波器设计函数 fir2()

f = [0 0.2 0.25 1]; % 滤波器频率
m = [1 1 0.2 0]; % 滤波器幅值
% 归一化频率
b = fir2(20,f,m); % 20个点
[h,w] = freqz(b,1,128);
subplot(211);
plot(f,m,w/pi,abs(h),'linewidth',3);
legend('设计参数','实际滤波器')
ylim([-0.1 1.1])
subplot(212);
impz(b,1)
axis tight
ylim([-0.1 0.4])

• 先定义一个滤波器，再给一个激励，对脉冲响应进行离散采样，得到滤波器系数
  
• 减小截断带来的误差——窗函数设计法
  

N = 61;
w = window(@hamming,N);
subplot(411);
plot(w,'linewidth',2);
title('hamming')
w = window(@hanning,N);
subplot(412);
plot(w,'linewidth',2);
title('hanning')
w = window(@blackmanharris,N);
subplot(413);
plot(w,'linewidth',2);
title('blackmanharris')
w = window(@tukeywin,N,0.5);
subplot(414);
plot(w,'linewidth',2);
title('tukeywin')

% 脉冲响应 + 窗函数 设计 FIR 滤波器
b = fir1(48,[0.35 0.65],blackman(48+1));
% 系数长度 48 截止频率 频率通带为归一化频率 滤波器类型 带通滤波器 窗类型 blackman
% 窗函数长度比滤波器系数长度多1
freqz(b,1,512)
% 画出实际滤波器的频率特性

• 滤波器系数越长，暂态部分越大，直接调用 filter() 但还是有头部的暂态
  

Fs = 2048;
dt = 1.0/Fs;
T = 1;
N = T/dt;
t = linspace(0,T,N);
x1 = sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
subplot(311);
plot(t,x1);
axis([0,0.1,-2,2]);
P = fft(x1,N);
Pyy = 2*sqrt(P.*conj(P))/N;
f = linspace(0,Fs/2,N/2);
subplot(312)
plot(f,Pyy(1:N/2));
b = fir1(48,0.1);
% freqz(b,1,128)
x2 = filter(b,1,x1);
subplot(313)
plot(t,x2);
axis([0 0.1 -2 2])
figure
subplot(311);
plot(t,x1);
axis([0,0.1,-2,2]);
b = fir1(48,[0.2,0.4]);
x2 = filter(b,1,x1);
subplot(312)
plot(t,x2);
axis([0 0.1 -2 2])
b = fir1(48,0.4,'high');
x2 = filter(b,1,x1);
subplot(313)
plot(t,x2);
axis([0 0.1 -2 2])

• FIR2() 窗函数法，增加过渡点，性能好于FIR1()

f = [0 0.2 0.3 0.4 1]; % 滤波器频率
m = [1 1 0.7 0 0]; % 滤波器幅值
% 归一化频率，1代表采样频率的一半
b = fir2(48,f,m); % 20个点
[h,w] = freqz(b,1,128);
plot(f,m,'o-',w/pi,abs(h),'linewidth',3);
legend('设计参数','实际滤波器')

• 复杂频率通带滤波器设计，定义两个频率通带
  

f = [0 0.1 0.2 0.4 0.5 0.6 1]; % 滤波器频率
m = [1 1 0 0 1 1 0]; % 滤波器幅值
% 归一化频率，1代表采样频率的一半
b = fir2(48,f,m); % 20个点
[h,w] = freqz(b,1,128);
plot(f,m,'o-',w/pi,abs(h),'linewidth',3);
legend('设计参数','实际滤波器')

Fs = 2048;
dt = 1.0/Fs;
T = 1;
N = T/dt;
t = linspace(0,T,N);
x1 = sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
subplot(311);
plot(t,x1);
axis([0,0.1,-2,2]);
P = fft(x1,N);
Pyy = 2*sqrt(P.*conj(P))/N;
f = linspace(0,Fs/2,N/2);
subplot(312)
plot(f,Pyy(1:N/2));
f = [0 0.2 0.2 1];
m = [1 1 0 0];
b = fir2(48,f,m);
% freqz(b,1,128)
x2 = filter(b,1,x1);
subplot(313)
plot(t,x2);
axis([0 0.1 -2 2])

• 最小二乘法滤波器设计 FIRLS()，设计的滤波器特性曲线在指定点的误差平方和最小
  

F = [0 0.2 0.4 0.6 0.8 0.9]; 
A = [1 1 0.6 0.4 0.2 0];
b = firls(24,F,A); 
[h,f] = freqz(b,1,512,2);
plot(F,A,'o-',f,abs(h),'linewidth',3);
legend('设计参数','实际滤波器')

• 最小二乘法可能在未定义点有误差
  

F = [0 0.1 0.2 0.3 0.4 0.5 0.6 1]; 
A = [1 1 0 0 0 1 1 0];
b = firls(48,F,A); 
[h,w] = freqz(b,1,128);
subplot(211)
plot(F,A,'o-',w/pi,abs(h),'linewidth',3);
legend('设计参数','实际滤波器')
b = fir2(48,F,A);
subplot(212)
[h,w] = freqz(b,1,128);
plot(F,A,'o-',w/pi,abs(h),'linewidth',3);
legend('设计参数','实际滤波器')

• 鼓包可能导致频率无法滤除干净
  

Fs = 2048;
dt = 1.0/Fs;
T = 1;
N = T/dt;
t = linspace(0,T,N);
x1 = sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
subplot(411);
plot(t,x1);
axis([0,0.1,-2,2]);
P = fft(x1,N);
Pyy = 2*sqrt(P.*conj(P))/N;
f = linspace(0,Fs/2,N/2);
subplot(412)
plot(f,Pyy(1:N/2));
f = [0 0.2 0.2 1];
%-------设计滤波器--------%
F = [0 0.1 0.4 0.6 0.8 1]; 
A = [1 1 0 0 0 0];
b = firls(48,F,A);
[h,f] = freqz(b,1,512,2);
subplot(413)
plot(F,A,'o-',f,abs(h),'linewidth',3);
legend('设计参数','实际滤波器')
x2 = filter(b,1,x1);
subplot(414)
plot(t,x2);
axis([0 0.1 -2 2])