DFT

FFT

#include "stdio.h"

#include "math.h"

#include "windows.h"

//FFT参数

#define PI   3.1415926

#define FFT_N    128

float   fs=400;     //采样率 hz

float     Ts=1/fs;    //采样周期 s

//定义复数结构体

struct complex

{

    float Re ; //实部

    float Im ; //虚部

};

//DFT数据结构体

struct FFT_Data

{

    float In;     //输入   

    float OutRe;  //输出实部

    float OutIm;  //输出虚部

    float OutAmp; //幅值

    float OutPha; //相位     

};

FFT_Data     gsta_FFT_Data[FFT_N] ;        //FTT数据数组

complex      gsta_FFT_DataTemp[FFT_N] ; //FFT计算数组 计算前保存f[n] 计算后保存F[k]

void FFT_WriteData(void)

{

    for(int n=0;n<FFT_N;n++)

    {

        gsta_FFT_Data[n].In = 2 + 3*cos(50*2*PI*Ts*n - PI/6) + 0.8*cos(100*2*PI*Ts*n + PI/3);

        gsta_FFT_Data[n].OutRe = 0;

        gsta_FFT_Data[n].OutIm = 0;

        gsta_FFT_Data[n].OutAmp = 0;

        gsta_FFT_Data[n].OutPha = 0;

        gsta_FFT_DataTemp[n].Re = gsta_FFT_Data[n].In;

        gsta_FFT_DataTemp[n].Im = 0;

    }   

}

void FFT_ReadData(void)

{

    //读取输入数据

    printf("FFT输入f[n]:\n");

    for(int n=0;n<FFT_N;n++)

    {

        printf("%f\t",gsta_FFT_Data[n].In);

        if((n % 8) == 7) //一行显示8个点

        {

            printf("\n");

        }       

    }

    //读取输出数据

    printf("\nFFT输出F[k]:\n");

    for(int k=0;k<FFT_N;k++)

    {

        printf("F[%d]\t %fHz\t Re:%f\t Im:%f\t Amp:%f\t Pha:%f \n",k,(float)k*fs/FFT_N,gsta_FFT_Data[k].OutRe,gsta_FFT_Data[k].OutIm,gsta_FFT_Data[k].OutAmp,gsta_FFT_Data[k].OutPha);

    }   

}

// 快速傅里叶变换

void FFT_Calculate(complex *Data, int N)

{

    if (N <= 1) return;

    // 分割奇偶序列

    complex even[N/2];

    complex odd[N/2];

    for (int i = 0; i < N/2; i++)

    {

        even[i] = Data[2*i];

        odd[i]  = Data[2*i + 1];

    }

    // 递归调用--对奇偶序列分别FFT

    FFT_Calculate(even,N/2);

    FFT_Calculate(odd,N/2);

    float w_N = 2*PI/N;  //N点FFT相邻采样点间角度

    //合并结果

    for (int k = 0; k < N/2; k++)

    {

        complex efujkwN_Mult_oddk;

        efujkwN_Mult_oddk.Re =  odd[k].Re * cos(k*w_N) + odd[k].Im * sin(k*w_N) ;

        efujkwN_Mult_oddk.Im = - odd[k].Re * sin(k*w_N) + odd[k].Im * cos(k*w_N);

        //Data[k + N/2] = even[k] + efujkwN_Mult_oddk*odd[k];

        Data[k].Re = even[k].Re + efujkwN_Mult_oddk.Re;

        Data[k].Im = even[k].Im  + efujkwN_Mult_oddk.Im;

        //Data[k + N/2] = even[k] - efujkwN_Mult_oddk*odd[k];

        Data[k + (N/2)].Re =even[k].Re - efujkwN_Mult_oddk.Re;

        Data[k + (N/2)].Im =even[k].Im - efujkwN_Mult_oddk.Im;

    }

}

void FFT_AnalyzeData()

{

    //直流分量

    gsta_FFT_Data[0].OutRe = gsta_FFT_DataTemp[0].Re;

    gsta_FFT_Data[0].OutIm = gsta_FFT_DataTemp[0].Im;

    gsta_FFT_Data[0].OutAmp = sqrt(gsta_FFT_Data[0].OutRe*gsta_FFT_Data[0].OutRe + gsta_FFT_Data[0].OutIm*gsta_FFT_Data[0].OutIm) / FFT_N ;       

    gsta_FFT_Data[0].OutPha =  atan2f(gsta_FFT_Data[0].OutIm,gsta_FFT_Data[0].OutRe);

    //fs k次谐波分量

    for(int k=1;k<FFT_N;k++)

    {       

        gsta_FFT_Data[k].OutRe = gsta_FFT_DataTemp[k].Re;

        gsta_FFT_Data[k].OutIm = gsta_FFT_DataTemp[k].Im;

        gsta_FFT_Data[k].OutAmp = sqrt(gsta_FFT_Data[k].OutRe*gsta_FFT_Data[k].OutRe + gsta_FFT_Data[k].OutIm*gsta_FFT_Data[k].OutIm) / FFT_N * 2;       

        gsta_FFT_Data[k].OutPha =  atan2f(gsta_FFT_Data[k].OutIm,gsta_FFT_Data[k].OutRe);

    }       

}

//运行时间测试程序

void FFT_TestRunTime(void)

{

    double         TimeRun;     //运行期间定时器计数值

    _LARGE_INTEGER TimeStart;    //开始时定时器计数值

    _LARGE_INTEGER TimeEnd;        //结束时定时器计数值

    _LARGE_INTEGER TimeFreq;    //定时器频率

    QueryPerformanceFrequency(&TimeFreq);

    FFT_WriteData();

    QueryPerformanceCounter(&TimeStart);    //计时开始

    FFT_Calculate(gsta_FFT_DataTemp, FFT_N);

    FFT_AnalyzeData();

    QueryPerformanceCounter(&TimeEnd);    //计时结束

    TimeRun = 1000000*(TimeEnd.QuadPart-TimeStart.QuadPart)/TimeFreq.QuadPart;

    printf("\nN=%d FFT_RunTime: %fus\n",FFT_N,TimeRun);       

}

int main()

{ 

    FFT_WriteData();       

    FFT_Calculate(gsta_FFT_DataTemp, FFT_N);

    FFT_AnalyzeData();

    FFT_ReadData();

    FFT_TestRunTime();

    return 0;

}