论坛 › 技巧探讨 › 【翻译】最简单的DFT算法的C代码

0xAA55 发表于 2018-10-24 10:10:30

【翻译】最简单的DFT算法的C代码

翻译原文：https://batchloaf.wordpress.com/2013/12/07/simple-dft-in-c/
原文作者：Ted Burke

我最近花了一些时间研究离散傅立叶变换（DFT），以及快速傅里叶变换（FFT）能进行更高效的DFT变换计算。

DFT定义如下。指定离散时间序列x［n］，其中n = 0,1,2，...，N-1


De Moivre的定理表明：


因此，我们可以将DFT求和重写为X［k］的实部和虚部的表达式，如下所示（假设x［n］为实数）。


其中：


且


作为我学习过程的一部分，我一直在尝试编写C和MATLAB / Octave中DFT和FFT的不同实现。这让我感觉很爽，因为我写了个非常简单的DFT暴力求解算法。//
// dft.c - 简单暴力求解DFT
// 作者：Ted Burke
// 最后更新：7-12-2013
//
// 编译：（用mingw或者cygwin，或者直接上Linux。在Linux可以不用加上最后的“.exe”后缀）
//    gcc dft.c -o dft.exe
//
// 运行：
//    dft.exe
// （Linux上则用“./dft”来运行）
//

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

#define N 16
#define PI 3.1415926535897932384626433832795
#define PI2 (PI * 2)

int main()
{
        // 时域和频域数据的数组
        int n, k;             // 时域和频域的索引
        float x;         // 离散时间信号x
        float Xre, Xim; // DFT变换后的x（实部和虚部）
        float P;         // x的功率图
       
        // 生成在(-1,+1)区间内的随机离散时间信号x
        srand(time(0));
        for (n=0 ; n<N ; ++n) x = ((2.0 * rand()) / RAND_MAX) - 1.0;
       
        // 用暴力算法计算x的DFT变换
        for (k=0 ; k<N ; ++k)
        {
                // X的实部
                Xre = 0;
                for (n=0 ; n<N ; ++n) Xre += x * cos(n * k * PI2 / N);
               
                // X的虚部
                Xim = 0;
                for (n=0 ; n<N ; ++n) Xim -= x * sin(n * k * PI2 / N);
               
                // 第k个频率点的功率
                P = Xre*Xre + Xim*Xim;
        }
       
        // 输出结果到MATLAB / Octave的M文件来绘制图表
        FILE *f = fopen("dftplots.m", "w");
        fprintf(f, "n = ;\n", N-1);
        fprintf(f, "x = [ ");
        for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x);
        fprintf(f, "];\n");
        fprintf(f, "Xre = [ ");
        for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xre);
        fprintf(f, "];\n");
        fprintf(f, "Xim = [ ");
        for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xim);
        fprintf(f, "];\n");
        fprintf(f, "P = [ ");
        for (k=0 ; k<N ; ++k) fprintf(f, "%f ", P);
        fprintf(f, "];\n");
        fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
        fprintf(f, "xlim()\n", N-1);
        fprintf(f, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
        fprintf(f, "xlim()\n", N-1);
        fprintf(f, "subplot(3,1,3)\nstem(n,P)\n");
        fprintf(f, "xlim()\n", N-1);
        fclose(f);
       
        // 正常退出
        return 0;
}用下图的方式来编译并运行上面的代码：
https://batchloaf.files.wordpress.com/2013/12/dftconsole.png
当你运行编译好的程序dft.exe后，它会生成一个名为dftplots.m的M文件，然后可以在MATLAB或Octave中运行来生成结果图。这是dft.exe生成的M文件：n = ;
x = [ 0.672957 -0.453061 -0.835088 0.980334 0.972232 0.640295 0.791619 -0.042803 0.282745 0.153629 0.939992 0.588169 0.189058 0.461301 -0.667901 -0.314791 ];
Xre = [ 4.358686 -2.627209 -2.558252 2.144204 1.888348 3.210599 2.147089 -1.166725 0.332542 -1.166852 2.146972 3.210663 1.888432 2.144283 -2.558007 -2.627225 ];
Xim = [ 0.000000 -0.959613 -0.352430 1.534066 0.408726 -0.478590 -0.390162 0.160532 0.000022 -0.160326 0.389984 0.478336 -0.408804 -1.534269 0.352723 0.959997 ];
P = [ 18.998148 7.823084 6.668859 6.950971 3.732913 10.536995 4.762218 1.387017 0.110584 1.387247 4.761577 10.537162 3.733295 6.951930 6.667812 7.823905 ];
subplot(3,1,1)
plot(n,x)
xlim()
subplot(3,1,2)
plot(n,Xre,n,Xim)
xlim()
subplot(3,1,3)
stem(n,P)
xlim()我在Octave中打开了dftplots.m，如下所示：
https://batchloaf.files.wordpress.com/2013/12/dftoctave.png
它产生了如下图表：
https://batchloaf.files.wordpress.com/2013/12/dftplots1.png
最上面的是原始随机离散时间序列x［n］。第二个图是X［k］的实部和虚部。注意X［k］的对称性，正如我们对实数x［n］所预想的那样。最后那个图是P［k］也就是x［n］的功率图。P［k］的每个值是对应的复数值X［k］的大小的平方。

作为上述示例的后续，我稍微修改了程序来进行更快的DFT变换，并转换更长（64个样本）的数据帧。我还稍微修改了随机时域信号生成过程以添加主频分量（在k = 5.7，即在区间5和6之间），用于测试算法的正确性。

