Sunday, November 18, 2012

Filter design: From Matlab to C code

Note 1:  For this, I am using not only Matlab but the Signal Processing Toolbox. Got to study the other tools that Matlab has to deal with filters...
Note 2: I should write something like this but with some free tool out there... Anyhow, for the moment just use Matlab, as I have it available in my company...

Quick reading, it looks like there are couple of methods:
  1. A=ellip(specs) --> H=dfilt.df1(A) where A are the numerator and denominator coeff. This is kind of going analog to then move to digital.
  2. D=fdesign.bandpass(specs) --> H=design(D, ellip)
H is the object describing the digital filter (not exactly the coefficients, see more below). We will look in detail and extend the method 2 following Getting started with discrete time filters (dfilt) or fdesign objects.
1/ First got to get the object describing the specification of the filter we want. So, we use fdesign. Just type  "fdesign.bandpass" for a very simple syntax, or go to fdesign for more info.
In our case, we want to create a bandpass filter around 2.5KHz. So, we will do:

>> D=fdesign.bandpass(2000,2200,2800,3000,40,1,40,44100)
And we get:
D =
               Response: 'Bandpass'                     
          Specification: 'Fst1,Fp1,Fp2,Fst2,Ast1,Ap,Ast2'
            Description: {7x1 cell}                     
    NormalizedFrequency: false                          
                     Fs: 44100                          
                 Fstop1: 2000                           
                 Fpass1: 2200                           
                 Fpass2: 2800                           
                 Fstop2: 3000                           
                 Astop1: 40                             
                  Apass: 1                              
                 Astop2: 40

2/ To get the real design, we use:
>> h=design(D,'ellip')
And get:
h =
         FilterStructure: 'Direct-Form II, Second-Order Sections'
              Arithmetic: 'double'                              
               sosMatrix: [4x6 double]                          
             ScaleValues: [5x1 double]                          
     OptimizeScaleValues: true                                  
        PersistentMemory: false
For help, type designmethods(h) to get a list of filters.

3/ At this moment, we may actually like a different structure. Notice it used the default one (direct form II)


but we may want the more standard one. I am going to use direct form I as it is easier for me to identify what output number below matches what coefficient (notice that I am learning while doing this), so, it'll be easier to program:

Note: notice the error on the - sign in the added of b(nb) path.

So, we do:
>> h1  = convert(h,'df1')

h1 =

     FilterStructure: 'Direct-Form I'
          Arithmetic: 'double'      
           Numerator: [1x9 double]  
         Denominator: [1x9 double]  
    PersistentMemory: false

To see what we get, we can do:
>> info(h1)
Discrete-Time IIR Filter (real)   
-------------------------------   
Filter Structure    : Direct-Form I
Numerator Length    : 9           
Denominator Length  : 9           
Stable              : Yes         
Linear Phase        : No

4/ At this point, if we want to see what we got, it is as simple as using the Filter Visualization Tool:
hfvt=fvtool(h1);


5/ To get the coefficients, use:
num = get(h1,'Numerator');
den = get(h1,'Denominator');

In our case:
num =     0.0099   -0.0729    0.2418   -0.4681    0.5789   -0.4681  0.2418   -0.0729    0.0099
den =     1.0000   -7.4222   24.5776  -47.3746   58.1077  -46.4331  23.6104   -6.9884    0.9229

Where numerator and denominator are (when a0=1, like in our case):


In this kind of structure

Notice that a0 in Matlab is den(1) as the index starts at 1. So, we create in Matlab a piece of code that will do this (equivalent to running y=filter(h1,rx)):

rx2=[zeros(length(num)-1,1);rx];
y=zeros(length(rx)+length(num)-1,1);
for R=length(num):length(rx2)
    data=0;
    for N=1:length(num)
        data=data+rx2(R-N+1)*num(N);
    end
    for D=1:(length(den)-1)
        data=data-y(R-D)*den(D+1);
    end
    y(R)=data/den(1);
end

Translating this into C is straightforward. Notice that the initialization is important. I.e., getting both, input and output past values equal to zero. We can do that adding those in front of our sequences/arrays, or, just make zero the first set of values on those arrays. Of course, this is just for test and assuming you don't care about that data being lost... Either way, you get the point :)
 
The plot on the left is the result using "filter" while the one on the right is the one from the "C" code.


Other useful information:

For an automatic way to do it: Forum describing how to do it
From Simulink to C
From Matlab to HDL
Real Time Workshop Embedded Coder

Matlab filter design
Designing Low Pass FIR filters
FDATool

To use your filter (DUH!): y=filter(num,den,x);
To see the impulse response of your filter use [h,t] = impz(num,den) and plot(t,h)
To check if it is stable: isstable(h) - Note: I didn't get it to work where I can check it with filter running in a given precision

Nice write up on filter theory: http://www.mikroe.com/chapters/view/73/chapter-3-iir-filters/
And another one...
And one more!

No comments:

Post a Comment