C语言做离散傅里叶变换和逆变换
快速傅里叶换要求做傅里叶变换的数据点数是2的整数次幂，即点数是2,4,8,16,32,64,128,256,512,1024,…。实际使用中很多数据点数并不是2的整数次幂，这时可以利用数据变化规律将其拓展为2的整数次幂个数据点，当然，还可以直接使用离散傅里叶变化，离散傅里叶变换对数据点数没有要求，本文分享自己使用的C语言做离散傅里叶变换和逆变换的代码。需要注意的是，快速傅里叶变化的时间复杂度是O(nlogn),而离散傅里叶变换的时间复杂度是O(n*n),当数据点数在1000以上时，离散傅里叶变化花费的时间明显大于快速傅里叶变换。
#include <stdio.h>#include <stdlib.h>#include <math.h>#define PI 3.1415926typedef struct{double real,imag;} complex; //定义复数结构void DFT_Calculate_Point(complex Input_Squence[], complex Result_Point[], int k, int len){int n = 0;complex Sum_Point;complex One_Point[len];Sum_Point.real = 0;Sum_Point.imag = 0;for (n = 0; n < len; n++){One_Point[n].real = (Input_Squence[n].real) * cos(2 * PI * k * n / len)+ Input_Squence[n].imag * sin(2 * PI * k * n / len); //复数的实部One_Point[n].imag = Input_Squence[n].imag * cos(2 * PI * k * n / len)- (Input_Squence[n].real) * sin(2 * PI * k * n / len); //复数的虚部Sum_Point.real += One_Point[n].real; //对实部求和Sum_Point.imag += One_Point[n].imag; //对虚部求和}Result_Point[k].real = Sum_Point.real;Result_Point[k].imag = Sum_Point.imag;}void DFT_Calculate(complex Input_Squence[], complex Result_Point[], int len){int i = 0;for (i = 0; i < len; i++){DFT_Calculate_Point(Input_Squence, Result_Point, i, 14);printf("%lf %lf\n", Result_Point[i].real, Result_Point[i].imag);}}void IDFT_Calculate_Point(complex Input_Squence[], complex Result_Point[], int k, int len){int n = 0;complex Sum_Point;complex One_Point[len];Sum_Point.real = 0;Sum_Point.imag = 0;for (n = 0; n < len; n++){One_Point[n].real = (Input_Squence[n].real) * cos(2 * PI * k * n / len)- Input_Squence[n].imag * sin(2 * PI * k * n / len); //复数的实部One_Point[n].imag = Input_Squence[n].imag * cos(2 * PI * k * n / len)+ (Input_Squence[n].real) * sin(2 * PI * k * n / len); //复数的虚部Sum_Point.real += One_Point[n].real; //对实部求和Sum_Point.imag += One_Point[n].imag; //对虚部求和}Result_Point[k].real = (Sum_Point.real) / len;Result_Point[k].imag = (Sum_Point.imag) / len;}void IDFT_Calculate(complex Input_Squence[], complex Result_Point[], int len){int i = 0;for (i = 0; i < len; i++){IDFT_Calculate_Point(Input_Squence, Result_Point, i, 14);printf("%lf %lf\n", Result_Point[i].real, Result_Point[i].imag);}}int main(int argc, char *argv[]){int i = 0;complex Input_Squence[14], Result_Point[100];// Input_Squence[0].real=105, Input_Squence[0].imag=0;// Input_Squence[1].real=-7, Input_Squence[1].imag=30.669;// Input_Squence[2].real=-7, Input_Squence[2].imag=14.53565;// Input_Squence[3].real=-7, Input_Squence[3].imag=8.77772;// Input_Squence[4].real=-7, Input_Squence[4].imag=5.58231;// Input_Squence[5].real=-7, Input_Squence[5].imag=3.37102;// Input_Squence[6].real=-7, Input_Squence[6].imag=1.59770;// Input_Squence[7].real=-7, Input_Squence[7].imag=0;// Input_Squence[8].real=-7, Input_Squence[8].imag=-1.59770;// Input_Squence[9].real=-7, Input_Squence[9].imag=-3.37102;// Input_Squence[10].real=-7, Input_Squence[10].imag=-5.58231;// Input_Squence[11].real=-7, Input_Squence[11].imag=-8.77772;// Input_Squence[12].real=-7, Input_Squence[12].imag=-14.53565;// Input_Squence[13].real=-7, Input_Squence[13].imag=-30.669;for (i = 0; i < 14; i++) //产生输入序列{Input_Squence[i].real = i + 1;Input_Squence[i].imag = 0;}// IDFT_Calculate(Input_Squence, Result_Point, 14); //进行DFT计算DFT_Calculate(Input_Squence, Result_Point, 14);return 0;}
上面给出了做离散傅里叶变换和逆变换的示例代码，包含了测试代码。分享一个在线做傅里叶变换的网站，可以验证自己做傅里叶变换的结果。http://www.yunsuan.info/signalprocessing/compute_fft_2d.html
