experiment for fft convolution using overlap-add method.c

采用FFT和IFFT的快速卷积的C程序实现源代码。

/* 
    exp7d.c - Experiment for FFT convolution 
               using Overlap-Add method
*/

#include "icomplex.h"

#include "input7.dat"       /* This is the data file used in Chapter 5 */

#define L   48               /* Filter length */
#define N   1024             /* FFT size, N >= L+M-1 */   
#define EXP 10               /* MM = log2(N) */
#define M   (N-L+1)          /* Data sample length */

#pragma DATA_SECTION(U, "fft_twd");
#pragma DATA_SECTION(X, "fft_temp");  
#pragma DATA_SECTION(H, "fft_coef");  
#pragma DATA_SECTION(LP_h, "fft_coef");  
#pragma DATA_SECTION(OVRLAP, "fft_in"); 
#pragma DATA_SECTION(re1, "fft_temp");  
#pragma DATA_SECTION(re2, "fft_temp");  
#pragma CODE_SECTION(main, "fft_code"); 

complex X[N];               /* Signal buffer */
complex H[N];               /* Frequency domain Filter */
complex U[N];               /* Twiddle e^(-j2pi/N) table */    
int OVRLAP[L-1];            /* Frequency domain overlap buffer */  
int re1[N],re2[N]; 

/* Low-pass FIR filter coefficients */
#include "firlp48.dat"      /* This is the same filter used in Chaper 5 */
//#include "firlp8.dat"
//#include "firlp16.dat"
//#include "firlp32.dat"
//#include "firlp64.dat" 
//#include "firlp128.dat"
//#include "firlp256.dat"
//#include "firlp512.dat"

extern void fft(complex *, unsigned int, complex *, int);
extern void bit_rev(complex *, unsigned int);
extern void w_table(complex *, unsigned int);                   
extern void olap_add(complex *, int *, unsigned int, unsigned int, unsigned int);
extern void freqflt(complex *, complex *, unsigned int);

main()
{ 
    unsigned int i,j; 
    
    /* Initialization */ 
    for (i=0; i<L-1; i++)  /* Initialize overlap buffer[L-1] */
        OVRLAP[i] =0;                                    
    
    for (i=0; i<L; i++)    /* Copy filter coefficient to work buffer */   
    {
        X[i].re = LP_h[i];                                             
        X[i].im = 0;
    }

    for (i=i; i<N; i++)    /* Fill zeros to the remaining work buffer */
    {
        X[i].re = 0;                                               
        X[i].im = 0;
    }
    
    w_table(U,EXP);        /* Create Twiddle lookup table for FFT */
    bit_rev(X,EXP);        /* Bit reversal arrange the coefficient */
    fft(X,EXP,U,1);        /* FFT to the filter coefficients */

    for (i=0; i<N; i++)    /* Save frequency domain coefficients */    
    {
	    H[i].re = X[i].re <<EXP;
	    H[i].im = X[i].im <<EXP;
    }
        	
    /* Start FFT Convolution test */   
    j=0;
    for (;;)
    {
        for (i=0; i<M; i++)  
        {
            X[i].re = input[j++]; /* Generate input samples */
            X[i].im = 0;
            re1[i] = X[i].re;     /* Copy to reference buffer */
            if (j==160)
                j=0;
    	  }  
    	
    	  for (i=i; i<N; i++)       /* Fill zeros to data buffer */
        {
            X[i].re = 0;
            X[i].im = 0;
        }    	
    	
        /* Start FFT convolution */
        bit_rev(X,EXP);           /* Arrange sample in bit reversal order */   
        fft(X,EXP,U,1);           /* Perform FFT */ 
        
        freqflt(X,H,N);           /* Perform frequency domain filtering */
        
        bit_rev(X,EXP);           /* Arrange sample in bit reversal order */
        fft(X,EXP,U,0);           /* Perform IFFT */   
        olap_add(X,OVRLAP,L,M,N); /* Overlap and add algorithm */
        /* end of FFT convolution */
                
        for (i=0; i<M; i++)       /* Save filtered data for reference */
        {
            re2[i] = X[i].re;
        }                  
    }    
}