这是修改过的C代码：//
// dft2.c - 基本，但效率更高的DFT
// 作者：Ted Burke
// 最后更新：10-12-2013
//
// 编译：（用mingw或者cygwin，或者直接上Linux。在Linux可以不用加上最后的“.exe”后缀）
//    gcc dft2.c -o dft2.exe
//
// 运行：
//    dft2.exe
// （Linux上则用“./dft2”来运行）
//

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

// 假设N大于4并且是2的幂
#define N 64
#define PI 3.1415926535897932384626433832795
#define PI2 (PI * 2)

// 旋转因子（64次单位根）
// https://zh.wikipedia.org/wiki/%E6%97%8B%E8%BD%89%E5%9B%A0%E5%AD%90
// https://zh.wikipedia.org/wiki/%E5%8D%95%E4%BD%8D%E6%A0%B9
const float W[] = {
       1.00000, 0.99518, 0.98079, 0.95694, 0.92388, 0.88192, 0.83147, 0.77301,
       0.70711, 0.63439, 0.55557, 0.47139, 0.38268, 0.29028, 0.19509, 0.09801,
        -0.00000,-0.09802,-0.19509,-0.29029,-0.38269,-0.47140,-0.55557,-0.63440,
        -0.70711,-0.77301,-0.83147,-0.88192,-0.92388,-0.95694,-0.98079,-0.99519,
        -1.00000,-0.99518,-0.98078,-0.95694,-0.92388,-0.88192,-0.83146,-0.77300,
        -0.70710,-0.63439,-0.55556,-0.47139,-0.38267,-0.29027,-0.19508,-0.09801,
       0.00001, 0.09803, 0.19510, 0.29030, 0.38269, 0.47141, 0.55558, 0.63440,
       0.70712, 0.77302, 0.83148, 0.88193, 0.92388, 0.95694, 0.98079, 0.99519
};

int main()
{
        // 时域和频域数据的数组
        int n, k;             // 时域和频域的索引
        float x;         // 离散时间信号x
        float Xre, Xim; // DFT变换后的x（实部和虚部）
        float P;         // x的功率图
       
        // 生成在(-1,+1)区间内的随机离散时间信号x
        srand(time(0));
        for (n=0 ; n<N ; ++n) x = ((2.0 * rand()) / RAND_MAX) - 1.0 + sin(PI2 * n * 5.7 / N);
       
        // Calculate DFT and power spectrum up to Nyquist frequency
        // 计算DFT和功率图，上至奈奎斯特频率（在不引入误差的情况下对信号进行采样的最小速率，这是信号中存在的最高频率的两倍）
        int to_sin = 3*N/4; // sin的索引偏移
        int a, b;
        for (k=0 ; k<=N/2 ; ++k)
        {
                Xre = 0; Xim = 0;
                a = 0; b = to_sin;
                for (n=0 ; n<N ; ++n)
                {
                        Xre += x * W;
                        Xim -= x * W;
                        a += k; b += k;
                }
                P = Xre*Xre + Xim*Xim;
        }
       
        // 输出结果到MATLAB / Octave的M文件来绘制图表
        FILE *f = fopen("dftplots.m", "w");
        fprintf(f, "n = ;\n", N-1);
        fprintf(f, "x = [ ");
        for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x);
        fprintf(f, "];\n");
        fprintf(f, "Xre = [ ");
        for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xre);
        fprintf(f, "];\n");
        fprintf(f, "Xim = [ ");
        for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xim);
        fprintf(f, "];\n");
        fprintf(f, "P = [ ");
        for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", P);
        fprintf(f, "];\n");
        fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
        fprintf(f, "xlim()\n", N-1);
        fprintf(f, "subplot(3,1,2)\nplot(,Xre,,Xim)\n", N/2, N/2);
        fprintf(f, "xlim()\n", N/2);
        fprintf(f, "subplot(3,1,3)\nstem(,P)\n", N/2);
        fprintf(f, "xlim()\n", N/2);
        fclose(f);
       
        // 正常退出
        return 0;
}（译者注：有关“旋转因子”请看下图）

当我在Octave中打开生成的M文件时，生成了以下图表：
https://batchloaf.files.wordpress.com/2013/12/dft2plots.png
引入随机信号的主要频率分量，让上图中产生了单个突出峰值，证明了算法的正确性。

顺便说一句，我在上面的例子中使用以下C程序（在PC而不是dsPIC上运行）生成全局旋转因子数组的值：//
// twiddle.c - 打印旋转因子数组
// 作者：Ted Burke
// 最后更新：10-12-2013
//
// 编译：（用mingw或者cygwin，或者直接上Linux。在Linux可以不用加上最后的“.exe”后缀）
//    gcc twiddle.c -o twiddle.exe
//
// 运行：
//    twiddle.exe
// （Linux上则用“./twiddle”来运行）
//

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

#define N 64
#define PI 3.1415926535897932384626433832795
#define PI2 (PI * 2)

int main()
{
    int n;
   
    for (n=0 ; n<N ; ++n)
    {
      printf("%8.5lf", cos(n*PI2/N));
      if (n<N-1) printf(",");
      if ((n+1)%8==0) printf("\n");
    }
    return 0;
}

青山彦樹 发表于 2018-10-24 18:34:16

说到快速傅里叶变换算法，我知道傅里叶这个热力学大佬其实是被热死的。那个年代大概没有中暑的概念吧。

0xAA55 发表于 2018-10-25 09:14:42

根据DFT算法的实现，对于大量数据的DFT变换似乎是可以进行GPU加速的。毕竟它就是针对每个单独的数据做一个求和的运算。

顺带一提谷歌浏览器的拼音插件真是个好东西，可以让我在只有日文输入法的机器上用Chrome输入中文。

查看完整版本: 【翻译】最简单的DFT算法的C代码