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

Question

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

最初只是 brute-force 方法。首先，我让程序打开一个数据文件（振幅）并将数据放入一个数组中（只有一个，因为我将自己限制为 real-valued 输入）。

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

这是我的代码，去掉了花里胡哨的东西：

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#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;
}

程序编译、完成，但我得到的不是功率（频率）图上的几个尖峰，而是：。

单个频率或不同频率给出完全相同的inverted-bathtub-curve。

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

dftreal

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

Answer 1

警告：我对此有点生疏

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

此外，我将从更简单的数据输入开始。

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

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

我还发现除法 i / linecount 有时会被截断，所以我为 illustrate/correct 创建了一个 FRAG 宏。如果没有这个，输入 cos(2 * pi * t) 的尖峰最终会出现在 outreal[0]（DC 部分？）而不是 outreal[1]

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

无论如何，这是代码[请原谅不必要的样式清理]：

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

//#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 2

好消息是您的 dftreal 实施没有任何问题。

这个问题，如果可以这么说的话，是你正在使用的测试波形包含相对于你的采样率来说频率非常低的频率分量 linecount。相应地，DFT 显示能量集中在前几个 bin 中，一些频谱泄漏到更高的 bin 中，因为波形频率分量不是采样率的精确倍数。

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

//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 3

您提到“...转换看起来不对，...”

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

您要输入的原始振幅数组有多大？如果你post这里的样本数据，我愿意运行通过我自己的程序，看看我的程序输出的结果是否和你的一样。

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

Discrete Fourier transform giving incorrect results

c

dft