傅立叶变换—DFT
背景:最近看到实验室其他同学在用傅立叶变换解决问题,我也想通过并行来解决这个问题,所以看了一下傅立叶变换的东西,感觉涵盖的东西还能多,我只是初步做了一下了解(一定很片面,但是我主要是为了应用它,主要了解它的实现原理),据我理解:傅立叶分析使用一系列sin函数和cos函数表示一个连续函数。傅立叶变化就是一个将时间域的函数序列f(k)映射到频率域上的函数序列F(j)。序列f(k)表示一组信号采样的时间函数,序列F(j)表示傅立叶系数的分布,这个分布是关于频率的函数。
傅立叶原理表明:任何连续测量的时序或信号,都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号,以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。
离散傅立叶变换(DFT)是一个矩阵-向量的乘积FNx,其中fij=wnij,(i,j从0到n),wn是单位元的第n个本原根。
复数知识回顾:
(1)单位1的n次复根是复数w,即wn=1;
(2)单位1的n次复根共有n个,即e2∏*i*k/n,其中k=1,2,...,n。
(3)用wn来表示复数e2∏*i/n,它是单位1的主n次方根。
(4)消去原理:对于任意整数n>=0,k>=0和d>0,有wdndk=wnk。
(5)对于任意偶数n>0,有wnn/2 = w2= -1
(6)等分引理:如果n为正偶数,单位1的n个n次方根的平方与它的(n/2)个n/2次方根是相同的。
DFT代码:
1 #include "stdio.h" 2 #include "math.h" 3 4 typedef struct DFT 5 { 6 float real; //实部 7 float img; //inscribrer 8 }DFT; 9 10 #define PI 3.1415926 11 #define LENGTH 4 12 int main() 13 { 14 DFT W[LENGTH][LENGTH],w[LENGTH+1]; 15 float var = 2 /(float)LENGTH; 16 17 w[0].real = 1; 18 w[0].img = 0; 19 w[1].real = cos(var * PI); 20 w[1].img = sin(var * PI); 21 22 int i,j; 23
//初始化w0到wn-1的值 24 for(i=2;i<=LENGTH;i++) 25 { 26 w[i].real = cos(var*i*PI); 27 w[i].img = sin(var*i*PI); 28 } 29 30 for(i=0;i<LENGTH;i++) 31 { 32 for(j=0;j<LENGTH;j++) 33 { 34 int m = i * j; 35 if((i*j)<=LENGTH) 36 { 37 W[i][j].real = w[m].real; 38 W[i][j].img = w[m].img; 39 } 40 else 41 { 42 int n = m % LENGTH; //w具有周期性 43 W[i][j].real = w[n].real; 44 W[i][j].img = w[n].img; 45 } 46 } 47 } 48 /* 49 for(i=0;i<WIDTH;i++) 50 { 51 for(j=0;j<WIDTH;j++) 52 { 53 printf("%f\t",W[i][j].real); 54 } 55 printf("\n"); 56 } 57 */ 58 // float a[WIDTH] = {1,2,3,4}; 59 float a[LENGTH] = {1,2,3,4}; 60 DFT result[LENGTH]; 61 62 DFT sum;
//w矩阵与初始一维矩阵a相乘 63 for(i=0;i<LENGTH;i++) 64 { 65 sum.real = 0; 66 sum.img = 0; 67 for(j=0;j<LENGTH;j++) 68 { 69 sum.real += W[i][j].real * a[j]; 70 sum.img += W[i][j].img * a[j]; 71 } 72 result[i].real = sum.real; 73 result[i].img = sum.img; 74 } 75 76 for(i=0;i<LENGTH;i++) 77 { 78 if(result[i].real != 0) 79 { 80 printf("%3.1f",result[i].real); 81 } 82 if(result[i].img != 0.0 ) 83 { 84 printf("+%3.1fi",result[i].img); 85 } 86 printf("\n"); 87 } 88 }
时间复杂度:O(n2)