离散傅立叶变换给出不正确的结果

我试图在C中进行离散傅立叶变换。

最初只是蛮力方法。 首先，我让程序打开一个数据文件（幅度）并将数据放入一个数组（只有一个，因为我限制自己使用实值输入）。

但是转换看起来不对，所以我尝试生成一个简单的波函数并检查它是否正确转换。

这是我的代码，剥夺了花里胡哨的声音：

#include  #include  #include  #include  #define M_PI 3.14159265358979323846 //the test wavefunction double theoretical(double t) { double a = sin(M_PI * t) + 2 * sin(2 * M_PI * t) + 4 * sin(4 * M_PI * t); return a; } //------------------------------------------------------------------------- void dftreal(double inreal[], double outreal[], double outimag[], int linecount) { int n, k; for (k = 0; k < linecount; k++) { double sumreal = 0; double sumimag = 0; for (n = 0; n < linecount; n++) { double angle = 2 * M_PI * n * ( k / (double) linecount); sumreal += inreal[n] * cos(angle); sumimag += inreal[n] * sin(angle); } outreal[k] = sumreal; outimag[k] = sumimag; } } //========================================================================= int main(void) { int linecount = 44100; //creates all necessary arrays double inreal[linecount], outreal[linecount], outimag[linecount], p[linecount]; FILE *fout = fopen("Output.txt", "w"); for (int i = 0 ; i < linecount ; ++i) { inreal[i] = theoretical( i / (double) linecount); } //actually computes the transform dftreal(inreal, outreal, outimag, linecount); for (int i = 0 ; i < linecount ; ++i) { p[i] = 2*(outreal[i] * outreal[i] + outimag[i] * outimag[i]); fprintf(fout, "%f %f \n", (i / (double) linecount), p[i]); } fclose(fout); printf("\nEnd of program"); getchar(); return 0; } 

程序编译，完成，但在功率（频率）图上没有几个尖峰，我得到这个： 。

单个频率或不同频率给出完全相同的倒置浴盆曲线。

我检查了几个关于DFT的消息来源，我仍然不知道出了什么问题，这个函数似乎没有任何明显的错误：

 dftreal 

本身。 我想就可能导致问题的原因寻求帮助。 我在Windows 7上使用MinGW编译器。谢谢！

Answer 1

好消息是你的dftreal实现没有任何问题。

如果可以称之为问题，那么您使用的测试波形包括频率分量，这些频率分量的频率相对于采样率linecount数非常低。 相应地，DFT表明能量集中在前几个箱中，一些谱泄漏到更高的箱中，因为波形频率分量不是采样率的精确倍数。

如果通过使频率相对于采样频率增加测试波形频率，例如：

 //the test wavefunction double theoretical(double t) { double f = 0.1*44100; double a = sin(2 * M_PI * f * t) + 2 * sin(4 * M_PI * f * t) + 4 * sin(8 * M_PI * f * t); return a; } 

你应该得到一个情节，如：

Answer 2

警告：我对此有点生疏

我记得它，cos / real部分产生频谱，sin / imag部分产生功率谱。 所以，我认为你想要打印的是频谱（它只是outreal[i] ）而不是你所做的（即添加outreal和outimag是错误的）。 你可以用不同的颜色绘制，但我会简单地开始。

另外，我从简单的数据输入开始。

我改变了theoretical来做一个单余弦波[以及你原来的]。 这是可预测的，你应该得到一个尖峰，你这样做，所以我认为dftreal是正确的。

我还添加了一些命令行选项，以便您可以进行实验。

我还发现划分i / linecount有时会被截断，所以我创建了一个FRAG宏来说明/纠正这个问题。 没有这个， cos(2 * pi * t)输入的尖峰最终出现在outreal[0] （DC部分？）而不是outreal[1]

此外，出于测试目的，您可以使用更少的样本（例如1000）。 这也可能有助于横向展开你的情节，以便更好地看到事物。

