我使用这个DFT实现:
/*
Direct fourier transform
*/
int DFT(int dir,int m,double *x1,double *y1)
{
long i,k;
double arg;
double cosarg,sinarg;
double *x2=NULL,*y2=NULL;
x2 = malloc(m*sizeof(double));
y2 = malloc(m*sizeof(double));
if (x2 == NULL || y2 == NULL)
return(FALSE);
for (i=0;i<m;i++) {
x2[i] = 0;
y2[i] = 0;
arg = - dir * 2.0 * 3.141592654 * (double)i / (double)m;
for (k=0;k<m;k++) {
cosarg = cos(k * arg);
sinarg = sin(k * arg);
x2[i] += (x1[k] * cosarg - y1[k] * sinarg);
y2[i] += (x1[k] * sinarg + y1[k] * cosarg);
}
}
/* Copy the data back */
if (dir == 1) {
for (i=0;i<m;i++) {
x1[i] = x2[i] / (double)m;
y1[i] = y2[i] / (double)m;
}
} else {
for (i=0;i<m;i++) {
x1[i] = x2[i];
y1[i] = y2[i];
}
}
free(x2);
free(y2);
return(TRUE);
}放置在这里的http://paulbourke.net/miscellaneous/dft/
第一个问题是为什么在应用直接变换(dir=1)之后我们应该标度值?我读了一些关于DFT实现的想法,但是没有发现任何关于它的东西。
作为输入,我使用cos和1024采样频率。
#define SAMPLES 2048
#define ZEROES_NUMBER 512
double step = PI_2/(SAMPLES-2*ZEROES_NUMBER);
for(int i=0; i<SAMPLES; i++)
{
/*
* Fill in the beginning and end with zeroes
*/
if(i<ZEROES_NUMBER || i > SAMPLES-ZEROES_NUMBER)
{
samplesReal[i] = 0;
samplesImag[i] = 0;
}
/*
* Generate one period cos with 1024 samples
*/
else
{
samplesReal[i] = cos(step*(double)(i-ZEROES_NUMBER));
samplesImag[i] = 0;
}
}对于绘图,我删除了上面提到的缩放,因为输出值变得非常小,并且不可能绘制图表。
我得到了这样的振幅和相位图:


正如你所看到的,相位总是0,振幅谱是相反的。为什么?
下面是我的更具可读性的版本,不需要缩放就会产生相同的结果:
void DFT_transform(double complex* samples, int num, double complex* res)
{
for(int k=0; k<num; k++)
{
res[k] = 0;
for(int n=0; n<num; n++)
{
double complex Wkn = cos(PI_2*(double)k*(double)n/(double)num) -
I*sin(PI_2*(double)k*(double)n/(double)num);
res[k] += samples[n]*Wkn;
}
}
}发布于 2015-10-08 20:02:51
好了伙计们。我很高兴地说,这个实现是有效的。问题是绘制图形的方法错误,缺乏理解公式。
它是如何工作的

如您所见,有k变量,用于改变频率。所以频率是ν=k/ T,其中T是获取样本所需的一段时间。T = N/S,其中S是你的采样频率。然后你可以发现你的频率是v = S*k/N
因此,当你得到你的结果时,你应该计算每个点的频率,去掉所有在S/2上的东西,然后才画出图形震级=震级(频率)。这是我以前不明白的。希望它能对某人有所帮助。
我有一些图表。
Sin 100赫兹

Sin 100 Cos +Cos 200 Cos

Sin 100/2+ (Cos 200/2)/2

正如你所看到的,频率和相关的震级显示出来了。有一个问题,缩放,但这并不重要,如果我们想要确定的频率,在一个信号。
感谢@PaulR
https://stackoverflow.com/questions/33013558
复制相似问题