///
// author: hjjdebug
// date : 2024年 06月 24日 星期一 15:59:53 CST
// descpripton:
// 用c 代码来研究 dft(discrete fourier transform)
///
文章目录
甲: 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.141593
double 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) {
//幅度我加了/ndft
amp[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*fm
int 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 /fs
for (j = 0; j < nb_samples; j++) //每个frame 固定为1024点
{
sample_data[channels*j] = (int)(sin(omiga*t) * 30000); //正负值为-30000 - +30000
for (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*fm
int 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 /fs
for (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不再是一头雾水,也算可以理解和使用了吧.
548
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
