CLA IIR Filter implementation

In C2000 DSP there is no floating point IIR library support. But you can find the fixed point version of IIR implementation in controlSUITE. It is under directory FixedPointLib/v101/source/iir5biq16 or iir5biq32.  Not pretty sure how many orders IIR filter are you going to implement and how do you arrange your coefficient buffer and delay buffer. What I recommend you to do is starting debugging from 2 order, as long as 2 order is done, then you can use biquad method to increase the order.

You can also check the associated document of that two fixed point functions (pg. 39 and pg. 44). The coefficient arrangement of IIR filter is different from FIR filter if you would like to get an efficient biquad solution.

Hi!

Sure, here it is: (in another task I initialize X and Y arrays with 0). I've also tried with MMOVD32 like for the feedforward part.

Many thanks,

Monica

_Cla1Task1:
        .if CLA_DEBUG == 1
        MDEBUGSTOP
       .endif
   
;//==============================================
;// CLA Task 1
;//
; float32 B [FILTER_LEN] = {0.0625L, 0.25L, 0.0375L, 0.25L, 0.0625L};   
; float32 A [FILTER_LEN] = { 1.0L, 0.02L, 0.25L, 0.25L, 0.02L};   

; y(k) = b(0)*x(k) + b(1)*x(k-1) + ... + b(n-1)*x(k-n+1)-a(1)*y(k-1)...-a(n-1)*y(k-n+1)
;  Y0 = B0 * X0 + B1 * X1 + B2 * X2 + B3 * X3 + B4 * X4  - A1 * Y1 - A2 * Y2 - A3 * Y3 - A4 * Y4
;
; Coefficients A[0, 1, 2, 3, 4] ; B[0, 1, 2, 3, 4]
; Data         X[0, 1, 2, 3, 4] (Delay Line - X[0] is newest value)
;              Y[0, 1, 2, 3, 4] (Delay Line - Y[0] is newest filtered value - shifted on Y1 at the end of the TASK)
;
_X4  .set _X+8 
_X3  .set _X+6 
_X2  .set _X+4 
_X1  .set _X+2 
_X0  .set _X+0 

_Y4  .set _Y+8 
_Y3  .set _Y+6 
_Y2  .set _Y+4 
_Y1  .set _Y+2 
_Y0  .set _Y+0 

_A4  .set _A+8 
_A3  .set _A+6 
_A2  .set _A+4 
_A1  .set _A+2 
_A0  .set _A+0

