DM642上信号的FFT变换

只看该作者 · 2007-9-27 18:45

#include <stdio.h>
#include <stdlib.h>
#include <csl.h>
#include <csl_timer.h>
#include <math.h>

/* ======================================================================== */
/* Include the DSPLIB header file                                           */
/* ======================================================================== */
#include "dsp_fft32x32.h"

/* ======================================================================== */
/* Macro definition                                                         */
/* ======================================================================== */
#define PI  (3.141592654)
#define NN  (512)

/* ======================================================================== */
/* Input and output arrays                                                  */
/* ======================================================================== */
#pragma DATA_ALIGN(x, 8)
#pragma DATA_ALIGN(xx, 8)
#pragma DATA_ALIGN(y, 8)
#pragma DATA_ALIGN(w, 8)
int x[2*NN];
int w[2*NN];
int xx[2*NN];
int y[2*NN];
unsigned int z[2*NN];
short ch[NN];
/* ======================================================================== */
/*  D2I -- Truncate a 'double' to a 'int',   with clamping.                 */
/* ======================================================================== */
static int d2i(double d)
{
    if (d >=  2147483647.0) return (int)0x7FFFFFFF;
    if (d <= -2147483648.0) return (int)0x80000000;
    return (int)d;
}

int gen_twiddle_fft32x32(int *w, int n, double scale)
{
    int i, j, k, s=0, t;

    for (j = 1, k = 0; j < n >> 2; j = j << 2, s++)
    {
        for (i = t=0; i < n >> 2; i += j, t++)
        {
            w[k +  5] = d2i(scale * cos(6.0 * PI * i / n));
            w[k +  4] = d2i(scale * sin(6.0 * PI * i / n));

            w[k +  3] = d2i(scale * cos(4.0 * PI * i / n));
            w[k +  2] = d2i(scale * sin(4.0 * PI * i / n));

            w[k +  1] = d2i(scale * cos(2.0 * PI * i / n));
            w[k +  0] = d2i(scale * sin(2.0 * PI * i / n));

            k += 6;
        }
    }

    return k;
}

void main(void)
{
    FILE *file;

    short iRead;
    int i;

    int tempI,tempR;

    /* ==================================================================== */
    /* Initialize Chip Support Library                                      */
    /* ==================================================================== */
    CSL_init();

    file = fopen("test1.PCM", "rb");

    iRead = fread(ch, 1, NN, file);
    //iRead = fread(ch, 1, NN, file);

    fclose(file);

    /* ==================================================================== */
    /* Generate Q.31 input data                                             */
    /* ==================================================================== */
    for(i=0; i<NN; i++)
    {
        x[2*i] = xx[2*i] = (int)ch;
        x[2*i+1] = xx[2*i+1] = 0;
    }

    /* ==================================================================== */
    /* Generate twiddle factors                                             */
    /* ==================================================================== */
    gen_twiddle_fft32x32(w, NN, 2147483647.);
    DSP_fft32x32(w, NN, x, y);

    for(i=0; i<NN; i++)
    {    tempI=y[i*2];
        tempR=y[i*2+1];
        z = tempI * tempI + tempR * tempR;

    }


    printf("END\n" );

}

/* ======================================================================== */
/*  End of file: fft32x32_main.c                                            */
/* ======================================================================== */

test1.pcm是8K采样的16BIT语音信号文件。用cool edit pro软件产生。实际是一个440hz正弦信号。
问题描述(用的是CCS3.1环境）：
1、当信号幅值很小的时候，运行程序到末尾后用view/graph看z数组。发现会出现两个对称的峰值。第一个落在z[28]附件，另外一个落在z[480]附件。为什么这样？这是一个正弦波
   。虽然z没有开平方，但是应该是单一峰值啊。
2、当增大波形的幅值，频率不变后，再看z。发现峰值的波形变宽变乱。为什么这个算法和输入的动态范围有关？如何让FFT变换和输入幅值没有关系？
3、把iRead = fread(ch, 1, NN, file);语句重复，也就是ch的值取文件的第二次读写结果进行变换，z的结果完全乱了。按照道理，波形文件足够长，又是一个正弦波形。
   应该没有变化才对。我怀疑什么地方错了。但是我用view/graph看xx的FFT，没有改变。这么说变换前数据是一样的。这是为什么？