原创 MATLAB信号处理仿真-FIR滤波器之基于样点处理的模型

 2012-10-9 22:47  4207 17 19 分类: 处理器与DSP 文集: MATLAB信号处理仿真

FIR滤波器是数字信号处理中频频出镜的算法之一，之所以如此露脸，是因为具有某些优良品质。
品质一、（广义）线性相位，就是说信号的任何频率成分经过这个滤波器，延时都一样（群延时恒定），这个性质对于通信系统比较重要，对于人耳来说不是很重要，因为人耳对相位不敏感。
品质二、可以较为精确的设定幅频响应的形状，相比之下，IIR或是模拟滤波器都是用零极点来控制频响，比较不灵活，fir可以设定一个幅频曲线，然后用算法生成抽头系数，当然，越精确，需要的抽头系数就越多，运算量越大。
品质三、相对于带有反馈回路的IIR而言，FIR不会自激振荡，因为它没有反馈的数据回路，IIR则不是，数据回路意味着存在极点，它对系数的量化噪声非常敏感，量化误差造成的极点位置偏移很容易使得滤波器出现自激振荡，而FIR仅有零点，系数和数据的量化误差仅会让输出数据的噪声增加而不会自激振荡。

天下没有免费的午餐，这一切的优良品质都是建立在一个事实的基础上，就是FIR比IIR的运算量庞大的多，同样的频响，fir的运算量是IIR的十几倍很常见。不过，话说回来，各种FPGA、DSP处理器拼命的往里面塞乘法器、RAM、浮点单元，还不是为了做乘加么，如果不是为了发论文申专利，没必要死命的省运算量，算法和电路结构的优化是没有尽头的，差不多就好，有时候时间成本更重要。

另一方面，FIR代表了最经典的一种信号处理的算法过程，向量的点积，就两个向量的对应分量进行两两相乘，然后把所有的乘积加在一起。这个运算过程是信号处理里面最常见的，比如在DFT中会出现，在FFT中有它的变种会出现。说的务虚一点，这就是计算数据向量在基函数上面投影的过程。

在信号处理系统中，FIR滤波器的实现方法有两种形式，软件算法或是数字电路，在通信系统的数字中频部分，通常用数字电路来实现，因为这样能够在单位时钟周期的时间精度上保证数据处理的实时性。在这样的系统中，滤波器的工作时序很简单，每N个时钟周期里，给电路输入一个样点，电路输出一个样点，同时滤波器更新自己缓存。

上面说到软件算法和数字电路，通常信号处理的软件使用C语言编写，数字电路用verilog或是vhdl来设计，但是在硬件的环境中，调试是一件比较麻烦的事情，因为捕获和观察数据的手段比较有限，所以，通常先用matlab来做个仿真，然后对着这个仿真来编写C程序或是verilog代码，这样硬件上的信号处理开发就变成了“大家来找碴”的游戏，即：对着matlab仿真里面的正确数据来判断硬件环境里面的数据是否正确。

这里首先给出一个基于样点的FIR算法的仿真代码，这个仿真分为两部分，首先是用matlab的滤波器设计工具fir2设计出滤波器的系数，然后用freqz绘出幅频和相频响应曲线。需要注意的是freqz函数绘出的曲线是一种静态分析的方法，即：根据系数计算出零极点，然后对单位圆进行角频率扫描，计算出幅频相频曲线。当然，上述过程仅存在于概念中，实际干活时，是用FFT算的，还记得么？“DFT是对Z变换单位圆上的等分圆周采样”噢。另一方面，我们可以用动态的方法进行测试，即：生成某种形式的信号，然后用得到的滤波器系数和信号进行运算，然后再观察运算结果的频谱。实际上，对于复杂的通信系统，我们都是用动态仿真的方法，因为算法已经复杂到无法静态分析的程度。

在动态仿真的过程中，为了实现这个滤波器的算法模型，我们可以有两种做法，一是用matlab提供的函数，此处用了conv函数，另外还有诸如filter函数，也可以进行滤波运算。使用matlab函数的优势在于，只要正确的使用，结果通常是正确的，这种做法可以让我们快速的对算法进行理论评估，比如，为了达到要求，在某个信号节点上到底应该用31阶还是63阶的FIR滤波器。
另一种方法，我们可以自己编写算法的计算过程函数，这种方法好处是算法内部细节的高度可见性，尤其是在设计数字电路的时候，当我们面对着从某个电路内部节点捕获出来的几千个数据样点波形时，如果没有一个实现计算出来的可以比对的正确结果，我们会非常的头疼。