_B4  .set _B+8 
_B3  .set _B+6 
_B2  .set _B+4 
_B1  .set _B+2 
_B0  .set _B+0

    MMOV32     MR0,@_X4                      ;1 Load MR0 with X4
    MMOV32     MR1,@_B4                      ;2 Load MR1 with B4
    MNOP                                     ;3 Wait till I8 to read result
    MNOP                                     ;4 Wait till I8 to read result
    MNOP                                     ;5 Wait till I8 to read result
    MNOP                                     ;6 Wait till I8 to read result
    MNOP                                     ;7 Wait till I8 to read result
    MUI16TOF32 MR2,  @_AdcResult.ADCRESULT0  ;8 Read ADCRESULT0 and convert to float


    MMPYF32    MR2, MR1, MR0                 ; MR2 (Y) = MR1 (B4) * MR0 (X4)
 || MMOV32     @_X0, MR2

    MMOVD32    MR0,@_X3                      ; Load MR0 with X3, Load X4 with X3
    MMOV32     MR1,@_B3                      ; Load MR1 with B3
                                         
    MMPYF32    MR3, MR1, MR0                 ; MR3 (Y) = MR1 (B3) * MR0 (X3)
 || MMOV32     MR1,@_B2                      ; Load MR1 with B2
    MMOVD32    MR0,@_X2                      ; Load MR0 with X2, Load X3 with X2

    MMACF32    MR3, MR2, MR2, MR1, MR0       ; MR3 = B3*X3 + B4*X4
 || MMOV32     MR1,@_B1                      ; MR2 = MR1 (B2) * MR0 (X2)
    MMOVD32    MR0,@_X1                      ; Load MR0 with X1, Load X2 with X1

    MMACF32    MR3, MR2, MR2, MR1, MR0       ; MR3 = B2*X2 + (B3*X3 + B4*X4)
 || MMOV32     MR1,@_B0                      ; MR2 = MR1 (B1) * MR0 (X1)
    MMOVD32    MR0,@_X0                      ; Load MR0 with X0, Load X1 with X0

    MMACF32    MR3, MR2, MR2, MR1, MR0       ; MR3 = B1*X1 + (B2*X2 +B3*X3 + B4*X4)
 || MMOV32     MR1,@_B0                      ; MR2 = MR1 (B0) * MR0 (X0)

    MADDF32    MR3, MR3, MR2                 ; MR3 = B0*X0 + (B1*X1 + B2*X2 +B3*X3 + B4*X4)

    MMOV32    @_R1, MR3                      ; Output R1=B0 * X0 + B1 * X1 + B2 * X2 + B3 * X3 + B4 * X4
   
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;   
    MMOV32     MR0,@_Y4                      ;1 Load MR0 with Y4
    MMOV32     MR1,@_A4                      ;2 Load MR1 with A4
   
    MMPYF32    MR2, MR1, MR0                 ; MR2 (Y) = MR1 (A4) * MR0 (Y4)

    MMOV32     MR0,@_Y3                      ; Load MR0 with Y3
    MMOV32     MR1,@_A3                      ; Load MR1 with A3
                                         
    MMPYF32    MR3, MR1, MR0                 ; MR3 (Y) = MR1 (A3) * MR0 (Y3)
 || MMOV32     MR1,@_A2                      ; Load MR1 with A2
    MMOV32     MR0,@_Y2                      ; Load MR0 with Y2

    MMACF32    MR3, MR2, MR2, MR1, MR0       ; MR3 = A3*Y3 + A4*Y4
 || MMOV32     MR1,@_A1                      ; MR2 = MR1 (A2) * MR0 (Y2)
    MMOV32     MR0,@_Y1                      ; Load MR0 with Y1

    MMACF32    MR3, MR2, MR2, MR1, MR0       ; MR3 = A2*Y2 + (A3*Y3 + A4*Y4)
 || MMOV32     MR1,@_A0                      ; MR2 = MR1 (A1) * MR0 (Y1) ; dummy load
   
   
    MADDF32    MR3, MR3, MR2                 ; MR3 = A1*Y1 + (A2*Y2 +A3*Y3 + A4*Y4)
    MMOV32     @_R2, MR3                     ; Output R2= A1 * Y1 + A2 * Y2 + A3 * Y3 + A4 * Y4
   
    MMOV32     MR0, @_R1                  
    MMOV32     MR1, @_R2
    MSUBF32    MR0,MR0,MR1
    MMOV32     @_Y0, MR0
   
    MMOV32     MR0,@_Y3                      ; shift the filtered values
    MMOV32     @_Y4,MR0                     
    MMOV32     MR0,@_Y2                     
    MMOV32     @_Y3,MR0                     
    MMOV32     MR0,@_Y1                     
    MMOV32     @_Y2,MR0                     
    MMOV32     MR0,@_Y0
     MMOV32     @_Y1,MR0                  

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;   
   

   
   
   
    MSTOP                                    ; End task  
  
  
_Cla1T1End:

Monica,

I cant exactly figure out whats going wrong in your code in the feedback part but i did write up my own version of the IIR and crosschecked it with a MATLAB script and it works fine. It might be helpful to you

ASSUMPTIONS:

1. I negated all the A coefficients so i could continue using the MMACF32 without having to resort to an MSUBF32

