技术宅的结界

 找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 1114|回复: 2
收起左侧

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

[复制链接]

1045

主题

2338

帖子

6万

积分

用户组: 管理员

一只技术宅

UID
1
精华
218
威望
294 点
宅币
18348 个
贡献
37770 次
宅之契约
0 份
在线时间
1753 小时
注册时间
2014-1-26
发表于 2018-10-24 10:10:30 | 显示全部楼层 |阅读模式

欢迎访问技术宅的结界,请注册或者登录吧。

您需要 登录 才可以下载或查看,没有帐号?立即注册→加入我们

x
翻译原文: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
1.png

De Moivre的定理表明:
2.png

因此,我们可以将DFT求和重写为X[k]的实部和虚部的表达式,如下所示(假设x[n]为实数)。
3.png

其中:
4.png


5.png

作为我学习过程的一部分,我一直在尝试编写C和MATLAB / Octave中DFT和FFT的不同实现。这让我感觉很爽,因为我写了个非常简单的DFT暴力求解算法。[mw_shl_code=c,true]//
// 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[N];           // 离散时间信号x
        float Xre[N], Xim[N]; // DFT变换后的x(实部和虚部)
        float P[N];           // x的功率图
         
        // 生成在(-1,+1)区间内的随机离散时间信号x
        srand(time(0));
        for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0;
         
        // 用暴力算法计算x的DFT变换
        for (k=0 ; k<N ; ++k)
        {
                // X[k]的实部
                Xre[k] = 0;
                for (n=0 ; n<N ; ++n) Xre[k] += x[n] * cos(n * k * PI2 / N);
                 
                // X[k]的虚部
                Xim[k] = 0;
                for (n=0 ; n<N ; ++n) Xim[k] -= x[n] * sin(n * k * PI2 / N);
                 
                // 第k个频率点的功率
                P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
        }
         
        // 输出结果到MATLAB / Octave的M文件来绘制图表
        FILE *f = fopen("dftplots.m", "w");
        fprintf(f, "n = [0:%d];\n", N-1);
        fprintf(f, "x = [ ");
        for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
        fprintf(f, "];\n");
        fprintf(f, "Xre = [ ");
        for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xre[k]);
        fprintf(f, "];\n");
        fprintf(f, "Xim = [ ");
        for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xim[k]);
        fprintf(f, "];\n");
        fprintf(f, "P = [ ");
        for (k=0 ; k<N ; ++k) fprintf(f, "%f ", P[k]);
        fprintf(f, "];\n");
        fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
        fprintf(f, "xlim([0 %d])\n", N-1);
        fprintf(f, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
        fprintf(f, "xlim([0 %d])\n", N-1);
        fprintf(f, "subplot(3,1,3)\nstem(n,P)\n");
        fprintf(f, "xlim([0 %d])\n", N-1);
        fclose(f);
         
        // 正常退出
        return 0;
}[/mw_shl_code]用下图的方式来编译并运行上面的代码:

当你运行编译好的程序dft.exe后,它会生成一个名为dftplots.m的M文件,然后可以在MATLAB或Octave中运行来生成结果图。这是dft.exe生成的M文件:[mw_shl_code=text,true]n = [0:15];
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([0 15])
subplot(3,1,2)
plot(n,Xre,n,Xim)
xlim([0 15])
subplot(3,1,3)
stem(n,P)
xlim([0 15])[/mw_shl_code]我在Octave中打开了dftplots.m,如下所示:

它产生了如下图表:

最上面的是原始随机离散时间序列x[n]。第二个图是X[k]的实部和虚部。注意X[k]的对称性,正如我们对实数x[n]所预想的那样。最后那个图是P[k]也就是x[n]的功率图。P[k]的每个值是对应的复数值X[k]的大小的平方。

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

这是修改过的C代码:[mw_shl_code=c,true]//
// 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[N];           // 离散时间信号x
        float Xre[N], Xim[N]; // DFT变换后的x(实部和虚部)
        float P[N];           // x的功率图
         
        // 生成在(-1,+1)区间内的随机离散时间信号x
        srand(time(0));
        for (n=0 ; n<N ; ++n) x[n] = ((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[k] = 0; Xim[k] = 0;
                a = 0; b = to_sin;
                for (n=0 ; n<N ; ++n)
                {
                        Xre[k] += x[n] * W[a%N];
                        Xim[k] -= x[n] * W[b%N];
                        a += k; b += k;
                }
                P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
        }
         
        // 输出结果到MATLAB / Octave的M文件来绘制图表
        FILE *f = fopen("dftplots.m", "w");
        fprintf(f, "n = [0:%d];\n", N-1);
        fprintf(f, "x = [ ");
        for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
        fprintf(f, "];\n");
        fprintf(f, "Xre = [ ");
        for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xre[k]);
        fprintf(f, "];\n");
        fprintf(f, "Xim = [ ");
        for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xim[k]);
        fprintf(f, "];\n");
        fprintf(f, "P = [ ");
        for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", P[k]);
        fprintf(f, "];\n");
        fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
        fprintf(f, "xlim([0 %d])\n", N-1);
        fprintf(f, "subplot(3,1,2)\nplot([0:%d],Xre,[0:%d],Xim)\n", N/2, N/2);
        fprintf(f, "xlim([0 %d])\n", N/2);  
        fprintf(f, "subplot(3,1,3)\nstem([0:%d],P)\n", N/2);
        fprintf(f, "xlim([0 %d])\n", N/2);
        fclose(f);
         
        // 正常退出
        return 0;
}[/mw_shl_code](译者注:有关“旋转因子”请看下图)
TwiddleFactor%20Diagram.jpg
当我在Octave中打开生成的M文件时,生成了以下图表:

引入随机信号的主要频率分量,让上图中产生了单个突出峰值,证明了算法的正确性。

顺便说一句,我在上面的例子中使用以下C程序(在PC而不是dsPIC上运行)生成全局旋转因子数组的值:[mw_shl_code=c,true]//
// 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;
}[/mw_shl_code]

本帖被以下淘专辑推荐:

0

主题

3

帖子

27

积分

用户组: 初·技术宅

UID
4330
精华
0
威望
0 点
宅币
24 个
贡献
0 次
宅之契约
0 份
在线时间
2 小时
注册时间
2018-9-29
发表于 2018-10-24 18:34:16 | 显示全部楼层
说到快速傅里叶变换算法,我知道傅里叶这个热力学大佬其实是被热死的。那个年代大概没有中暑的概念吧。

1045

主题

2338

帖子

6万

积分

用户组: 管理员

一只技术宅

UID
1
精华
218
威望
294 点
宅币
18348 个
贡献
37770 次
宅之契约
0 份
在线时间
1753 小时
注册时间
2014-1-26
 楼主| 发表于 2018-10-25 09:14:42 | 显示全部楼层
根据DFT算法的实现,对于大量数据的DFT变换似乎是可以进行GPU加速的。毕竟它就是针对每个单独的数据做一个求和的运算。

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

本版积分规则

QQ|申请友链||Archiver|手机版|小黑屋|技术宅的结界 ( 滇ICP备16008837号 )

GMT+8, 2019-12-7 05:33 , Processed in 0.083868 second(s), 34 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2020, Tencent Cloud.

快速回复 返回顶部 返回列表