FFT和DFT完整算法C语言实现

#include <stdio.h> #include <math.h> #include <time.h> #define N 1024 //设置抽样点数 #define PI 3.1514926535897932384626433832795028841971 //定义圆周率 typedef struct //定义复数结构体变量 { double real; double imag; }complex; void c_plus(complex, complex, complex *); //复数加运算 void c_sub(complex, complex, complex *); //复数减运算 void c_mul(complex, complex, complex *); //复数乘运算 void Wn_i(int, int, complex *); //FFT旋转因子 void Wn_ik(int, int, int, complex *); //DFT旋转因子 int main() { complex f[N], A[N], a[N]; //f[N]为输出的快速傅里叶变换序列 A[N]为DFT变换输出序列 a[N]为初始序列 double x[N] = { 1, 2, 3, 4, 5, 6, 7, 8 }; //x[N]为进行抽样运算的序列 int LH, K, J, B, L, k, N1, P, M, K1; //L表示第L级蝶形 p旋转因子指数 B两序列间隔点数k 第k个序列 double T; clock_t begin, end; double cost1,cost2; M = (int)(log2(N)); LH = N / 2; J = LH; N1 = N - 1; for (int i = 0; i < N; i++) //为DFT运算提供初始序列 { a[i].real = x[i]; a[i].imag = 0; } printf("******************************* 级数为：%d *******************************\n", M); //....................................................................................................................... for (int I = 1; I < N1; I++) //定义倒序序列函数 { if (I < J) { T = x[I]; x[I] = x[J]; x[J] = T; } K = LH; if (J >= K) { do { J = J - K; K = K / 2; } while (J >= K); } J = J + K; } //....................................................................................................................... /* for (int i = 0; i < N; i++) //输出倒序后的序列 { printf("x[%d]为：%lf\n", i, x[i]); }*/ //....................................................................................................................... for (int i = 0; i < N; i++) //将序列赋给结构体函数进行运算 { f[i].real = x[i]; f[i].imag = 0; } /* for (int i = 0; i < N; i++) { printf("f[%d]为：%lf + j%lf\n", i, f[i].real, f[i].imag); }*/ //....................................................................................................................... begin = clock(); //开始记录时间 for (L = 1; L <= M; L++) //FFT运算 { B = (int)(pow(2, L - 1)); for (J = 0; J < B; J++) { P = (int)(J*pow(2, M - L)); for (k = J; k < N; k = (int)(k + pow(2, L))) { K1 = k + B; complex wn, t; Wn_i(N, P, &wn); c_mul(f[K1], wn, &t); //。。。。。。。。。。。。 c_sub(f[k], t, &(f[K1])); //蝶形运算 c_plus(f[k], t, &(f[k])); //。。。。。。。。。。。。 } } } end = clock(); //结束记录时间 cost1 = (double)(end - begin) / CLOCKS_PER_SEC; //....................................................................................................................... printf("*****************************快速傅里叶变换输出*****************************\n"); for (int i = 0; i < N; i++) //快速傅里叶变换输出 { printf("f[%d]为： %lf + j(%lf)\n", i, f[i].real, f[i].imag); } //....................................................................................................................... begin = clock(); //开始记录时间 for (int k = 0; k < N; k++) //DFT运算 { complex sum = { 0, 0 }; for (int n = 0; n < N; n++) { complex wn, t; Wn_ik(N, n, k, &wn); c_mul(a[n], wn, &t); c_plus(sum, t, &sum); } A[k].real = sum.real; A[k].imag = sum.imag; } end = clock(); //结束记录时间 cost2 = (double)(end - begin) / CLOCKS_PER_SEC; //....................................................................................................................... printf("*****************************DFT变换输出*****************************\n"); for (int i = 0; i < N; i++) //DFT变换输出 { printf("A[%d]为： %lf + j(%lf)\n", i, A[i].real, A[i].imag); } //....................................................................................................................... printf("FFT算法一共花费时间为：%lf秒\n", cost1); //CLOCKS_PER_SEC 单位1000ms printf("DFT算法一共花费时间为：%lf秒\n", cost2); //CLOCKS_PER_SEC 单位1000ms return 0; } void c_plus(complex a, complex b, complex *c) //复数加法 { c->real = a.real + b.real; c->imag = a.imag + b.imag; } void c_sub(complex a, complex b, complex *c) //复数减法 { c->real = a.real - b.real; c->imag = a.imag - b.imag; } void c_mul(complex a, complex b, complex *c) //复数乘法 { c->real = a.real*b.real - a.imag*b.imag; c->imag = a.real*b.imag + a.imag*b.real; } void Wn_i(int n1, int i, complex *Wn) //定义FFT旋转因子 { Wn->real = cos(2 * PI*i / n1); Wn->imag = -sin(2 * PI*i / n1); } void Wn_ik(int n1, int i, int k, complex *Wn) //定义DFT旋转因子 { Wn->real = cos(2 * PI*i*k / n1); Wn->imag = -sin(2 * PI*i*k / n1); }