[Matlab]FIR滤波器设计:(FIR滤波器的结构)
FIR(Finite Impulse Response)滤波器:有限长单位冲激响应滤波器,又称为非递归型滤波器,是一种在数字信号领域应用非常广泛的基础型滤波器。FIR滤波器的特点是能够在输入具体任意幅频特性的数字信号后,保证输出数字信号的相频特性仍然保持严格线性。另外FIR滤波器是数字信号处理系统中最基本的元件,它可以在保证任意幅频特性的同时具有严格的线性相频特性,同时具有有限长的脉冲采样响应特性,因而滤波器是稳定的系统。因此,FIR滤波器的应用要远胜于IIR滤波器,在通信、图像处理、模式识别等领域都有着广泛而举足轻重的作用。
FIR滤波器的结构:
尽管IIR数字滤波器的设计理论已经非常成熟、经典,但是IIR滤波器有一个自身的中要缺陷:相位特性通常是非线性的。
- IIR滤波器的非线性原理:在IIR滤波器的滤波器设计中,只是对滤波器的幅频特性进行了研究,并获得了良好的幅频特性。但对相频特性却没有考虑,所以IIR滤波器的相频特性通常是非线性的。
- IIR滤波器的应用问题:由于IIR滤波器的相位特性通常是非线性的,所以它适合用于在对系统相位特性要求不严格的场合。但是另一方面,由于一个具有线性相位特性的滤波器可以保证在滤波器通带内信号传输不失真,所以在许多领域需要滤波器具有严格得线性相位特性,比如图像处理以及数据传输等。但是如果要使用IIR滤波器实现线性的相位特性,则必须对其相位特性用全通滤波器进行校正,其结果使得滤波器设计变得非常复杂,实现困难成本提高。所以,要实现具有线性相位的数字滤波器,必须另辟蹊径。
- 有限脉冲响应系统的单位脉冲响应 h ( s ) h(s) h(s)为有限长序列,系统函数 H ( z ) H(z) H(z)在有限z平面上不存在几点,其运算结构中没有反馈支路,即没有环路。所以,有限脉冲响应滤波器可以设计成在整个频率范围内均可以提供精确的线性相位,而且总是可以独立于滤波器系数保持有限长输入有限长输出的稳定。因此,这样的滤波器在很多领域是首选的。
FIR滤波器的特点
有限长单位冲激响应(FIR)滤波器有以下特点:
- 系统的单位冲激响应h (n)在有限个n值处不为零
- 系统函数H(z)在|z|>0处收敛,极点全部在z = 0处(因果系统)
- 结构上主要是非递归结构,没有输出到输入的反馈
FIR滤波器的基本结构:
- 直接型:
- 级联型:
特点:
- 每个基本节控制一对零点便于控制滤波器的传输零点
- 系数比直接型多所需的乘法运算多。
直接型与级联型:
clear all;
n=0:10;
N=30;
b=0.9.^n;
delta=impseq(0,0,N);%脉冲系列生成器
h=filter(b,1,delta);%其中b是分子系数向量,a是分母系数向量。
x=[ones(1,5),zeros(1,N-5)];
y=filter(b,1,x);
subplot(2,2,1);stem(h);%绘制火柴图
title('直接型h(n)');
subplot(2,2,2);stem(y);
title('直接型y(n)');
[b0,B,A]=dir2cas(b,1);
h=casfilter(b0,B,A,delta);%级联型单位脉冲响应
y=casfilter(b0,B,A,x);
subplot(2,2,3);stem(h);
title('级联型h(n)');
subplot(2,2,4);stem(y);
title('级联型y(n)');
%%impseq函数设计生成m文件放入当前文件夹:
function [x,n] = impseq(n0,n1,n2)
% Generates x(n) = delta(n-n0); n1 <= n,n0 <= n2
% ----------------------------------------------
% [x,n] = impseq(n0,n1,n2)
%
if ((n0 < n1) | (n0 > n2) | (n1 > n2))
error('arguments must satisfy n1 <= n0 <= n2')
end
n = [n1:n2];
%x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))];
x = [(n-n0) == 0];
%%dir2cas函数设计生成m文件放入当前文件夹:
function [b0,B,A] = dir2cas(b,a)
% 直接型到级联型的型式转换(复数对型)
% [b0,B,A] = dir2cas(b,a)
% b = 直接型的分子多项式系数
% a = 直接型的分母多项式系数
% b0 = 增益系数
% B = 包含各bk的K乘3维实系数矩阵
% A = 包含各ak的K乘3维实系数矩阵
% compute gain coefficient b0
b0 = b(1); b = b/b0;
a0 = a(1); a = a/a0;
b0 = b0/a0;
%
M = length(b); N = length(a);
if N > M
b = [b zeros(1,N-M)];
elseif M > N
a = [a zeros(1,M-N)]; N = M;
else
NM = 0;
end
%
K = floor(N/2); B = zeros(K,3); A = zeros(K,3);
if K*2 == N;
b = [b 0];
a = [a 0];
end
%
broots = cplxpair(roots(b));
aroots = cplxpair(roots(a));
for i=1:2:2*K
Brow = broots(i:1:i+1,:);
Brow = real(poly(Brow));
B(fix((i+1)/2),:) = Brow;
Arow = aroots(i:1:i+1,:);
Arow = real(poly(Arow));
A(fix((i+1)/2),:) = Arow;
end
%%casfilter函数设计生成m文件放入当前文件夹:
function y = casfilter(b0,B,A,x)
% CASCADE form realization of IIR and FIR filters
% -----------------------------------------------
% y = casfiltr(b0,B,A,x);
% y = output sequence
% b0 = gain coefficient of CASCADE form
% B = K by 3 matrix of real coefficients containing bk's
% A = K by 3 matrix of real coefficients containing ak's
% x = input sequence
%
[K,L] = size(B);
N = length(x);
w = zeros(K+1,N);
w(1,:) = x;
for i = 1:1:K
w(i+1,:) = filter(B(i,:),A(i,:),w(i,:));
end
y = b0*w(K+1,:);
- 快速卷积型:
运用FFT实现有限长序列 x ( n ) x(n) x(n)与 h ( x ) h(x) h(x)的线性卷积,则可以得到FIR滤波器的快速卷积结构如图所示,对于 x ( n ) x(n) x(n)的无限长序列的一般情况,可以采用重叠相加法或者重叠保留法实现FIR滤波器的快速卷积结构
- 频率采样型:
FIR系统直接型转换为频率采样型结构的Matlab实现源代码:
function [C,B,A] = dir2fs(h)
% 直接型到频率采样型的转换
% [C,B,A] = dir2fs(h)
% C = 包含各并行部分增益的行向量
% B = 包含按行排列的分子系数矩阵
% A = 包含按行排列的分母系数矩阵
% h = FIR滤波器的脉冲响应向量
M = length(h);
H = fft(h,M);
magH = abs(H); phaH = angle(H)';
% check even or odd M
if (M == 2*floor(M/2))
L = M/2-1; % M为偶数
A1 = [1,-1,0;1,1,0];
C1 = [real(H(1)),real(H(L+2))];
else
L = (M-1)/2; % M is odd
A1 = [1,-1,0];
C1 = [real(H(1))];
end
k = [1:L]';
% 初始化 B 和 A 数组
B = zeros(L,2); A = ones(L,3);
% 计算分母系数
A(1:L,2) = -2*cos(2*pi*k/M); A = [A;A1];
% 计算分子系数
B(1:L,1) = cos(phaH(2:L+1));
B(1:L,2) = -cos(phaH(2:L+1)-(2*pi*k/M));
% 计算增益系数
C = [2*magH(2:L+1),C1]';
案例分析:
%[例6-1]
%%利用频率采样法设计一个低通FIR数字低通滤波器,其理想频率特性是矩形的
%pi=3.14
%给定抽样频率为Ωs=2*pi*1.5*10^4(rad/sec),通带截至频率为Ωp=2*pi*1.6*10^3(rad/sec),
%阻带截至频率为Ωst=2*pi*3.1*10^3(rad/sec),通带波动<=1dB ,阻带衰减>=50dB。
%通带的截至频率为:w_p = Ωp/f_s = 2*pi Ωp/Ωs = 0.213*pi
%阻带的起始频率为:w_st = Ωst/f_s = 2*pi Ωfs/Ωs = 0.413*pi
%理想低通截至频率:Ωc = 1/2(Ωp+Ωst) = 2*pi*2.35*10^3(rad/sec)
%其对应的数字频率:w_c = 2*pi Ωc/Ωs 0.313*pi
%过渡带带宽为:△w = w_st - w_p =0.2*pi 由于△w = (2*pi/N) *3 所以抽样频率为30rad/sec
% 直接型到频率采样型的转换
close all;
clear;
N=30;
H=[ones(1,4),zeros(1,22),ones(1,4)];
H(1,5)=0.5886;H(1,26)=0.5886;H(1,6)=0.1065;H(1,25)=0.1065;
k=0:(N/2-1);k1=(N/2+1):(N-1);k2=0;
A=[exp(-j*pi*k*(N-1)/N),exp(-j*pi*k2*(N-1)/N),exp(j*pi*(N-k1)*(N-1)/N)];
HK=H.*A;
h=ifft(HK);
fs=15000;
[c,f3]=freqz(h,1);
f3=f3/pi*fs/2;
subplot(221);
plot(f3,20*log10(abs(c)));
title('频谱特性');
xlabel('频率/HZ');ylabel('衰减/dB');
grid on;
subplot(222);
title('输入采样波形');
stem(real(h),'.');
line([0,35],[0,0]);
xlabel('n');ylabel('Real(h(n))');
grid on;
t=(0:100)/fs;
W=sin(2*pi*t*750)+sin(2*pi*t*3000)+sin(2*pi*t*6500);
q=filter(h,1,W);
[a,f1]=freqz(W);
f1=f1/pi*fs/2;
[b,f2]=freqz(q);
f2=f2/pi*fs/2;
subplot(223);
plot(f1,abs(a));
title('输入波形频谱图');
xlabel('频率');ylabel('幅度')
grid on;
subplot(224);
plot(f2,abs(b));
title('输出波形频谱图');
xlabel('频率');ylabel('幅度')
grid on;