C语言直接DFT和FFTW两种方式实现离散傅里叶变换

  /*
   * dft.c
   *
   *  Created on: 2018年5月16日
   *      Author: ChiWang
   */
  
  #include <stdio.h>
  #include <math.h>
  #include "fftw3.h"
  
  
  #define N   16             //采样个数
  #define PI  3.1415926535
  typedef struct              //complex数据类型，用于实现傅里叶运算
  {
    float r;
    float i;
  }complex;
  complex dft_out[N];      //complex out data
  float  dft_in[N];  	//sample in data
  float  amp[N]; 	//real_amp=amp*2/N
  
  void DFT_Cal(void)
  {
      int i = 0;
      int n = 0;
      complex part[N];
  
  
      for(i=0; i<N; i++)
      {
	  dft_out[i].r = 0;
	  dft_out[i].i = 0;
	  for(n=0; n<N; n++)
	  {
	      //欧拉公式 cos(x)-jsin(x)
	      part[n].r =  cos(2*PI*i/N*n)*dft_in[n];
	      part[n].i = -sin(2*PI*i/N*n)*dft_in[n];
  
	      dft_out[i].r += part[n].r;
	      dft_out[i].i += part[n].i;
	  }
	  amp[i] = sqrt(dft_out[i].r * dft_out[i].r + dft_out[i].i * dft_out[i].i);
      }
  }
  int main(int argc, char *argv[])
  {
      int i = 0;
      //
      for(i=0; i<N; i++)//假设周期是工频50Hz，那么T=0.02,那么根据采样定律，一个周波采样2个点就行，这里取一个周波采样N=16个点，那么采样率Fs=N*50=50N，
      {
	  //220sin(wt) -> f=50Hz,w=2*Pi/T -> w=2*pi*f -> 220sin(2*pi*f*t)
	  dft_in[i] = 220*sin(2*PI*50*i*0.02/N);//一个周期采样N个点，那么一个点的时间就是T/N=0.02/N
	  printf("%lf ",dft_in[i]);
      }
  
      /* fftw related init */
      float *fft_in_buffer=NULL;
      fftwf_complex *fft_out_buffer = NULL;
      fftwf_plan fft_plan;
      size_t FFT_len = N;
      fft_in_buffer  = (float *)fftwf_malloc(sizeof(float) * FFT_len);
      fft_out_buffer = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * FFT_len);
  
  #define FFTW 1
  #if FFTW
      fft_plan = fftwf_plan_dft_r2c_1d(FFT_len, dft_in, fft_out_buffer, FFTW_ESTIMATE);
      fftwf_execute(fft_plan);
  #else
      DFT_Cal();//进行DFT计算
  #endif
  
      printf("\n\n");
      for(i=0; i<N; i++)
	{
	  //x轴 = 0: 1/T : (n-1)/T, 或者x=i*Fs/N, amp=out/2*N
	  //根据香农定律，只有Fs/2有效，后半部分的频谱是对称的共轭部分
  #if FFTW
	  double AMP = sqrt(fft_out_buffer[i][0]*fft_out_buffer[i][0] + fft_out_buffer[i][1]*fft_out_buffer[i][1]);
	  printf("(freq=%f,amp=%f)\n",i*1./0.02,AMP*2/N);
  #else
	  printf("(freq=%f,amp=%f)\n",i*1./0.02,amp[i]*2/N);
  #endif
	}
      return 0;
  }

1，使能FFTW，运行输出（#define FFTW 1，并下载libfftw3f-3.dll，然后编译时连接该库）：

0.000000 84.190353 155.563492 203.253494 220.000000 203.253494 155.563492 84.190353 0.000000 -84.190353 -155.563492 -203.253494 -220.000000 -203.253494 -155.563492 -84.190353 
(freq=0.000000,amp=0.000000)
(freq=50.000000,amp=220.000000)
(freq=100.000000,amp=0.000000)
(freq=150.000000,amp=0.000002)
(freq=200.000000,amp=0.000000)
(freq=250.000000,amp=0.000002)
(freq=300.000000,amp=0.000000)
(freq=350.000000,amp=0.000006)
(freq=400.000000,amp=0.000000)
(freq=450.000000,amp=0.000000)
(freq=500.000000,amp=0.000000)
(freq=550.000000,amp=0.000000)
(freq=600.000000,amp=0.000000)
(freq=650.000000,amp=0.000000)

(freq=700.000000,amp=0.000000)

2，直接DFT计算运行输出（#define FFTW 0）：

0.000000 84.190353 155.563492 203.253494 220.000000 203.253494 155.563492 84.190353 0.000000 -84.190353 -155.563492 -203.253494 -220.000000 -203.253494 -155.563492 -84.190353 


(freq=0.000000,amp=0.000002)
(freq=50.000000,amp=220.000000)
(freq=100.000000,amp=0.000002)
(freq=150.000000,amp=0.000001)
(freq=200.000000,amp=0.000000)
(freq=250.000000,amp=0.000001)
(freq=300.000000,amp=0.000000)
(freq=350.000000,amp=0.000004)
(freq=400.000000,amp=0.000000)
(freq=450.000000,amp=0.000004)
(freq=500.000000,amp=0.000000)
(freq=550.000000,amp=0.000001)
(freq=600.000000,amp=0.000000)
(freq=650.000000,amp=0.000001)
(freq=700.000000,amp=0.000002)
(freq=750.000000,amp=220.000000)