psd_floatpt.c

用dsp解压mp3程序的算法

/*****************************************************************
*  psd_floatpt.c - C program for FFT-based PSD estimation     
******************************************************************
*  System configuration:
*
*          _____     ____     __________     __________
*         |     |   |    |   |          |   |          |
*  x(n)-->| Buf |-->| BR |-->| (1/N)*FFT|-->| |X(k)|^2 |--> P(k)
*         |_____|   |____|   |__________|   |__________|                         
*           
******************************************************************
*  System simulation configuration:
*
*     x(n) is the input data from data file "in.dat"
*     out(n) is the power spectrum data to data file "out.dat"
*
*****************************************************************/

#include <stdio.h>        /* header files                       */
#include <stdlib.h>
#include <math.h>   
#include "def_complex.h"  /* complex.h header file              */

/*    external functions used by this program                   */
extern void ditr2fft(complex *, unsigned int, complex *, unsigned int);
extern void bit_reversal(complex *, unsigned int);

/*    define variables and constants                            */
#define N 512             /* number of FFT points               */
#define EXP 9             /* EXP=log2(N)                        */
#define pi 3.1415926535897  
                       
complex X[N];             /* declare input array                */
complex W[EXP];           /* twiddle e^(-j2pi/N) table          */    
complex temp;
float spectrum[N];
float xn;

void main()
{
  unsigned int i,L,LE,LE1; 
  /***************************************************************
  *    Declare file pointers
  ***************************************************************/

  FILE *xn_in;                   /* file pointer of x(n)        */
  FILE *yn_out;                  /* file pointer of y(n)        */
  xn_in = fopen("in.dat","r");   /* open file for input x(n)    */
  yn_out = fopen("out.dat","w"); /* open file for output y(n)   */

  /* Step 1: Create a twiddle factor table                      */

  for (L=1; L<=EXP; L++)         /* create twiddle factor table */
  {
    LE=1<<L;                     /* LE=2^L=points of sub DFT    */
    LE1=LE>>1;     	    /* number of butterflies in sub-DFT */
    W[L-1].re = cos(pi/LE1);
    W[L-1].im = -sin(pi/LE1);
  }    
    
  /* Step 2: Enter input data to reference buffer               */

  for(i=0; i<N; i++)
  {			/* input signal */
    fscanf(xn_in,"%f",&xn);
    {	
      X[i].re = xn;
      X[i].im = 0.0;
    }
  }

  /* Step 3: Perform bit reversal, follow by FFT                */

  bit_reversal(X,EXP);    /* arrange X[] in bit-reverse order   */
  ditr2fft(X,EXP,W,1);    /* FFT with scaling 0.5 in each stage */

  /* Step 4: Perform magnitude-square                           */

  for (i=0; i<N; i++)     /* verify FFT result                  */
  {
    temp.re = X[i].re*X[i].re;
    temp.im = X[i].im*X[i].im;        
    spectrum[i] = (temp.re+temp.im); /* perform scaling  in fft */
    fprintf(yn_out,"%f\n",spectrum[i]);
  }
  fcloseall();
}