2014年6月3日 星期二

FIR Filter Design in Verilog

源自於
http://www.vlsitechnology.com/fir-filter-design-in-verilog/
http://vaplab.ee.ncu.edu.tw/chinese/pcchang/course2005b/comsp/audio/Algorithm.htm

濾波器又分為兩種不同的系統:有限脈衝響應系統(Finite Impulse Response Filter)與無限脈衝響應系統(Infinite Impulse Response Filter),簡稱FIR系統與IIR系統。
 
Ø       FIR Filter
FIR 濾波器的響應會在有限的取樣時間內衰減,對於有限的輸入,其輸出響應必為有限,因此可視為是輸入值的加權平均,故又可稱為移動平均(Moving Average)濾波器,簡稱MA 濾波器。
FIR濾波器的脈衝響應係數h(n)是一個有限的項,而其運算量跟濾波器階數有很大的相關性,如果FIR濾波器的階數變大,運算量也會跟著變大,因運算量與硬體相關性過高,因此一般我們不採用此系統。


Ø       IIR Filter
IIR濾波器的脈衝響應h(n)為無限長度的項,但無限大的階數在實現上是一件不可能的事,因此我們採用遞迴的方式來處理訊號,在處理未知系統時,可因為其無限的特性達成使用較少階數濾波器的效果。又因為加入回授,所以輸出訊號將不僅跟輸入訊號有關,也跟上一筆輸出訊號有關連。
IIR 濾波器的輸出是由過去時間的輸出訊號以遞迴的方式而得,故也稱自動回歸(Autoregressive)濾波器,簡稱AR 濾波器。一般而言,一個IIR 濾波器,由於包含AR MA兩部分,故此類IIR 濾波器又稱為自動回歸移動平均(ARMA)濾波器。

FIR Filter Design in Verilog





module fir(data_out,data_in,clock,reset);
parameter order=8;
parameter word_size_in=8;
parameter word_size_out=2*word_size_in+2;
parameter b0=8'd1;
parameter b1=8'd1;
parameter b2=8'd1;
parameter b3=8'd1;
parameter b4=8'd1;
parameter b5=8'd1;
parameter b6=8'd1;
parameter b7=8'd1;
parameter b8=8'd1;
output[word_size_out-1:0] data_out;
input[word_size_in-1:0] data_in;
input clock,reset;
reg[word_size_in-1:0] samples[0:order];
integer k;
assign data_out=b0*samples[0]+b1*samples[1]+b2*samples[2]+b3*samples[3]+b4*samples[4]+b5*samples[5]+b6*samples[6]+b7*samples[7]+b8*samples[8];
always@(posedge clock)
if(reset==0)
begin
for(k=0;k<=order;k=k+1)
samples[k]<=0;
end
else
begin
samples[0]<=data_in;
for(k=1;k<=order;k=k+1)
samples[k]<=samples[k-1];
end
endmodule




Testbench


module fir_tb;
    // Inputs
    reg [7:0] data_in;
    reg clock;
    reg reset;
    // Outputs
    wire [17:0] data_out;
    // Instantiate the Unit Under Test (UUT)
    fir uut (
        .data_out(data_out),
        .data_in(data_in),
        .clock(clock),
        .reset(reset)
    );
    initial
    clock=1'b0;
    always #50 clock=~clock;
    initial begin
        // Initialize Inputs
        data_in = 8'b00000000;
        reset = 0;
        // Wait 100 ns for global reset to finish
        #100;
        data_in = 8'b00000001;
        reset = 1;
        #100;
        data_in = 8'b00000010;
        reset = 1;
        #100;
        data_in = 8'b00000011;
        reset = 1;
        #100;
        data_in = 8'b00000100;
        reset = 1;
        #100;
        data_in = 8'b00000101;
        reset = 1;
        #100;
      data_in = 8'b00000110;
        reset = 1;
        #100;
        data_in = 8'b00000111;
        reset = 1;
        #100;
        data_in = 8'b00000111;
        reset = 1;
        #100;
        data_in = 8'b11111111;
        reset = 1;
        #100;
        data_in = 8'b00000010;
        reset = 1;
        #100;
        // Add stimulus here
    end
