来自:http://blog.csdn.net/calcular/article/details/46804779
DFT的执行效率是O(n^2),FFT为O(log2n),但是它对点数没有限制。
/*
离散傅立叶算法V1.0
含有:DFT,IDFT
made by xyt
2015/7/5
*/
#ifndef _DFT_H
#define _DFT_H
#include<iostream>
#include<math.h>
using namespace std;
#define PI 3.14159265354
struct complex{
double r,i;
};
complex multi(complex a,complex b){
complex tmp;
tmp.r=a.r*b.r-a.i*b.i;
tmp.i=a.r*b.i+a.i*b.r;
return tmp;
}
int fi(double in){
if((in-(int)in)>0.5) return (int)in+1;
else return (int)in;
}
/* 离散傅立叶正变换,输出[][2]数组实部在前,采样容量n可以任意正整数 */
void DFT(int *in,double **out,const int &n)
{
int i,j;
complex **W=new complex*[n];
for(i=0;i<n;i++){
W[i]=new complex[n];
}
complex *lis=new complex[(n-1)*(n-1)+1];
lis[0].r=1;lis[0].i=0;
lis[1].r=cos(2.0*PI/n);
lis[1].i=-1.0*sin(2.0*PI/n);
for(i=2;i<=(n-1)*(n-1);i++){
lis[i]=multi(lis[1],lis[i-1]);
}
for(i=0;i<n;i++){
for(j=0;j<n;j++){
W[i][j]=lis[i*j];
}
}
complex sum;
for(i=0;i<n;i++){
sum.r=0;sum.i=0;
for(j=0;j<n;j++){
sum.r+=in[j]*W[i][j].r;
sum.i+=in[j]*W[i][j].i;
}
out[i][0]=sum.r;
out[i][1]=sum.i;
}
for(i=0;i<n;i++) delete []W[i];
delete []W;
delete []lis;
}
/* 离散傅立叶逆变换 */
void IDFT(double **in,int *out,const int &n)
{
int i,j;
complex **W=new complex*[n];
for(i=0;i<n;i++){
W[i]=new complex[n];
}
complex *lis=new complex[(n-1)*(n-1)+1];
lis[0].r=1;lis[0].i=0;
lis[1].r=cos(2.0*PI/n);
lis[1].i=sin(2.0*PI/n);
for(i=2;i<=(n-1)*(n-1);i++){
lis[i]=multi(lis[1],lis[i-1]);
}
for(i=0;i<n;i++){
for(j=0;j<n;j++){
W[i][j]=lis[i*j];
}
}
complex sum;
for(i=0;i<n;i++){
sum.r=0;sum.i=0;
for(j=0;j<n;j++){
sum.r+=W[i][j].r*in[j][0]-W[i][j].i*in[j][1];
sum.i+=W[i][j].i*in[j][0]+W[i][j].r*in[j][1];
}
out[i]=fi(sum.r/n);
}
for(i=0;i<n;i++) delete []W[i];
delete []W;
delete []lis;
}
#endif