C语言做快速傅里叶变换和快速逆傅里叶变换
快速傅里叶变换(FFT)和快速逆傅里叶变换(IFFT)要求做傅里叶变换的数据点数只能是2的整数次幂,比如2,4,8,16,32,64,128,256,512,1024,2048,.......如果是2000个数据,那么用快速傅里叶变换(FFT)的结果就不对了,就需要使用对数据点数不限制的离散傅里叶变化(DFT)了,据说也可以使用补0的方法凑够2048个点来使用FFT。
这里只分享C语言做快速傅里叶变化的代码,以后再分享C语言做傅里叶变换的代码。C语言做快速傅里叶变化(FFT)的时间复杂度是O(nlogn),傅里叶变化(DFT)的时间复杂度是O(n*n).
以下是做快速傅里叶变换和快速逆傅里叶变换的代码(包含测试代码):
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct
{
double real;
double img;
}complex;
void initW(complex W[], int size_x, double PI);
void change(int size_x, complex x[]);
void add(complex a, complex b, complex *c);
void mul(complex a, complex b, complex *c);
void sub(complex a, complex b, complex *c);
void divi(complex a, complex b, complex *c);
void output(int size_x, complex x[]);
void fft(int size_x, complex x[], complex W[]);
void ifft(int size_x, complex x[], complex W[]);
#define N 1000
int main()
{
int size_x = 0;
complex x[N], W[N];
double PI;
PI = atan(1) * 4;
size_x = 16;
x[0].real = 1, x[0].img=0;
x[1].real = 2, x[1].img=0;
x[2].real = 3, x[2].img=0;
x[3].real = 4, x[3].img=0;
x[4].real = 5, x[4].img=0;
x[5].real = 6, x[5].img=0;
x[6].real = 7, x[6].img=0;
x[7].real = 8, x[7].img=0;
x[8].real = 9, x[8].img=0;
x[9].real = 10, x[9].img=0;
x[10].real = 11, x[10].img=0;
x[11].real = 12, x[11].img=0;
x[12].real = 13, x[12].img=0;
x[13].real = 14, x[13].img=0;
x[14].real = 15, x[14].img=0;
x[15].real = 16, x[15].img=0;
initW(W, size_x, PI);
fft(size_x, x, W);
output(size_x, x);
return 0;
}
void initW(complex W[], int size_x, double PI)
{
for(int i = 0; i < size_x; i++)
{
W[i].real = cos(2 * PI / size_x * i);
W[i].img = -1 * sin(2 * PI / size_x * i);
}
}
void change(int size_x, complex x[])
{
complex temp;
unsigned short i =0, j = 0, k = 0;
double t;
for(i = 0; i < size_x; i++)
{
k = i;
j = 0;
t = (log(size_x) / log(2));
while((t--)>0)
{
j = j << 1;
j |= (k&1);
k = k >> 1;
}
if(j > i)
{
temp = x[i];
x[i] = x[j];
x[j] = temp;
}
}
}
void add(complex a, complex b, complex *c)
{
c->real = a.real + b.real;
c->img = a.img + b.img;
}
void mul(complex a, complex b, complex *c)
{
c->real = a.real * b.real - a.img * b.img;
c->img = a.real * b.img + a.img * b.real;
}
void sub(complex a, complex b, complex *c)
{
c->real = a.real - b.real;
c->img = a.img - b.img;
}
void divi(complex a, complex b, complex *c)
{
c->real = (a.real * b.real + a.img * b.img) / (b.real * b.real + b.img * b.img);
c->img = (a.img * b.real - a.real * b.img) / (b.real * b.real + b.img * b.img);
}
void output(int size_x, complex x[])
{
int i;
printf("The result are as follows\n");
for(i = 0; i < size_x; i++)
{
printf("%.4f", x[i].real);
if(x[i].img >= 0.0001)
{
printf("+%.4fj\n",x[i].img);
}
else if(fabs(x[i].img) < 0.0001)
{
printf("\n");
}
else
{
printf("%.4fj\n",x[i].img);
}
}
}
void fft(int size_x, complex x[], complex W[])
{
int i = 0, j = 0, k = 0, l = 0;
complex up, down, product;
change(size_x, x);
for(i = 0; i < (int)(log10(size_x)/log10(2)); i++)
{
l = 1 << i;
for(j = 0; j < size_x; j += 2 * l)
{
for(k = 0; k < l; k++)
{
mul(x[j+k+l], W[size_x*k/2/l], &product);
add(x[j+k], product, &up);
sub(x[j+k], product, &down);
x[j+k] = up;
x[j+k+l] = down;
}
}
}
}
void ifft(int size_x, complex x[], complex W[])
{
int i = 0, j = 0, k = 0, l = size_x;
complex up, down;
for(i = 0; i < (int)(log(size_x) / log(2)); i++)
{
l/=2;
for(j = 0; j < size_x; j += 2 * l)
{
for(k = 0; k < l; k++)
{
add(x[j + k], x[j + k +l], &up);
up.real /= 2;
up.img /= 2;
sub(x[j+k], x[j+k+l], &down);
down.real /= 2;
down.img /=2;
divi(down, W[size_x * k / 2 / l], &down);
x[j + k] = up;
x[j + k + l] = down;
}
}
}
change(size_x, x);
}
分享一个在线做傅里叶变换的网站,可以验证自己做傅里叶变换的结果。http://www.yunsuan.info/signalprocessing/compute_fft_2d.htm