离散傅里叶变换C++代码

/*
* myfft.h
*/

#ifndef __MYFFT_H__
#define __MYFFT_H__

#include <windows.h>

typedef struct _my_complex
{
double r; //复数实部
double i; //复数虚部

_my_complex(){}
_my_complex(double _r, double _i){r = _r; i = _i;}
~_my_complex(){}

_my_complex operator + (const _my_complex & c) //重载加号运算符
{
return _my_complex(r + c.r, i + c.i);
}

_my_complex operator - (const _my_complex & c) //重载减号运算符
{
return _my_complex(r - c.r, i - c.i);
}

_my_complex operator * (const _my_complex & c) //重载乘号运算符
{
return _my_complex(r * c.r - i * c.i, r * c.i + i * c.r);
}
}my_complex;

int myfft_by_define_c2c_1d(my_complex * in, int n1, my_complex * out, int n2, int line_width = 1, int is_inv = 0); //按照定义进行一维复数离散傅里叶变换
int myfft_by_define_c2c_2d(my_complex * in, int m, int n, int is_inv = 0); //按照定义进行二维复数离散傅里叶变换

int myfft_c2c_1d_base2(my_complex * in, int n, int line_width = 1, int is_inv = 0); //进行一维复数离散快速傅里叶变换(基2)

int print_complex_array_1d(my_complex * in, int n);
int print_complex_array_2d(my_complex * in, int m, int n);

#endif //__MYFFT_H__

/*
* myfft.cpp
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "myfft.h"

#define PI 3.1415926

/*------------按照定义进行一维复数离散傅里叶变换----------------------------------------
* 原理：
*     一维离散傅里叶变换公式为： X[k] = ∑x
*e^(-2πi * nk/N)，其中n=0,1,2,...,N-1
*     一维离散傅里叶逆变换公式为： x
= 1/N * ∑X[k]*e^(2πi * nk/N)，其中k=0,1,2,...,N-1
*     欧拉公式：e^ix = cos(x) + isin(x)
*
* 参数1：in            输入的一维复数数组
* 参数2：n1            输入的一维复数数组的长度
* 参数3：out           输出的一维复数数组
* 参数4：n2            输出的一维复数数组的长度(必须和n1相等)
* 参数5：line_width    数组下标的间隔(默认为1)，对于行变换等于1，对于列变换等于一行的宽度
* 参数6：is_inv        是否是逆变换(默认不是)
*
* 过客 & 386520874@qq.com & 2016.01.25
*/
int myfft_by_define_c2c_1d(my_complex * in, int n1, my_complex * out, int n2, int line_width/*=1*/, int is_inv/*=0*/)
{
if(n1 != n2){return 0;}

double _2PI_N = 2.0 * PI / n1;

if(is_inv == 0) //正变换
{
for(int m = 0; m < n1; m++)
{
int M = m * line_width;
for(int k = 0; k < n1; k++ )
{
int K = k * line_width;
out[M].r += in[K].r * cos(_2PI_N * m * k) + in[K].i * sin(_2PI_N * m * k); //计算实部
out[M].i += -in[K].r * sin(_2PI_N * m * k) + in[K].i * cos(_2PI_N * m * k); //计算虚部
}
}
}else if(is_inv == 1) //逆变换
{
for(int m = 0; m < n1; m++)
{
int M = m * line_width;
for(int k = 0; k < n1; k++ )
{
int K = k * line_width;
out[M].r += in[K].r * cos(_2PI_N * m * k) - in[K].i * sin(_2PI_N * m * k); //计算实部
out[M].i += in[K].r * sin(_2PI_N * m * k) + in[K].i * cos(_2PI_N * m * k); //计算虚部
}
out[M].r /= n1; //归一化
out[M].i /= n1; //归一化
}
}

return 1;
}

//按照定义进行二维复数离散傅里叶变换
int myfft_by_define_c2c_2d(my_complex * in, int m, int n, int is_inv/*=0*/)
{
my_complex * temp = new my_complex[m * n];
memset(temp, 0, sizeof(my_complex) * m * n);

for(int y = 0; y < n; y++) //先行变换
{
myfft_by_define_c2c_1d(in + y * n, n, temp + y * n, n, 1, is_inv);
}

memset(in, 0, sizeof(my_complex) * m * n);
for(int x = 0; x < m; x++) //再列变换
{
myfft_by_define_c2c_1d(temp + x, m, in + x, m, n, is_inv);
}

delete [] temp; temp = NULL;

return 1;
}