无论如何，这是代码[请原谅无偿风格的清理]：

 #include  #include  #include  #include  #include  //#define M_PI 3.14159265358979323846 #define M_2PI (M_PI * 2.0) int opt_A; // algorithm to use int linecount; // number of samples int opt_a; // use absolute value on output int opt_s; // use squared value on output int opt_i; // use integer output index int opt_j; // output imaginary part int opt_c; // .csv compatibility time_t todlast; // the first was your original and would truncate #if 0 #define FRAG(_i) ((_i) / linecount) #else #define FRAG(_i) ((double) (_i) / linecount) #endif void pgr(int k,int n,int count) { time_t todnow = time(NULL); if ((todnow - todlast) >= 1) { todlast = todnow; printf("\r%d %d ",count,k); fflush(stdout); } } // the test wavefunction double theoretical(double t) { double a; switch (opt_A) { case 0: a = 0.0; a += sin(M_PI * t); a += 2.0 * sin(2.0 * M_PI * t); a += 4.0 * sin(4.0 * M_PI * t); break; default: a = cos(M_2PI * t * opt_A); break; } return a; } //------------------------------------------------------------------------- void dftreal(double inreal[], double outreal[], double outimag[], int linecount) { int n; int k; for (k = 0; k < linecount; k++) { double sumreal = 0; double sumimag = 0; double omega_k = (M_2PI * k) / (double) linecount; for (n = 0; n < linecount; n++) { // omega k: #if 0 double angle = M_2PI * n * FRAG(k); #else double angle = omega_k * n; #endif sumreal += inreal[n] * cos(angle); sumimag += inreal[n] * sin(angle); pgr(k,n,linecount); } outreal[k] = sumreal; outimag[k] = sumimag; } } //========================================================================= int main(int argc,char **argv) { char *cp; --argc; ++argv; for (; argc > 0; --argc, ++argv) { cp = *argv; if (*cp != '-') break; switch (cp[1]) { case 'a': // absolute value opt_a = ! opt_a; break; case 's': // square the output value opt_s = ! opt_s; break; case 'c': // .csv output opt_c = ! opt_c; break; case 'i': // integer output opt_i = ! opt_i; break; case 'j': // imaginary output opt_j = ! opt_j; break; case 'A': opt_A = atoi(cp + 2); break; case 'N': linecount = atoi(cp + 2); break; } } if (linecount <= 0) linecount = 44100 / 2; // creates all necessary arrays double inreal[linecount]; double outreal[linecount]; double outimag[linecount]; #if 0 double p[linecount]; #endif FILE *fout = fopen("Output.txt", "w"); for (int i = 0; i < linecount; ++i) inreal[i] = theoretical(FRAG(i)); todlast = time(NULL); // actually computes the transform dftreal(inreal, outreal, outimag, linecount); for (int i = 0; i < linecount; ++i) { #if 0 p[i] = 2 * (outreal[i] * outreal[i] + outimag[i] * outimag[i]); fprintf(fout, "%f %f\n", (i / (double) linecount), p[i]); #endif #if 0 fprintf(fout, "%f %f %f\n", (i / (double) linecount), outreal[i] * outreal[i], outimag[i] * outimag[i]); #endif #if 0 fprintf(fout, "%f %f\n", (i / (double) linecount),outreal[i] * outreal[i]); #endif if (opt_a) { outreal[i] = fabs(outreal[i]); outimag[i] = fabs(outimag[i]); } if (opt_s) { outreal[i] *= outreal[i]; outimag[i] *= outimag[i]; } if (! opt_c) { if (opt_i) fprintf(fout, "%d ",i); else fprintf(fout, "%f ",FRAG(i)); } fprintf(fout, "%f",outreal[i]); if (opt_j) fprintf(fout, "%f",outimag[i]); fprintf(fout, "\n"); } fclose(fout); printf("\nEnd of program\n"); //getchar(); return 0; } 

Answer 3

你提到“……变换看起来不对，……”

您是否将结果与其他程序或软件包进行比较，以确认结果确实是错误的？ 我最近用C ++和JavaScript编写了一个DFT，并将输出与MATLAB，Maple和MathCAD的结果进行了比较。 有时，差异只是缩放或格式化的问题。

您想要输入的原始振幅数量有多大？ 如果你在这里发布样本数据，我愿意通过我自己的程序运行它，看看我的程序是否输出与你相同的结果。