窗函数法设计滤波器并且matlab/C语言实现

本文介绍滤波器的设计方法与滤波器的实现过程。

基于窗函数的方法是最为简单的设计方法，并且实际性能也还可以，并且有着计算量小的优点，在工程中应用广泛。本文介绍两种：多个矩形窗/三角窗及联/并联的方法设计滤波器。

本文中统一采样率为4Hz，另外采用截止波长lambda代替滤波器设计参数截止频率f，两者的关系为lambda = 1/f。

矩形窗的系统函数可以写成：
$$W(z)=\frac{1}{2M+1}\frac{z^M-z^{-M-1}}{1-z^{-1}}$$

$$W(\Omega )=\frac{1}{N}\frac{\sin (N\Omega /2)}{\sin (\Omega /2)},N=2M+1$$

对应三角窗表示为：
$$\begin{array}{l}W(z)=\frac{z^M-2+z^{-M}}{M^2{\left(1-z^{-1}\right)}^2}\\ W(\Omega )=\frac{1}{M^2}{\left(\frac{\sin (M\Omega /2)}{\sin (\Omega /2)}\right)}^2\end{array}$$
下面设计一个两个矩形窗及联的低通滤波器Fir1，两个滤波器的窗长分别是101和81。分析幅频特性的方法参考下面：

lamda = 1:0.1:10000;
pesi  = 1./lamda;           %空间域频率
omega = 2*pi*pesi*0.25;
figure1 = figure('Color',[1 1 1]);
semilogx(lamda, abs(...
  1.*wincau(101,omega).*wincau(81,omega)),'LineWidth',1);hold on;
function Tri = triCau(N,omega)
M1 = (N-1)/2;
for i=1:length(omega)
Tri(i)= ((sin(omega(i)*M1/2)/sin(omega(i)/2))/M1)*((sin(omega(i)*M1/2)/sin(omega(i)/2))/M1);
end
end
function Rec = wincau(N3 , omega)
for i=1:length(omega)
Rec(i) = ((sin(omega(i)*N3/2)/sin(omega(i)/2))/N3);
end
end

当然也可以用这种方法设计高通滤波器：

semilogx(lamda, abs(...
  1-1.*wincau(101,omega).*wincau(81,omega)),'LineWidth',1);hold on;

下面介绍在程序中的实现方法，设计滤波器形式为：3个矩形窗及联为fir1，2个矩形窗及联为fir2，再将fir1与fir2并联，得到一个滤波器形式如下：

在滤波器的实现过程中，需要考虑数组存放的问题，以及并联的滤波器的对齐的问题。实现代码参考如下：

clear all;
x = 0:0.25:1e4;
lambda = 62;
amcol = sin(2*pi/lambda*x);

%% 参数配置
NRec_1 = 201;
NRec_2 = 101;
NRec_3 = 251;
TriM1 = (NRec_1-1)/2;
TriM2 = (NRec_2-1)/2;
TriM3 = (NRec_3-1)/2;
N = 2048;
data1 = zeros(N,1);
data2 = zeros(N,1);
data3 = zeros(N,1);
center_pos1 = 1000;
center_pos2 = center_pos1 - TriM2;
center_pos3 = center_pos2 - TriM3;
RecSum1 = 0;
RecSum2 = 0;
RecSum3 = 0;

%% 并联的滤波器的参数配置
SecN1 = 45;
SecN2 = 89;
SecM1 = (SecN1-1)/2;
SecM2 = (SecN2-1)/2;
Seccenter_pos1 = center_pos3 + SecM2;
Seccenter_pos2 = center_pos3;

Seccenter_pos1 = 1000;
Seccenter_pos2 = Seccenter_pos1 - SecM2;
Secdata1 = zeros(N,1);
Secdata2 = zeros(N,1);
SecRecSum1 = 0;
SecRecSum2 = 0;


%% 
for i = 1:length(amcol)
   %% 第一层
   input_index = center_pos1 + TriM1;
   data1( find_index(input_index)) = amcol(i);
   RecSum1 = RecSum1 + (data1(find_index(center_pos1 + TriM1 )) - data1(find_index( center_pos1 - TriM1 - 1 )))/NRec_1 ;
   
   %% 第二层
   input_index = center_pos2 + TriM2;
   data2( find_index(input_index)) = RecSum1;
   RecSum2 = RecSum2 + (data2(find_index(center_pos2 + TriM2 )) - data2(find_index( center_pos2 - TriM2 - 1 )))/NRec_2;
   
   %% 第三层
   input_index = center_pos3 + TriM3;
   data3( find_index(input_index)) = RecSum2;
   RecSum3 = RecSum3 + (data3(find_index(center_pos3 + TriM3 )) - data3(find_index( center_pos3 - TriM3 - 1 )))/NRec_3;
   
   %% second 第一层
   input_index = Seccenter_pos1 + SecM1;
   Secdata1( find_index(input_index) ) = data1(find_index( center_pos3 + SecM1 + SecM2 ));
   SecRecSum1 = SecRecSum1 + (  Secdata1(find_index(Seccenter_pos1 + SecM1 )) - Secdata1(find_index( Seccenter_pos1 - SecM1 - 1 )) )/SecN1;
   
   %% second 第二层
   input_index = Seccenter_pos2 + SecM2;
   Secdata2( find_index(input_index)) = SecRecSum1;
   SecRecSum2 = SecRecSum2 + ( Secdata2(find_index(Seccenter_pos2 + SecM2 )) - Secdata2(find_index( Seccenter_pos2 - SecM2 - 1 )) )/SecN2;
   %% 计算结果
   yout(i) = data1(center_pos3) - 1.5 * RecSum3 + 0.5 * SecRecSum2;%%center_pos3象征着延时的点数
   
   %% 更新
   center_pos1 = center_pos1+1;
   center_pos2 = center_pos2+1;
   center_pos3 = center_pos3+1;
   center_pos1 = find_index(center_pos1);
   center_pos2 = find_index(center_pos2);
   center_pos3 = find_index(center_pos3);
   
   Seccenter_pos1 = Seccenter_pos1 + 1;
   Seccenter_pos2 = Seccenter_pos2 + 1;
   Seccenter_pos1 = find_index(Seccenter_pos1);
   Seccenter_pos2 = find_index(Seccenter_pos2);
end

%% 下面说明了延时
figure1 = figure('Color',[1 1 1]);plot(yout(TriM1+TriM2+TriM3+1:end));
hold on;plot(amcol);
function out = find_index(in)

%% 用在滤波器的队列中
Num = 2048;%%这个队列不会改变延时的特性
if mod(in,Num) == 0
   out = Num;
else
   out = mod(in,Num);
end
end

模拟58.9m的波长信号输入滤波器，发现与幅度谱表现一致，说明滤波器的正确性。