2. For simulation purposes I fed the task a wave of 128 data points, 0.5(sin(2*pi*0.1) + sin(2*pi*0.46)) in the variable Temp

    MMOV32     MR1,@_X4                      ;1 Load MR1 with X4
    MMOV32     MR2,@_B4                      ;2 Load MR2 with B4
    MNOP                                     ;3 Wait till I8 to read result
    MNOP                                     ;4 Wait till I8 to read result
    MNOP                                     ;5 Wait till I8 to read result
    MNOP                                     ;6 Wait till I8 to read result
    MNOP                                     ;7 Wait till I8 to read result
    MMOV32    MR0,  @_Temp                   ;8 Read a temp value(simulating an
                                             ;ADC input) into MR0
    MMPYF32      MR3,MR1,MR2 ||                 ;Y(MR3) = X4(MR1)*B4(MR2)
    MMOV32      @_X0,MR0                         ;X0 = latest ADC input
    MMOV32      MR0,@_B3                         ;MR0 = B3
    MMOVD32      MR1,@_X3                         ;MR1 = X3, X4 = X3
   
    MMPYF32      MR2,MR0,MR1 ||                ;Y`(MR2) = B3(MR0) * X3(MR1)
    MMOV32      MR0,@_B2                        ;MR0 = BR2
    MMOVD32      MR1,@_X2                        ;MR1 = X2, X3 = X2
   
    MMACF32      MR3,MR2,MR2,MR0,MR1 ||        ;Y(MR3) = Y+Y`, Y"(MR2)=B2(MR0)*X2(MR1)
    MMOV32      MR0,@_B1                        ;MR0 = B1
    MMOVD32      MR1,@_X1                        ;MR1 = X1, X2 = X1
   
    MMACF32      MR3,MR2,MR2,MR0,MR1 ||        ;Y(MR3) = Y+Y", Y`(MR2)=B1(MR0)*X1(MR1)
    MMOV32      MR0,@_B0                        ;MR0 = B0
    MMOVD32      MR1,@_X0                        ;MR1 = X0, X1 = X0
   
    MMACF32      MR3,MR2,MR2,MR0,MR1 ||        ;Y(MR3) = Y+Y', Y"(MR2)=B0(MR0)*X0(MR1)
    MMOV32      MR0,@_A4                        ;MR0 = A4
    MMOV32      MR1,@_Y4                        ;MR1 = Y4
   
    MMACF32      MR3,MR2,MR2,MR0,MR1 ||        ;Y(MR3) = Y+Y", Y`(MR2)=A4(MR0)*Y4(MR1)
    MMOV32      MR0,@_A3                        ;MR0 = A3
    MMOVD32      MR1,@_Y3                        ;MR1 = Y3, Y4 = Y3
   
    MMACF32      MR3,MR2,MR2,MR0,MR1 ||        ;Y(MR3) = Y+Y`, Y"(MR2)=A3(MR0)*Y3(MR1)
    MMOV32      MR0,@_A2                        ;MR0 = A2
    MMOVD32      MR1,@_Y2                        ;MR1 = Y2, Y3 = Y2   

    MMACF32      MR3,MR2,MR2,MR0,MR1 ||        ;Y(MR3) = Y+Y", Y`(MR2)=A2(MR0)*Y2(MR1)
    MMOV32      MR0,@_A1                        ;MR0 = A1
    MMOVD32      MR1,@_Y1                        ;MR1 = Y1, Y2 = Y1
   
    MMACF32      MR3,MR2,MR2,MR0,MR1 ||        ;Y(MR3) = Y+Y`, Y"(MR2)=A1(MR0)*Y1(MR1)
    MMOV32      MR0,@_A0                        ;MR0 = A0 ; dummy read
    MMOVD32      MR1,@_Y0                        ;MR1 = Y0, Y1 = Y0; dummy read

    MMOV32      @_Y0,MR3                        ;Y0 = Y(MR3)
                   
    MSTOP                                    ; End task

The MATLAB script: I used the fdatool to get the coefficients for the C code

clc
clear all
close all

%%IIR Implementation on the CLA
SAMPLING_FREQ = 100e3;
FREQ1 = 10e3;
FREQ2 = 46e3;
n = linspace(0,128,128);
fd1= FREQ1/SAMPLING_FREQ;
fd2 = FREQ2/SAMPLING_FREQ;
s = 0.5*sin(2*pi*n*fd1)+ 0.5*sin(2*pi*n*fd2);;
figure,plot(s);
axis([0 128 -1.5 1.5]);
grid on;

%% Generating the filter coeffs
n = 4;
Wn = 0.5;
[b,a] = butter(n,Wn);

%% Filtering the input
y = filter(b,a,s);
figure,plot(y);
grid on;

Hi Vishal,

Thanks for the code, it's for sure much compact. I've tested it with a rectangular signal, and it works like mine.

I've also checked it with a MatLab script, and at the input of the ADC I use a PWM generated on the board 600Hz (9KHz sampling rate). On the 9kHz interrupt I do the IIR filtering in CLA. 

The MatLab script is:

clc;
clear all;
close all;

f=600;
T=1/f;
step=1/9000;
t=0:step:10*T;

tri=2.*3.23*asin(sin(2*pi*f.*t+pi))./pi;
threshold=0;


for i=1:length(t)
    if tri(i)>=threshold
        in1(i) = 3.23;
    else
        in1(i) = 0;
    end
end

figure;
plot(t,in1,'b','linewidth',2);hold on;

 
% b0=0.0625;
% b1=0.25;
% b2=0.0375;
% b3=0.25;
% b4=0.0625;

% a0 =  1.0000 
% a1=-0.0000
% a2=0.4860  
% a3=-0.0000  
% a4=0.0177

a0=0.5;
a1=0.5;
a2=0.25;
a3=0.25;
a4=0.5;

b0 = 0.0940; 
b1=0.3759;   
b2=0.5639;  
b3=0.3759; 
b4=0.0940;

in1_filt=zeros(1,length(t));

for i=2:length(t)
    if i>=5
        in1_filt(i)=b0*in1(i)+b1*in1(i-1)+b2*in1(i-2)+b3*in1(i-3)+b4*in1(i-4)- a1*in1_filt(i-1)-a2*in1_filt(i-2)-a3*in1_filt(i-3)-a4*in1_filt(i-4);

%  if in the line above I change:

%  -a1*in1_filt(i-2)-a2*in1_filt(i-3)-a3*in1_filt(i-4)-a4*in1_filt(i-5);

% and if>=6

% the result match perfectly (shape and levels)

    else
    in1_filt(i)=0;      
    end
end

figure;
plot(t,in1_filt,'*b-','linewidth',2);hold on;

Many thanks,

Monica

Hi!

I've found it  in the meantime. First time the task is called, it shifts twice the values in Y. So, the Y buffer contains 5 zeros at the beginning of the task, and at the end of the task the new filtered value is both on position 0 and 1, so instead of having 4 zero and the new value, I had only 3 zeros. I've changed the order of coefficients and it works great for any values of A[] and B[].

Here's the code in case somebody needs it:

_Cla1Task1:
        .if CLA_DEBUG == 1
        MDEBUGSTOP
       .endif
   
;//==============================================
;// CLA Task 1
;//
; float32 B [FILTER_LEN] = {B0, B1, B2, B3, B4};   
; float32 A [FILTER_LEN] = {A1, A2, A3, A4, 1/A0};   
;
;
;        (B0 * X0 + B1 * X1 + B2 * X2 + B3 * X3 + B4 * X4  - A1 * Y1 - A2 * Y2 - A3 * Y3 - A4 * Y4)
;   Y0 = -----------------------------------------------------------------------------------------------------------------
;                                                                                    A0
;
; Coefficients A[0, 1, 2, 3, 4] ; B[0, 1, 2, 3, 4]
; Data         X[0, 1, 2, 3, 4] (Delay Line - X[0] is newest value)
;              Y[0, 1, 2, 3, 4] (Delay Line - Y[0] is newest filtered value)
;
_X4  .set _X+8 
_X3  .set _X+6 
_X2  .set _X+4 
_X1  .set _X+2 
_X0  .set _X+0 

_Y4  .set _Y+8 
_Y3  .set _Y+6 
_Y2  .set _Y+4 
_Y1  .set _Y+2 
_Y0  .set _Y+0 

_A4  .set _A+8 
_A3  .set _A+6 
_A2  .set _A+4 
_A1  .set _A+2 
_A0  .set _A+0

_B4  .set _B+8 
_B3  .set _B+6 
_B2  .set _B+4 
_B1  .set _B+2 
_B0  .set _B+0

     MMOV32     MR1,@_X4                      ;1 Load MR1 with X4
     MMOV32     MR2,@_B4                      ;2 Load MR2 with B4
     MNOP                                     ;3 Wait till I8 to read result
     MNOP                                     ;4 Wait till I8 to read result
     MNOP                                     ;5 Wait till I8 to read result
     MNOP                                     ;6 Wait till I8 to read result
     MNOP                                     ;7 Wait till I8 to read result
     MUI16TOF32  MR0, @_AdcResult.ADCRESULT0  ;8 Read ADCRESULT0 and convert to float

    MMPYF32     MR3,MR1,MR2 ||                 ;Y(MR3) = X4(MR1)*B4(MR2)
    MMOV32      @_X0,MR0                       ;X0 = latest ADC input
    MMOV32      MR0,@_B3                       ;MR0 = B3
    MMOVD32     MR1,@_X3                       ;MR1 = X3, X4 = X3
  
    MMPYF32     MR2,MR0,MR1 ||                 ;Y`(MR2) = B3(MR0) * X3(MR1)
    MMOV32      MR0,@_B2                       ;MR0 = BR2
    MMOVD32     MR1,@_X2                       ;MR1 = X2, X3 = X2
  
    MMACF32     MR3,MR2,MR2,MR0,MR1 ||         ;Y(MR3) = Y+Y`, Y"(MR2)=B2(MR0)*X2(MR1)
    MMOV32      MR0,@_B1                       ;MR0 = B1
    MMOVD32     MR1,@_X1                       ;MR1 = X1, X2 = X1
  
    MMACF32     MR3,MR2,MR2,MR0,MR1 ||         ;Y(MR3) = Y+Y", Y`(MR2)=B1(MR0)*X1(MR1)
    MMOV32      MR0,@_B0                       ;MR0 = B0
    MMOVD32     MR1,@_X0                       ;MR1 = X0, X1 = X0
  
    MMACF32     MR3,MR2,MR2,MR0,MR1 ||         ;Y(MR3) = Y+Y', Y"(MR2)=B0(MR0)*X0(MR1)
    MMOV32      MR0,@_A3                       ;MR0 = A3
    MMOVD32     MR1,@_Y3                       ;MR1 = Y3, Y4 = Y3
  
    MMACF32     MR3,MR2,MR2,MR0,MR1 ||         ;Y(MR3) = Y+Y", Y'(MR2)=A3(MR0)*Y3(MR1)
    MMOV32      MR0,@_A2                       ;MR0 = A2
    MMOVD32     MR1,@_Y2                       ;MR1 = Y2, Y3 = Y2  

    MMACF32     MR3,MR2,MR2,MR0,MR1 ||         ;Y(MR3) = Y+Y', Y"(MR2)=A2(MR0)*Y2(MR1)
    MMOV32      MR0,@_A1                       ;MR0 = A1
    MMOVD32     MR1,@_Y1                       ;MR1 = Y1, Y2 = Y1
   
    MMACF32     MR3,MR2,MR2,MR0,MR1 ||         ;Y(MR3) = Y+Y", Y'(MR2)=A1(MR0)*Y1(MR1)
    MMOV32      MR0,@_A0                       ;MR0 = A0
    MMOVD32     MR1,@_Y0                       ;MR1 = Y0, Y1 = Y0
   
    MMACF32     MR3,MR2,MR2,MR0,MR1 ||         ;Y(MR3) = Y+Y', Y"(MR2)=A0(MR0)*Y0(MR1)
    MMOV32      MR0,@_A4                       ;MR0 = A4
   
    MADDF32     MR3, MR3, MR2                  ;Y(MR3) = Y+Y"
    MMPYF32     MR3,MR3,MR0                    ;Y(MR3) = Y(MR3) * A4(MR0)
   
    MMOV32      @_Y0,MR3                       ;Y0 = Y(MR3)
                
    MSTOP                                      ; End task
  
  
