psd_fixpt_intr.c

用dsp解压mp3程序的算法

/*****************************************************************
*  psd_fixpt_intr.c - C program computes FFT-based PSD using
*                     C5000 intrinsics on CCS
******************************************************************
*  System configuration:
*
*          _____     ____     __________     __________
*         |     |   |    |   |          |   |          |
*  x(n)-->| Buf |-->| BR |-->| (1/N)*FFT|-->| |X(k)|^2 |--> P(k)
*         |_____|   |____|   |__________|   |__________|                         
*           
*  (a) Periodigram
******************************************************************
*  System simulation configuration:
*
*     x(n) is the input data from data file "in_int_ccs.dat"
*
*****************************************************************/
/* Header files  */

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

int in_buffer[512];
#define TRUE 1 

/* Function */
static void dataIO(void);

/* External functions  */
extern void ditr2fft_fixpt_intr(complex *, unsigned int, complex *, unsigned int);
extern void ibit_reversal_intr(complex *, unsigned int);

/*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 */    
lcomplex ltemp;
int spectrum[N];
int spectrum_db[N];

void main()
{
    unsigned int i,L,LE,LE1; 
    int *input = &in_buffer[0]; 
	
    /* 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 = (int)((0x7fff*cos(pi/LE1))+0.5);
        W[L-1].im = -(int)((0x7fff*sin(pi/LE1))+0.5);
    }    
    
    /* Step 2: Enter input data to reference buffer */
    
    dataIO();
    for(i=0; i<N; i++)
    { 	 
	X[i].re = *input++;
	X[i].im = 0;	
    } 

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

    /* Start FFT */
    ibit_reversal_intr(X,EXP);       /* Arrange X[] in bit-reverse order */
    ditr2fft_fixpt_intr(X,EXP,W,1);  /* Perform FFT with scaling of 0.5 in each stage */


    /* Step 4: Perform Magnitude-Square and display results in MATLAB */

    for (i=0; i<N; i++) /* Verify FFT result */
    {
 	ltemp.re = _lsmpy(X[i].re,X[i].re);
  	ltemp.im = _lsmpy(X[i].im,X[i].im);  
 	spectrum[i] = (int)(_lsadd(ltemp.re,ltemp.im)>>13);/*Scaling is perform in fft*/       
        spectrum_db[i] = 10*log10(spectrum[i]+1);
    }
} 

/* Function for dataIO */

static void dataIO()
{
    return;
}