技术宅的结界

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

QQ登录

只需一步,快速开始

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

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

[复制链接]

1010

主题

2246

帖子

5万

积分

用户组: 管理员

一只技术宅

UID
1
精华
200
威望
265 点
宅币
16897 个
贡献
33725 次
宅之契约
0 份
在线时间
1608 小时
注册时间
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暴力求解算法。
[C] 纯文本查看 复制代码
//
// 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;
}
用下图的方式来编译并运行上面的代码:

当你运行编译好的程序dft.exe后,它会生成一个名为dftplots.m的M文件,然后可以在MATLAB或Octave中运行来生成结果图。这是dft.exe生成的M文件:
[Plain Text] 纯文本查看 复制代码
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])
我在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代码:
[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[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;
}
(译者注:有关“旋转因子”请看下图)
TwiddleFactor%20Diagram.jpg
当我在Octave中打开生成的M文件时,生成了以下图表:

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

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

0

主题

3

帖子

27

积分

用户组: 初·技术宅

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

1010

主题

2246

帖子

5万

积分

用户组: 管理员

一只技术宅

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

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

本版积分规则

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

GMT+8, 2018-12-15 10:58 , Processed in 0.125917 second(s), 40 queries , Gzip On.

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

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