/*------------进行一维复数离散快速傅里叶变换(基2)----------------------------------------
* 原理：
*     对于FFT算法，一般是基2型的，即要求被变换的数组长度必须是2的m幂次方，
*     这种要求过于苛刻，当然相比于按照定义进行计算来说，速度是很快的。此外，
*     还有基4的算法，以及混合基算法等，但都对数组长度N有严格的限制。
*     所以想自己开发一个适合任意N的FFT算法，还是比较难的，目前还是调用FFTW库暂时凑合着用吧。
*
* 参数1：in            输入的一维复数数组
* 参数2：n             输入/出的一维复数数组的长度
* 参数3：line_width    数组下标的间隔(默认为1)，对于行变换等于1，对于列变换等于一行的宽度
* 参数4：is_inv        是否是逆变换(默认不是)
*
* 过客 & 386520874@qq.com & 2016.01.25
*/
int myfft_c2c_1d_base2(my_complex * in, int n, int line_width/*=1*/, int is_inv/*=0*/)
{
int ttt = n;
int L = 1; //蝶形运算的层数
while(ttt = ttt >> 1){L++;}

my_complex * X = in;

int N = n;
int len = 1 << L; //将输入的一维数组对齐为 N=2^L
if(n * 2 == len) //刚好是2的幂
{
L -= 1;
}else if(n < len) //说明n不是2的幂
{
//		printf("Erorr: n is not 2^m format.\n"); //非基2算法还有待开发
//		return 0;

N = len;
X = new my_complex
; //扩充数组大小
int aaa = sizeof(my_complex);
memset(X, 0, sizeof(my_complex) * N); //全部初始化为0
memcpy(X, in, sizeof(my_complex) * n); //复制数据
}

//--------数组倒位序 二进制(n0n1n2) => (n2n0n1)-------------------
//例如： 3的二进制为 011，对应的倒位序为 110，即为6
my_complex * XX = new my_complex
;
for(int i = 0; i < N; i++)
{
int p = 0;
for(int f = 0; f < L; f++)
{
if(i & (1 << f))
{
p += 1 << (L - f - 1);
}
}
XX[i] = X[p];
}
memcpy(X, XX, sizeof(my_complex) * N); //复制数据
delete [] XX; XX = NULL;

//---------旋转因子-------------------
my_complex * W = new my_complex[N / 2];
memset(W, 0, sizeof(my_complex) * N / 2); //全部初始化为0
for(int i = 0; i < N / 2; i++)
{
W[i] = my_complex(cos(-2.0 * PI * i / N), sin(-2.0 * PI * i / N));
}

//---------FFT算法(基2)----------------
for(int f = 0; f < L; f++) //蝶形运算层数
{
int G = 1 << (L - f);
int num = 1 << f;

for(int u = 0; u < num; u++) //每组的元素个数
{
for(int g = 0; g < G; g += 2) //每层的组数
{
my_complex X1 = X[g * num + u];
my_complex X2 = X[(g + 1) * num + u];

X[g * num + u]       = X1 + W[u * (1 << (L - f - 1))] * X2; //蝶形运算
X[(g + 1) * num + u] = X1 - W[u * (1 << (L - f - 1))] * X2; //蝶形运算
}
}
}

if(X != NULL && X != in) //说明另外申请了内存空间
{
memcpy(in, X, sizeof(my_complex) * n);
if(X){delete [] X; X = NULL;}
}
if(W != NULL){delete [] W; W = NULL;}

return 1;
}

int print_complex_array_1d(my_complex * in, int n)
{
for(int i = 0; i < n; i++)
{
printf("%d:%f+i%f ", i, in[i].r, in[i].i);
}
printf("\n");

return 1;
}

int print_complex_array_2d(my_complex * in, int m, int n)
{
for(int y = 0; y < n; y++)
{
for(int x = 0; x < m; x++)
{
printf("%d_%d:%f+i%f ", y, x, (*(in + y * m + x)).r, (*(in + y * m + x)).i);
}
printf("\n");
}
printf("\n");

return 1;
}

/*
* main.cpp
*/

#include <stdio.h>
#include <stdlib.h>
#include "myfft.h"

