用c 代码来研究 dft(discrete fourier transform)

///
// author: hjjdebug
// date : 2024年 06月 24日星期一 15:59:53 CST
// descpripton:
// 用c 代码来研究 dft(discrete fourier transform)
///

文章目录

甲: DFT 的定义：
乙: 下面给出用c代码实现的dft 公式，验证了手工计算的正确性.
丙：理解还不够深刻? 再来一段程序，专门创建所需要的信号及采样数据.
丁：利用程序进行试验
- 试验1: 创建一个1k赫兹的频率信号,数据格式S16,幅度{-30000,30000}用16k赫兹的采用频率采样.
- 试验2: 用上边数据(1khz信号，16k采样频率)，取16个点进行dft 运算，查看其频谱输出
- 试验3: 仍然用1khz,16k采样频率,但窗口size取32点，计算其dft, 此时频率分辨率是0.5khz
- 试验4: 做一个1khz, 2khz合成的信号，用16k采样频率采样，用dft 分析.

甲: DFT 的定义：

对于长度为n的实数队列x, 它的DFT 定义如下:
X(k) = xigema(x(n)exp(-j2pikn/N)

其中:
k 就是频率的索引,(0<=k<N),
n 就是时间的索引，(0<=n<N）
j 表示虚数单位.
pi 就是3.1415926
exp 就是以e为底的指数函数
exp(-j)表示顺时针旋转的意思？

xigema 就是把它们的值累加的意思.
X(k) 就是以k为自变量的复数数组

假设有一个4点的实数序列x=[1,2,3,4],我们计算它的DFT
则需要计算每个频率点的DFT.
对于k=0,由公式得到：
X(0)= 1exp(-j2pi00/4)+2exp(-j2pi01/4)+3exp(-j2pi02/4)+4exp(-j2pi03/4)=1+2+3+4=10
对于k=1,由公式得到：
X(1)= 1exp(-j2pi10/4)+2exp(-j2pi11/4)+3exp(-j2pi12/4)+4exp(-j2pi13/4)=1+2exp(-jpi/2)+3exp(-jpi)+4exp(-j3pi/2)=1-2j-3+4j=-2+2j

对于k=2,由公式得到：
X(2)= 1exp(-j2pi20/4)+2exp(-j2pi21/4)+3exp(-j2pi22/4)+4exp(-j2pi23/4)=1+2exp(-jpi)+3exp(-j2pi)+4exp(-j3pi)=1-2+3-4=-2

对于k=3,由公式得到：
X(3)= 1exp(-j2pi30/4)+2exp(-j2pi31/4)+3exp(-j2pi32/4)+4exp(-j2pi33/4)=1+2exp(-j3pi/2)+3exp(-j3pi)+4exp(-j9pi/2)=1+2j-3-4j=-2-2j
于是,其结果为: (说实话，我好几个数都计算错了,总是正负号搞错，是看了计算机结果才更正的!)
X[]={10+0j,-2+2j,-2+0j,-2-2j

乙: 下面给出用c代码实现的dft 公式，验证了手工计算的正确性.

$ cat main.cpp 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>#define PI 3.141593double realComput(double x[], int ndft, int k);
double imageComput(double x[], int ndft, int k);typedef struct {double real;double image;
} Complex;Complex* dft(double x[], int ndft) {Complex* dftRes = (Complex*)malloc(ndft * sizeof(Complex));if (dftRes == NULL) {return NULL;}for (int i = 0; i < ndft; ++i) {dftRes[i].real = realComput(x, ndft, i);dftRes[i].image = imageComput(x, ndft, i);
//		printf("%lf + %lfi\n", dftRes[i].real, dftRes[i].image);}return dftRes;
}double realComput(double x[], int ndft, int k) {double realPart = 0;double omiga=2*PI/ndft*k;printf("k:%d,ndft:%d,omiga:%lf\n",k,ndft, omiga);double theta;double delta;double shadow;for (int i = 0; i < ndft; ++i) {
//		delta = x[i] * cos(2 * PI / ndft * k * i);theta = omiga*i;double theta_m = ((float)((k*i)%ndft))/ndft * (2*PI); //不超过2*PI 的角度shadow = cos(theta);delta = x[i] * shadow;printf("k:%d,i:%d,theta:%lf,theta_m:%lf,x:%lf,shadow:%lf,delta:%lf\n",k,i,theta,theta_m,x[i],shadow,delta);realPart += delta;}return realPart;
}//因为image 是累减,所以是顺时针方向旋转
double imageComput(double x[], int ndft, int k) {double imagePart = 0;double theta;double shadow;for (int i = 0; i < ndft; ++i) {imagePart -= x[i] * sin(2 * PI / ndft * k * i);theta=2*PI/ndft*k*i;shadow=sin(theta);printf("--k:%d,i:%d,theta:%lf,x:%lf,shadow:%lf,delta:%lf\n",k,i,theta,x[i],shadow,x[i]*shadow);}return imagePart;
}double* ampSpectrum(Complex* dftRes, int ndft) {double* amp = (double*)malloc(sizeof(double) * ndft);if (amp == NULL) {return NULL;}for (int i = 0; i < ndft; ++i) {//幅度我加了/ndftamp[i] = sqrt(dftRes[i].real * dftRes[i].real+ dftRes[i].image* dftRes[i].image)/ndft;}return amp;
}double* phaseSpectrum(Complex* dftRes, int ndft) {double* phase = (double*)malloc(sizeof(double) * ndft);if (phase == NULL) {return NULL;}for (int i = 0; i < ndft; ++i) {phase[i] = atan2(dftRes[i].image, dftRes[i].real);}return phase;
}#define  WIN_S 16
void testDFT() {
//	double x[] = { 1, 2, 3, 4};double x[WIN_S];short sht[WIN_S];int res;FILE *fp=fopen("audio.data","rb");if(!fp){printf("error open file\n");exit(1);}//读取32个符号整数，进行dft变换res=fread(sht,sizeof(short),WIN_S,fp);if(res!=WIN_S){printf("error read file\n");exit(1);}int ndft = sizeof(x) / sizeof(double);for(int i=0;i<ndft;i++){x[i]=sht[i];}Complex* dftRes = dft(x, ndft);for (int i = 0; i < ndft; ++i) {printf("%lf + %lfi\n", dftRes[i].real, dftRes[i].image);}printf("\n");double* ampRes = ampSpectrum(dftRes, ndft);for (int i = 0; i < ndft; ++i) {printf("%lf, ", ampRes[i]);}printf("\n");double* phaRes = phaseSpectrum(dftRes, ndft);for (int i = 0; i < ndft; ++i) {
//		printf("%lf, ", phaRes[i]);}printf("\n");free(dftRes);free(ampRes);free(phaRes);
}int main() {testDFT();return 0;
}

输出复数值：
10.000000 + 0.000000i
-1.999998 + 2.000001i
-2.000000 + 0.000003i
-2.000005 + -1.999997i

由此你还会想到更多的问题：例如:

输出频率是什么? 答: 输出的频率是Fs/N, Fs/N2, Fs/N3… Fs/N*(N-1), 其中N 是采样点数,就是说输出频率是采样频率N等分.
推论：要想使输出频率
采样频率是什么? 答: 采样频率Fs 就是数据采样采到1,2,3,4 数据时使用的频率. 要求采样频率至少比信号频率大2倍
信号频率是什么? 答: 就是输入的正弦（或余弦）波频率,可想象成匀速圆周运动在y轴上的投影.
dft输出的幅度是什么意思?

丙：理解还不够深刻? 再来一段程序，专门创建所需要的信号及采样数据.

create_audio.cpp

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(int argc, char **argv)
{if (argc != 4) {fprintf(stderr, "Usage: %s <output_file> <fm> <fs> \n""fm is signal frequency, fs is sample frequency\n""Example: %s audio.data 1000 16000 \n", argv[0],argv[0]);exit(1);}char *filename = argv[1];FILE *file = fopen(filename, "wb");if (!file) {fprintf(stderr, "Could not open destination file %s\n", filename);exit(1);}int fm=0,fs=0;sscanf(argv[2],"%d",&fm);if(fm==0) exit(2);sscanf(argv[3],"%d",&fs);if(fs==0) exit(3);/* encode a single tone sound */int i,j,k;int nb_samples = 1024;
//	int channels = 2;int channels = 1;int sample_rate=fs;float t = 0;//这里由2个频率, 信号频率,每秒钟转过一个omiga角,每次采样时在y轴投影都不同//采样频率,决定了下一次的采样时间
//对sin函数,2pi 弧度是一个周期, 均匀圆周运动,其角频率为2*pi*fmint		omiga = 2 * M_PI * fm;//每次采样的时间递长值float tincr = 1.0 / sample_rate; typedef unsigned short uint16_t;//分配内存，存储一个frame 的数据uint16_t *sample_data = (uint16_t*)malloc(nb_samples * channels *sizeof(uint16_t)); //内存容量要考虑上通道个数for (i = 0; i < 200; i++) { // create 200 个 frame, 可持续时间200 * 1024 /fsfor (j = 0; j < nb_samples; j++) //每个frame 固定为1024点{sample_data[channels*j] = (int)(sin(omiga*t) * 30000); //正负值为-30000 - +30000for (k = 1; k < channels; k++)sample_data[channels*j + k] = sample_data[channels*j];t += tincr;}fwrite(sample_data, 1, nb_samples * channels *sizeof(uint16_t), file);}fprintf(stderr, "play the output file with the command:\n""ffplay -f s16le -ac %d -ar %d %s\n",channels,sample_rate,filename);		//指明数据采样格式，通道数，采样率就可以播放原始数据了fclose(file);return 0;
}

丁：利用程序进行试验

试验1: 创建一个1k赫兹的频率信号,数据格式S16,幅度{-30000,30000}用16k赫兹的采用频率采样.

采样过程可想象成如下图示:

采样的数据如下表所示:
0000h: 00 00 D8 2C DC 52 43 6C 30 75 45 6C DE 52 DA 2C …Ø,ÜRCl0uElÞRÚ,
0010h: 02 00 2B D3 26 AD BE 93 D0 8A BA 93 20 AD 23 D3 …+Ó&¾“ÐŠº“ #Ó
0020h: FB FF D3 2C D8 52 41 6C 30 75 47 6C E2 52 DF 2C ûÿÓ,ØRAl0uGlâRß,
0030h: 08 00 30 D3 2A AD C0 93 D1 8A B8 93 1C AD 1E D3 …0Ó*À“ÑŠ¸“..Ó
0040h: F5 FF CD 2C D4 52 3F 6C 2F 75 49 6C E6 52 E5 2C õÿÍ,ÔR?l/uIlæRå,
0050h: 0E 00 35 D3 2E AD C2 93 D1 8A B6 93 18 AD 19 D3 …5Ó.Â“ÑŠ¶“..Ó

试验2: 用上边数据(1khz信号，16k采样频率)，取16个点进行dft 运算，查看其频谱输出

 由于采样率是16K， 16个点其频率分辨率是1K， 能够把我们的1k信号给过滤出来.其频谱图(振幅图)为:

2.000000, 239998.747258, 9.897716, 5.690789, 4.487543, 3.588161, 2.606504, 3.878767, 2.000068, 3.879591, 2.581657, 3.547289, 4.425996, 5.653179, 9.610440, 239998.947999,

我们一下就看到了第2个值239998!! 显然它对应着1K频率，它这个幅度值代表了什么?

试验3: 仍然用1khz,16k采样频率,但窗口size取32点，计算其dft, 此时频率分辨率是0.5khz

结果如下:
6.000000, 20.862890, 479996.500013, 33.166659, 19.795433, 13.383284, 11.795656, 10.343337, 10.027572, 6.892723, 6.725493, 4.885620, 5.213012, 4.484607, 7.032958, 6.083789, 6.000091, 6.060954, 7.024984, 4.464189, 5.163309, 4.845143, 6.652847, 6.829602, 9.917541, 10.198174, 11.703062, 13.056493, 19.220876, 31.801808, 479996.901704, 18.721646,
分析: 耀眼的是第3项479996,(我们只分析前半部分频率<Fs/2, 因为后半部分与前半部分是共轭的.)
目测比16采样点又大了一倍.
479996/32=14999.9
239998/16=14999.9

研究代码,这个幅度值由实部，虚部用勾股计算而来.实部与虚部由时序数值累加而来，其数值只所以很大，
是对时序采样点累加的结果,可以对其归一化(除以采样点数),就消除了采样点数的影响.
归一化后的幅度谱图（16点):
0.125000, 14999.921704, 0.618607, 0.355674, 0.280471, 0.224260, 0.162907, 0.242423, 0.125004, 0.242474, 0.161354, 0.221706, 0.276625, 0.353324, 0.600652, 14999.934250,

那这个14999.9是什么东西? 这个描述就有点复杂了,采样点的最大幅值是30000.
16个采样点数据投射到1khz 的旋转转盘上,实部与虚部共同作用形成的合力是15000.
而投射到其它频率上，形成的合力几乎可以忽略(目测小于1).

试验4: 做一个1khz, 2khz合成的信号，用16k采样频率采样，用dft 分析.

这就需要修改代码了,不过很简单的. (注意修改幅值大小使其和不要超限制65536/2）
代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(int argc, char **argv)
{if (argc != 4) {fprintf(stderr,"Description: add a 5k frequecy to output.\n""Usage: %s <output_file> <fm> <fs> \n""fm is signal frequency, fs is sample frequency\n""Example: %s audio.data 1000 16000 \n", argv[0],argv[0]);exit(1);}char *filename = argv[1];FILE *file = fopen(filename, "wb");if (!file) {fprintf(stderr, "Could not open destination file %s\n", filename);exit(1);}int fm=0,fs=0;sscanf(argv[2],"%d",&fm);if(fm==0) exit(2);sscanf(argv[3],"%d",&fs);if(fs==0) exit(3);/* encode a single tone sound */int i,j;int nb_samples = 1024;int channels = 1;int sample_rate=fs;float t = 0;//这里由2个频率, 信号频率,每秒钟转过一个omiga角,每次采样时在y轴投影都不同//采样频率,决定了下一次的采样时间//对sin函数,2pi 弧度是一个周期, 均匀圆周运动,其角频率为2*pi*fmint		omiga = 2 * M_PI * fm;int omiga2 = 2 *M_PI * 5000; //5k 频率信号int val_5k;//每次采样的时间递长值float tincr = 1.0 / sample_rate; typedef unsigned short uint16_t;//分配内存，存储一个frame 的数据uint16_t *sample_data = (uint16_t*)malloc(nb_samples * sizeof(uint16_t)); for (i = 0; i < 200; i++) { // create 200 个 frame, 可持续时间200 * 1024 /fsfor (j = 0; j < nb_samples; j++) //每个frame{val_5k = (int)(sin(omiga2*t) * 10000);sample_data[j] = (int)(sin(omiga*t) * 10000)+val_5k; //正负值为-10000 - +10000, 2个值想加不会超过2万，在数据表达范围以内t += tincr;}fwrite(sample_data, 1, nb_samples * sizeof(uint16_t), file);}fprintf(stderr, "play the output file with the command:\n""ffplay -f s16le -ac %d -ar %d %s\n",channels,sample_rate,filename);		//指明数据采样格式，通道数，采样率就可以播放原始数据了fclose(file);return 0;
}

直接看dft输出吧.

0.187500, 4999.898941, 0.276690, 0.635727, 0.764330, 4999.632634, 0.739231, 0.515267, 0.187523, 0.512679, 0.728653, 4999.630904, 0.758657, 0.634583, 0.276951, 4999.905571,
看到1k,5khz处确实有频率输出了. 大于8k就不用看了,它们与前8k共轭.

于是这3段程序演示了如何对信号进行时域采样, 如何进行dft转换并验证输入频率等. 如此dtf不再是一头雾水,也算可以理解和使用了吧.