傅立叶变换—DFT

背景：最近看到实验室其他同学在用傅立叶变换解决问题，我也想通过并行来解决这个问题，所以看了一下傅立叶变换的东西，感觉涵盖的东西还能多，我只是初步做了一下了解（一定很片面，但是我主要是为了应用它，主要了解它的实现原理），据我理解：傅立叶分析使用一系列sin函数和cos函数表示一个连续函数。傅立叶变化就是一个将时间域的函数序列f(k)映射到频率域上的函数序列F(j)。序列f(k)表示一组信号采样的时间函数，序列F(j)表示傅立叶系数的分布，这个分布是关于频率的函数。

傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。

离散傅立叶变换（DFT）是一个矩阵-向量的乘积F_Nx,其中f_ij=w_n^ij,(i,j从0到n)，w_n是单位元的第n个本原根。

复数知识回顾：

（1）单位1的n次复根是复数w,即w_n=1;

（2）单位1的n次复根共有n个，即e^2∏*i*k/n,其中k=1,2,...,n。

（3）用w_n来表示复数e^2∏*i/n，它是单位1的主n次方根。

（4）消去原理：对于任意整数n>=0,k>=0和d>0,有w_dn^dk=w_n^k。

（5）对于任意偶数n>0,有w_n^n/2 = w₂= -1

（6）等分引理：如果n为正偶数，单位1的n个n次方根的平方与它的（n/2）个n/2次方根是相同的。

DFT代码：

#include "stdio.h"
#include "math.h"

typedef struct DFT
{
  float real;  //实部
  float img;  //inscribrer
}DFT;

#define PI 3.1415926
#define LENGTH  4
int main()
{
  DFT W[LENGTH][LENGTH],w[LENGTH+1];  
  float var = 2  /(float)LENGTH;
  
  w[0].real = 1;
  w[0].img = 0;
  w[1].real = cos(var * PI);
  w[1].img = sin(var * PI);
  
  int i,j;
  
    //初始化w0到wn-1的值
  for(i=2;i<=LENGTH;i++)
  {
    w[i].real = cos(var*i*PI);
    w[i].img = sin(var*i*PI);
  }
  
  for(i=0;i<LENGTH;i++)
  {
    for(j=0;j<LENGTH;j++)
    {
      int m = i * j;
      if((i*j)<=LENGTH)
      {
    W[i][j].real = w[m].real;
    W[i][j].img = w[m].img;
      }
      else 
      {
    int n = m % LENGTH;  //w具有周期性
    W[i][j].real = w[n].real;
    W[i][j].img = w[n].img;
      }
    }
  }
/*  
  for(i=0;i<WIDTH;i++)
  {
    for(j=0;j<WIDTH;j++)
    {
      printf("%f	",W[i][j].real);
    }
    printf("
");
  }
*/
// float a[WIDTH] = {1,2,3,4};
 float a[LENGTH] = {1,2,3,4};
 DFT result[LENGTH];
  
 DFT sum;
    //w矩阵与初始一维矩阵a相乘
 for(i=0;i<LENGTH;i++)
 {
   sum.real = 0;
   sum.img = 0;
   for(j=0;j<LENGTH;j++)
   {
     sum.real += W[i][j].real * a[j];
     sum.img += W[i][j].img * a[j];
   }
   result[i].real = sum.real;
   result[i].img = sum.img;
 }
 
 for(i=0;i<LENGTH;i++)
  {
    if(result[i].real != 0)
    {
      printf("%3.1f",result[i].real);
    }
    if(result[i].img != 0.0  )
    {
      printf("+%3.1fi",result[i].img);
    } 
    printf("
");
  }
}
   时间复杂度：O（n2）