_Cla1T1End:

MatLab script to check it:

clc;
clear all;
close all;

f=600;
T=1/f;
step=1/9000;
t=0:step:10*T;

tri=2.*3.23*asin(sin(2*pi*f.*t+pi))./pi;
threshold=0;


for i=1:length(t)
    if tri(i)>=threshold
        in1(i) = 3.23;
    else
        in1(i) = 0;
    end
end

figure;
plot(t,in1,'b','linewidth',2);hold on;

 


% a0= 1.0000 
% a1=-0.0000
% a2= 0.4860  
% a3=-0.0000  
% a4= 0.0177
% b0=0.0940; 
% b1=0.3759;   
% b2=0.5639;  
% b3=0.3759; 
% b4=0.0940;

a0=0.5;
a1=0.05;
a2=0.25;
a3=0.25;
a4=0.05;

b0=0.0625;
b1=0.25;
b2=0.0375;
b3=0.25;
b4=0.0625;

in1_filt=zeros(1,5);

for i=1:length(t)
    if i>=5
        in1_filt(i)=(1/a0)*(b0*in1(i)+b1*in1(i-1)+b2*in1(i-2)+b3*in1(i-3)+b4*in1(i-4)-a1*in1_filt(i-1)-a2*in1_filt(i-2)-a3*in1_filt(i-3)-a4*in1_filt(i-4));
