Internal Technical Memo

98-11-11

 

An Implementation of IIR and Zero Phase Filtering

 

Date: 1998.11.11

Presented by

Yangquan Chen and Huifang Dou

 

 

  1. Purposes: Develop an ANSI C version of the two useful scripts in MATLAB, i.e., FILTER() and FILTFILT().
  2. C Codes:

    3. /*

    FILTER.C

    An ANSI C implementation of MATLAB FILTER.M (built-in)

    Written by Chen Yangquan <elecyq@nus.edu.sg>

    1998-11-11

    */

     

    #include<stdio.h>

    #define ORDER 3

    #define NP 1001

     

    /*

    void filter(int,float *,float *,int,float *,float *);

    */

     

    filter(int ord, float *a, float *b, int np, float *x, float *y)

    {

    int i,j;

    y[0]=b[0]*x[0];

    for (i=1;i<ord+1;i++)

    {

    y[i]=0.0;

    for (j=0;j<i+1;j++)

    y[i]=y[i]+b[j]*x[i-j];

    for (j=0;j<i;j++)

    y[i]=y[i]-a[j+1]*y[i-j-1];

    }

    /* end of initial part */

    for (i=ord+1;i<np+1;i++)

    {

    y[i]=0.0;

    for (j=0;j<ord+1;j++)

    y[i]=y[i]+b[j]*x[i-j];

    for (j=0;j<ord;j++)

    y[i]=y[i]-a[j+1]*y[i-j-1];

    }

    } /* end of filter */

     

     

     

     

     

     

     

    main()

    {

    FILE *fp;

    float x[NP],y[NP],a[ORDER+1],b[ORDER+1];

    int i,j;

     

    /* printf("hello world \n"); */

     

    if((fp=fopen("acc1.dat","r"))!=NULL)

    {

    for (i=0;i<NP;i++)

    {

    fscanf(fp,"%f",&x[i]);

    /* printf("%f\n",x[i]); */

    }

    }

    else

    {

    printf("\n file not found! \n");

    exit(-1);

    }

    close(fp);

     

    /* test coef from

    [b,a]=butter(3,30/500); in MATLAB

    */

    b[0]=0.0007;

    b[1]=0.0021;

    b[2]=0.0021;

    b[3]=0.0007;

    a[0]=1.0000;

    a[1]=-2.6236;

    a[2]=2.3147;

    a[3]=-0.6855;

     

    filter(ORDER,a,b,NP,x,y);

    /* NOW y=filter(b,a,x);*/

     

    /* reverse the series for FILTFILT */

    for (i=0;i<NP;i++)

    { x[i]=y[NP-i-1];}

    /* do FILTER again */

    filter(ORDER,a,b,NP,x,y);

    /* reverse the series back */

    for (i=0;i<NP;i++)

    { x[i]=y[NP-i-1];}

    for (i=0;i<NP;i++)

    { y[i]=x[i];}

    /* NOW y=filtfilt(b,a,x); boundary handling not included*/

     

    if((fp=fopen("acc10.dat","w+"))!=NULL)

    {

    for (i=0;i<NP;i++)

    {

    fprintf(fp,"%f\n",y[i]);

    }

    }

    else

    {

    printf("\n file cannot be created! \n");

    exit(-1);

    }

    close(fp);

    }

    /* end of filter.c */

     

     

  3. Application Demo

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Fig. 1 Original acceleration signal. Fig. 2 Filtered acceleration signal.

 

 

Appendix: MATLAB .m file listing

 

A.1 FILTER.M

 

FILTER Digital filter.

Y = FILTER(B, A, X) filters the data in vector X with the

filter described by vectors A and B to create the filtered

data Y. The filter is a "Direct Form II Transposed"

implementation of the standard difference equation:

y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)

- a(2)*y(n-1) - ... - a(na+1)*y(n-na)

 

[Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final

conditions, Zi and Zf, of the delays.

See also FILTFILT in the Signal Processing Toolbox.

 

A.2 FILTFILT.M

function y = filtfilt(b,a,x)

%FILTFILT Zero-phase forward and reverse digital filtering.

% Y = FILTFILT(B, A, X) filters the data in vector X with the

% filter described by vectors A and B to create the filtered

% data Y. The filter is described by the difference equation:

%

% y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)

% - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

%

% After filtering in the forward direction, the filtered

% sequence is then reversed and run back through the filter.

% The resulting sequence has precisely zero-phase distortion

% and double the filter order. Care is taken to minimize

% startup and ending transients by matching initial conditions.

%

% The length of the input x must be more than three times

% the filter order, defined as max(length(b)-1,length(a)-1).

%

% See also FILTER.

% Author(s): L. Shure, 5-17-88

% revised by T. Krauss, 1-21-94

% initial conditions: Fredrik Gustafsson

% Copyright (c) 1984-94 by The MathWorks, Inc.

% $Revision: 1.9 $ $Date: 1994/01/25 17:59:07 $

error(nargchk(3,3,nargin))

if (isempty(b)|isempty(a)|isempty(x))

y = [];

return

end

[m,n] = size(x);

if (n>1)&(m>1)

error('Only works for vector input.')

end

if m==1

x = x(:); % convert row to column

end

len = size(x,1); % length of input

b = b(:).';

a = a(:).';

nb = length(b);

na = length(a);

nfilt = max(nb,na);

nfact = 3*(nfilt-1); % length of edge transients

if (len<=nfact), % input data too short!

error('Data must have length more than 3 times filter order.');

end

% set up filter's initial conditions to remove dc offset problems at the

% beginning and end of the sequence

if nb < nfilt, b(nfilt)=0; end % zero-pad if necessary

if na < nfilt, a(nfilt)=0; end

zi = ( eye(nfilt-1) - [-a(2:nfilt).' [eye(nfilt-2); zeros(1,nfilt-2)]] ) \ ...

( b(2:nfilt).' - a(2:nfilt).'*b(1) );

% Extrapolate beginning and end of data sequence using a "reflection

% method". Slopes of original and extrapolated sequences match at

% the end points.

% This reduces end effects.

y = [2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)];

% filter, reverse data, filter again, and reverse data again

y = filter(b,a,y,[zi*y(1)]);

y = y(length(y):-1:1);

y = filter(b,a,y,[zi*y(1)]);

y = y(length(y):-1:1);

% remove extrapolated pieces of y

y([1:nfact len+nfact+(1:nfact)]) = [];

if m == 1

y = y.'; % convert back to row if necessary

end

<End of this MEMO>