endmodule







FIR filter implementation

FIR (finite impulse response) filter is a digital filter without feedback. Absence of feedback makes impulse filter response to be finite. Digital filters can be with infinite impulse response (IIR filters).
The advantages of FIR filters compared to IIR filters:
1. FIR filter can be designed with linear phase.
2. FIR filter is simple to implement.

 3. FIR filter can be easily implemented on finite-precision arithmetic (a lot of microcontrollers can operate with 16-bit words, but for IIR filter correct working, in some cases, you need 32 bits to store "Y" coefficients. So there is much more problems with IIR filter implementation on 16-bit MCU, than with FIR filter implementation).

The disadvantages of FIR filters compared to IIR filters:
1. FIR sometimes requires more memory and calculations to achieve a desired response characteristic.
2. Some responses are not practical to implement.

The equation of the FIR filter is as follows:
Input x(n) and output y(n) signals are related via the convolution operation. h(k) - filter coefficients. N is a filter length.
Block diagram of the FIR filter:
Coefficients h(n):
Input data x(n):
Let N = 3 then:
Coefficients h(n):
 
Input data x(n): 

The equation in this case is as follows: 

A detailed description of the entry and modification of data in the filter is given below:


Consider the filter implementation in C language. The input signal is a unit impulse at zero timeThus, finding the Fourier transform of the output signal we obtain the amplitude-frequency characteristic of the filter. The length of the input signal is 1000 points. After filtering the signal, the result is stored in the file Y.dat using WinAPI functions.
int main()
{

    float X[1000] = {0};
    X[0] = 1;

    float Y[1000] = {0};


    float yn = 0;
    const int N = 3;
    float h[N] = {0.32590173795979338000,
                  0.34819652408041324000,
                  0.32590173795979338000};
    float x[N] = {0,0,0};


    for(int i=0; i<1000; i++)
    {
                                                 //  Alternative implementation
        for(int k=0; k < N-1; k++)               //  for(int k=N-1; k>0; k--)
        {                                        //  {                      
          x[N-k-1] = x[N-k-2];//shift the data   //    x[k] = x[k-1];
        }                                        //  }              

              
        x[0] = X[i]; // move input sample to buffer
        yn = 0; // clear output sample

        for(int k=0; k < N; k++)
        {
            yn += h[k]*x[k]; // multiply data on coefficients with accumulation
        }

        Y[i] = yn;
    }

    HANDLE hFile = CreateFile("Y.dat", GENERIC_READ | GENERIC_WRITE, 0, NULL, 
                               CREATE_ALWAYS, 0, NULL);

    DWORD dwCount;
    WriteFile(hFile,Y,sizeof(Y),&dwCount,NULL);
    CloseHandle(hFile);

    return 0;
}

Filter response can be plotted in MATLAB. To do this, open the file Y.dat, read the response of the filter and apply Fourier transform.
clc;
close all;
clear all;

fid = fopen('Y.dat','r'); %open file
F = fread(fid,'float');   %read file
fclose(fid);              %close file 

N = length(F);            %get signal length (number of samples)
Fs = 500;                 %set sample frequency
dt = 1/Fs;                %get sample period in time domain
df = Fs/N;                %get sample period in frequency domain
t = (1:N)*dt;             %get time vector 
f = (1:N)*df;             %get frequency vector 

plot(t,F);                %plot signal in time 

figure;
S = abs(fft(F));          %get frequency response of filter 
plot(f(1:N/2),20*log10(S(1:N/2)/max(S(1:N/2))));  %plot frequency response
xlabel('Hz');
ylabel('S, dB');

Frequency response of the filter:

Sources:
1. Emmanuel C. Ifeachor and Barrie W. Jervis Digital Signal Processing: A Practical Approach2nd Edition.Hardcover. from English. - Moscow: Publishing House "Williams", 2008. - 992 pp.

沒有留言:

張貼留言