%     else
%     in1_filt(i)=0;      
    end
end

figure;
plot(t,in1_filt,'*b-','linewidth',2);hold on;

Thanks Vishal for the code, it's very nice written, looks profi compared to mine ;) I don't think I would have spend time on this, the moment it was working I wouldn't have touch it again :)

Monica

hmmm, could you post the MATLAB code. Id like to take a stab at it

Sure. Here it is. I'm testing with both rectangular and sinewave sampled at 9KHz. What I get in MatLab it's what I need. Can't obtan it in reality...

Thanks again,

Monica

I hope it's readable :)

clc;
clear all;
close all;

%% input signal - rectangular 600Hz - it is generated on board -> ADC input

f=600;
T=1/f;
step=1/9000
t=0:step:2000*T;
tri=2.*3.23*asin(sin(2*pi*f.*t))./pi;
cmp=0;
for i=1:length(t)
    if tri(i)>=cmp
        dr(i) = 3.23;
    else
        dr(i) = 0;
    end
end


%% input signal - sinewave 50Hz - from the signal generator -> ADC input
% f=50;
% T=1/f;
% step=1/9000
% t=0:step:200*T;
% dr=(1.615.*sin(2.*pi.*f.*t)+1.615);


%% average amplitude detection
in1=abs(dr).*abs(dr);

figure;
plot(t,dr,'r','linewidth',2);hold on;
plot(t,in1,'b','linewidth',2);hold on;

in1_filt=zeros(1,length(t));
in2_filt=zeros(1,length(t));
 
b0=step*step
b1=2*step*step
b2=step*step

a0=step*step+0.4*step+0.04
a1=2*step*step-0.08
a2=step*step-0.4*step+0.04
a00=1/a0
 

for i=2:length(t)
 
    if i>=3
        in1_filt(i)=a00*(b0*in1(i)+b1*in1(i-1)+b2*in1(i-2)-a1*in1_filt(i-1)-a2*in1_filt(i-2));
    else
        in1_filt(i)=0;      
    end

end


figure;
plot(t,in1_filt,'*b-','linewidth',2);hold on;