typedef struct structIIR{    
    double*    b;
    int        lb;
    double* a;
    int        la;
}IIR, *PIIR;


int conv(double *a, int la, double *b, int lb, double *c){
    if((!a)||(!b)||(!c)||(!la)||(!lb))
        return 0;
    memset(c, 0, (la+lb-1)*sizeof(double));
    for(int i = 0; i < la; i++){
        for(int j = 0; j < lb; j++)
            c[i+j] += a[i]*b[j];
    }
    return 1;
}

int iir_filter(PIIR piir, double *s, int n){
    double* c = (double*)malloc(piir->la*sizeof(double));
    double* d = (double*)malloc(piir->lb*sizeof(double));
    double* v = (double*)malloc((piir->lb+n-1)*sizeof(double));

    for(int i =0; i<piir->la; i++)
        c[i] = piir->a[i]/piir->a[0];
    for(int i =0; i<piir->lb; i++)
        d[i] = piir->b[i]/piir->a[0];
    conv(d,piir->lb,s,n,v);
    s[0] = v[0];
    
    for(int i = 1; i<n; i++){
        int k = 1;
        s[i] = v[i];
        while(((i-k)>=0)&&(k<piir->la)){
            s[i]-=c[k]*s[i-k];
            k++;
        }
    }
    free(c);
    free(d);
    free(v);
    return 1;
}

/*
* Filter Coefficients (C Source) generated by the Filter Design and Analysis Tool
*
* Generated by MATLAB(R) 7.8 and the Signal Processing Toolbox 6.11.
*
* Generated on: 26-Apr-2010 15:26:51
*
*/

/*
* Discrete-Time IIR Filter (real)
* -------------------------------
* Filter Structure    : Direct-Form II, Second-Order Sections
* Number of Sections  : 3
* Stable              : Yes
* Linear Phase        : No
*/

/* General type conversion for MATLAB generated C-code  */
#include "tmwtypes.h"
/*
* Expected path to tmwtypes.h
* C:\Program Files\MATLAB\R2009a\extern\include\tmwtypes.h
*/
#define MWSPT_NSEC 7
const int NL[MWSPT_NSEC] = { 1,3,1,3,1,2,1 };
const real64_T NUM[MWSPT_NSEC][3] = {
  {
   0.001582012905215,                 0,                 0
  },
  {
                   1,   -1.997425031815,                 1
  },
  {
    0.06078064837324,                 0,                 0
  },
  {
                   1,   -1.998996713475,                 1
  },
  {
   0.007688806860118,                 0,                 0
  },
  {
                   1,                 1,                 0
  },
  {
                   1,                 0,                 0
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1,3,1,2,1 };
const real64_T DEN[MWSPT_NSEC][3] = {
  {
                   1,                 0,                 0
  },
  {
                   1,   -1.996305874613,   0.9963329403565
  },
  {
                   1,                 0,                 0
  },
  {
                   1,   -1.998581806418,   0.9986427868237
  },
  {
                   1,                 0,                 0
  },
  {
                   1,  -0.9976855336623,                 0
  },
  {
                   1,                 0,                 0
  }
};

#include "fdacoefs.h"

double WHistory[MWSPT_NSEC][2];

double IIR_filt (double X){
    for (int i=0; i<MWSPT_NSEC; i++){
        double W;
        W = X;
        for (int j=1; j<DL[i]; j++){
            W-= DEN[i][j]*WHistory[i][j-1];
        };
        W*=1/DEN[i][0];

        X = W * NUM[i][0];
        for (int j=1; j<NL[i]; j++){
            X+=NUM[i][j]*WHistory[i][j-1];
        };
        WHistory[i][1] = WHistory[i][0];
        WHistory[i][0] = W;
    };
    return X;
};

#include "fdacoefs.h"

typedef double    real;

#define Fs        3600    // Частота сэмплирования

#define IIR_hist_size        5    // Кол-во ячеек памяти в реализации фильтра

real History[IIR_hist_size];

real IIR_filter (real X, real* Hist){
    for (int i=0; i<MWSPT_NSEC; i++){
        int jmax = ((NL[i]>DL[i])?(NL[i]):(DL[i]))-1;
        if (jmax>0){
            // Блок фильтра из каскада
            real W = 0;
            for (int j=0; j<jmax; j++){
                W += DEN[i][j+1] * Hist[j];
            };
            W = X - W;
            W *= (1/DEN[i][0]);
    
            X = W * NUM[i][0];
            for (int j=0; j<jmax; j++){
                X += NUM[i][j+1] * Hist[j];
            };
            
            for (int j=1; j<jmax; j++){
                Hist[j] = Hist[j-1];
            };
            Hist[0] = W;
            Hist += jmax;
        } else if (jmax==0){
            // Блок "Gain" между каскадами фильтра
            X *= NUM[i][0] * (1/DEN[i][0]);
        };
    };
    return X;
};

int main(){
    // Обнуляем ячейки памяти
    for (int i=0; i<IIR_hist_size; i++){
        History[i]=0;
    };
    
    real X,Y,F;

    F = 10.0; // 10 Гц

    F=2*CONST_PI*F/(real)Fs;

    for (int i=0; i<1000000; i++){
        X=cos(F*(real)i);
        Y=IIR_filter (X,History);
    };

    return 0;
}