int main(int argc, char * argv[]) //DFT测试 一维
{
int n = 8;
int n1 = n;
int n2 = n;
my_complex * in = new my_complex[n1];
my_complex * out = new my_complex[n2];

double data[] = {4, 3, 2, 6, 7, 8, 9, 0};

for(int i = 0; i < n; i++)
{
in[i].r = data[i]; in[i].i = 0;
out[i].r = 0; out[i].i = 0;
}

print_complex_array_1d(in, n1);

//----------正变换----------------------
myfft_by_define_c2c_1d(in, n1, out, n2, 1, 0); //按照定义进行一维复数离散傅里叶变换

print_complex_array_1d(out, n2);

//----------逆变换----------------------
memset(in, 0, sizeof(my_complex) * n);

myfft_by_define_c2c_1d(out, n2, in, n1, 1, 1); //按照定义进行一维复数离散傅里叶变换

print_complex_array_1d(in, n1);

if(in){delete [] in; in = NULL;}
if(out){delete [] out; out = NULL;}

system("pause");

return 0;
}
/*
0:4.000000+i0.000000 1:3.000000+i0.000000 2:2.000000+i0.000000 3:6.000000+i0.000000 4:7.000000+i0.000000 5:8.000000+i0.000000 6:9.000000+i0.000000 7:0.000000+i0.000000
0:39.000000+i0.000000 1:-10.778175+i6.292892 2:0.000001+i-5.000001 3:4.778176+i-7.707106 4:5.000000+i0.000001 5:4.778172+i7.707108 6:-0.000002+i4.999998 7:-10.778169+i-6.292899
0:4.000000+i-0.000001 1:3.000000+i-0.000001 2:2.000000+i-0.000000 3:6.000000+i-0.000000 4:7.000000+i0.000000 5:8.000000+i0.000000 6:9.000000+i0.000001 7:0.000000+i0.000001
请按任意键继续. . .
*/

int main2(int argc, char * argv[]) //DFT测试 二维
{
int m = 3;
int n = 3;
my_complex * in = new my_complex[m * n];

in[0] = my_complex(0, 0); in[1] = my_complex(2, 0); in[2] = my_complex(4, 0);
in[3] = my_complex(6, 0); in[4] = my_complex(1, 0); in[5] = my_complex(3, 0);
in[6] = my_complex(5, 0); in[7] = my_complex(7, 0); in[8] = my_complex(4, 0);

print_complex_array_2d(in, m, n);

//----------正变换----------------------
myfft_by_define_c2c_2d(in, m, n, 0); //按照定义进行二维复数离散傅里叶变换

print_complex_array_2d(in, m, n);

//----------逆变换----------------------
myfft_by_define_c2c_2d(in, m, n, 1); //按照定义进行二维复数离散傅里叶变换

print_complex_array_2d(in, m, n);

if(in){delete [] in; in = NULL;}

system("pause");

return 0;
}
/*
0_0:0.000000+i0.000000 0_1:2.000000+i0.000000 0_2:4.000000+i0.000000
1_0:6.000000+i0.000000 1_1:1.000000+i0.000000 1_2:3.000000+i0.000000
2_0:5.000000+i0.000000 2_1:7.000000+i0.000000 2_2:4.000000+i0.000000

0_0:32.000000+i0.000000 0_1:0.500000+i0.866025 0_2:0.500001+i-0.866027
1_0:-7.000001+i5.196152 1_1:-1.000000+i-1.732051 1_2:-8.499999+i-6.062178
2_0:-6.999999+i-5.196154 2_1:-8.500001+i6.062177 2_2:-1.000000+i1.732051

0_0:0.000000+i-0.000000 0_1:2.000000+i-0.000000 0_2:4.000000+i-0.000000
1_0:6.000000+i-0.000000 1_1:1.000000+i-0.000000 1_2:3.000000+i0.000000
2_0:5.000000+i-0.000000 2_1:7.000000+i0.000000 2_2:4.000000+i0.000001

请按任意键继续. . .
*/

int main3(int argc, char * argv[]) //FFT测试 一维(基2)
{
int n = 8;
my_complex * in = new my_complex
;
double data[] = {4, 3, 2, 6, 7, 8, 9, 0};

for(int i = 0; i < n; i++)
{
in[i].r = data[i];
in[i].i = 0;
}

print_complex_array_1d(in, n);

//----------正变换----------------------
myfft_c2c_1d_base2(in, n, 1, 0); //进行一维复数离散快速傅里叶变换(基2)

print_complex_array_1d(in, n);

if(in){delete [] in; in = NULL;}

system("pause");

return 0;
}
/*
0:4.000000+i0.000000 1:3.000000+i0.000000 2:2.000000+i0.000000 3:6.000000+i0.000000 4:7.000000+i0.000000 5:8.000000+i0.000000 6:9.000000+i0.000000 7:0.000000+i0.000000
0:39.000000+i0.000000 1:-10.778175+i6.292893 2:0.000000+i-5.000000 3:4.778175+i-7.707106 4:5.000000+i0.000000 5:4.778174+i7.707107 6:-0.000000+i5.000000 7:-10.778175+i-6.292894
请按任意键继续. . .
*/