fsf.c

用于DSK5416的程序

/***************************************************************
*    FSF.C - FREQUENCY SAMPLING FILTER                         *
***************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define K 1001

void main(int argc, char *argv[])
{
    FILE  *fpxn;                   // file pointer of x(n) data file
    FILE  *fpyn;                   // file pointer of y(n) data file
    float d,d1[K],d2[K],Hk[K],y;
    float *x;
    float c1,c2,c3,c4[K],f,f1;
    float r,rN, pi = 3.14159;
    int k,i,iN,M;
    int n=0;
    int N=200;                     // filter order

/*  Get parameters    */

    if (argc < 2)
    {
        fprintf(stderr,"Usage: fsf infile outfile\n");
        exit(0);
    }
    fpxn = fopen(argv[1],"rb");    // open x(n) file
    fpyn = fopen(argv[2],"wb");    // open y(n) file

/*  Define constants    */

    r = 1.-pow(2.,-26.);           // radius
    rN=pow(r,(float)N);            // r^L
    c1 = 2./N;                     // 2/L in (5.3.42)
    c2=c1*rN;                      // 2*r^L/L in (5.3.42)
    c3=r*r;                        // r^2 in (5.3.43)
    for (k=0;k<N/2;k++)
    {
        c4[k] = r * cos(2.*pi*k/(float)N); // coefs in (5.3.43)
        Hk[k] = 0.;                        // gain Hk
    }
    Hk[20] = Hk[25] = 0.707;               // nonzero gains
    Kk[21] = Hk[22] = Hk[23] = Hk[24] = 1.;
    M = N + 2;                     // size of circular buffer
    
/*  Initialize parameters    */

    x = (float *)calloc(M,sizeof(float)); // define x(n) buffer
    i = N+1;                       // circular pointer for x(n)
    iN = 1;                        // circular pointer for x(n-L)
    f1 = 0.;                       // u(n-1)
    for (k=0;k<N/2;k++)
    {
	d1[k] = 0.;                // fk(n-1)
	d2[k] = 0.;                // fk(n-2)
    }

/*****************************************************************
*      MAIN FILTER LOOP                                          *
*****************************************************************/
    while( fread(&x[i],sizeof(float),1,fpxn) == 1)
    {
	y = 0.;

	f = c1*x[i] - c2*x[iN];  // u(n) in (5.3.42)
	for(k=1; k<N/2; k++)
	{
	    d = f + c4[k]*(2.*d1[k] - f1) - c3*d2[k]; // fk(n) in (5.3.43)
	    y += d*Hk[k];        // y(n) in (5.3.44)

/*  Update buffers  */
	    d2[k] = d1[k];       // update fk(n-2) buffer
	    d1[k] = d;           // update fk(n-1) buffer
	}
	f1 = f;                  // update u(n-1)
	fwrite(&y,sizeof(float),1,fpyn); // write output to data file

/*  Update circular pointers for x[n] buffer    */
	i = (i+1)%M;             // update circular pointer for x(n)
	iN = (iN+1)%M;           // update circular pointer for x(n-N)
    }
}