此处提供的my_fir()函数就是一个自己编写的滤波器算法过程，该函数对于根据的数据和FIR滤波器系数进行计算，生成滤波器输出的结果，函数内部是一个两层循环的结构，内层循环用于把滤波器的数据向量和系数向量进行点积运算，外层循环用于逐次向内层循环更新滤波器的缓存数据。同时该代码中还展示了一个称为“录像带”的小技巧，对于内层循环中的滤波器累加变量，该变量的值被刷新后，历史值就被抹去，如果调试过程中需要看到该值的变化过程，则可以建立一个二维数据来保存所有中间结果。当然，由于这种“录像带”需要消耗大量的内存资源，所以，代码中建立了
一个控制开关变量，用于开启这个“录像带”的内存申请和功能启用。代码中还展示了绘制一个3d的图形，这样做仅仅是为了展示，并不是有什么实际的意义。

使用freqz得到的滤波器幅频-相频曲线

输入的多音正弦信号频谱

matlab滤波后的信号频谱

my_fir滤波后的信号频谱

用3D形式绘制的中间结果曲面图形（仅供演示用）

仿真代码（生成信号和绘制频谱的代码见于“MATLAB信号处理仿真-使用函数封装，多音正弦和谱分析”一文）

%///////////////////////////////////////////////////////////
% test1.m Simulation of sample based FIR model
%///////////////////////////////////////////////////////////
close all;
clear;
clc;

% specify input multi-tone sine component frequency
sin_freq    = [1:8]*30E3;

data_len    = 2048 ;% signal data length
fs          = 512E3 ;% sample rate
quant_bits = 12    ;% signal quant bits
kaiser_beta = 8     ;% beta of kaiser win

% design the FIR coefficents use matlab toolbox
n_coef      = 32;    % number of coefficents
f = [0, 0.2, 0.2, 1];
m = [1, 1 , 0 , 0];
coeff = fir2(n_coef-1, f, m);
% draw the filter reponse curve
freqz(coeff);

filter_in = gen_quant_multi_sin(fs, sin_freq, data_len, quant_bits);
data_conv = conv(coeff, filter_in);
filter_out = data_conv(1:data_len);
my_out     = my_fir(filter_in, coeff);

kaiser_win_spectrum_plot(fs, filter_in, kaiser_beta);
ylim([-160,10])     ;   % set y-axis range
title('Filter In, Normalized Spectrum', 'fontsize', 14);

kaiser_win_spectrum_plot(fs, filter_out, kaiser_beta);
ylim([-160,10])     ;   % set y-axis range
title('Matlab Out, Normalized Spectrum','fontsize', 14);

kaiser_win_spectrum_plot(fs, my_out, kaiser_beta);
ylim([-160,10])     ;   % set y-axis range
title('My Filter Out, Normalized Spectrum','fontsize', 14);

%///////////////////////////////////////////////////////////
% file: my_fir.m
% perform FIR fitering process to input data with
% given coefficents.
% the input data and output are same length
%///////////////////////////////////////////////////////////
function out = my_fir(in, coeff);
en_debug   = 1                 ; % debug switch
data_len   = length(in)        ;
n_coeff    = length(coeff)     ;
out        = zeros(data_len, 1);
sample_buf = zeros(n_coeff,1) ;

% combine the filter buffer data and the input data
total_buf             = zeros(n_coeff + data_len, 1);
total_buf(1:n_coeff) = sample_buf                  ;
total_buf(n_coeff+1:n_coeff + data_len) = in        ;

out = zeros(data_len,1);
if(en_debug) % create the debug mem area
    sum_debug = zeros(data_len,n_coeff);
end
% loop by each input data sample
for(idx_sample = 1:data_len)
    % get current loop fitler data in buffer
    sample_buf = total_buf(idx_sample:idx_sample+n_coeff-1);
    sum_val = 0;
    % loop by each filter mult
    for(idx_coeff = 1:n_coeff)
      % multiply and accumulate
      mult = sample_buf(n_coeff - idx_coeff + 1) * coeff(idx_coeff);
      sum_val = sum_val + mult;
      if(en_debug) % record the change of sum value
        sum_debug(idx_sample, idx_coeff) = sum_val;
      end
    end
    % save to output buffer
    out(idx_sample) = sum_val;
end
if(en_debug) %3d-plot debug value
    figure; surf(sum_debug); colormap hot;
end
end