离散傅里叶变换给出不正确的结果
我正在尝试在C中进行离散傅立叶变换。
最初只是蛮力方法。 首先，我让程序打开一个数据文件（幅度），并将数据放入一个数组（仅限一个数据文件，因为我将自己限制为实值输入）。
但是转换看起来不正确，所以我尝试生成一个简单的波函数并检查它是否正确转换。
这是我的代码，没有了钟声和口哨声：
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define M_PI 3.14159265358979323846
//the test wavefunction
double theoretical(double t)
{
double a = sin(M_PI * t) + 2 * sin(2 * M_PI * t) + 4 * sin(4 * M_PI * t);
return a;
}
//-------------------------------------------------------------------------
void dftreal(double inreal[], double outreal[], double outimag[], int linecount)
{
int n, k;
for (k = 0; k < linecount; k++)
{
double sumreal = 0;
double sumimag = 0;
for (n = 0; n < linecount; n++)
{
double angle = 2 * M_PI * n * ( k / (double) linecount);
sumreal += inreal[n] * cos(angle);
sumimag += inreal[n] * sin(angle);
}
outreal[k] = sumreal;
outimag[k] = sumimag;
}
}
//=========================================================================
int main(void)
{
int linecount = 44100;
//creates all necessary arrays
double inreal[linecount], outreal[linecount], outimag[linecount], p[linecount];
FILE *fout = fopen("Output.txt", "w");
for (int i = 0 ; i < linecount ; ++i)
{
inreal[i] = theoretical( i / (double) linecount);
}
//actually computes the transform
dftreal(inreal, outreal, outimag, linecount);
for (int i = 0 ; i < linecount ; ++i)
{
p[i] = 2*(outreal[i] * outreal[i] + outimag[i] * outimag[i]);
fprintf(fout, "%f %f n", (i / (double) linecount), p[i]);
}
fclose(fout);
printf("nEnd of program");
getchar();
return 0;
}
程序编译，完成，但不是功率（频率）图上的几个尖峰，我得到这个： 
。
单个频率或不同的频率给出完全相同的反向浴盆曲线。
我查了几个有关DFT的资料，但我仍然不知道什么是错误的，这个函数似乎没有任何明显的错误：
dftreal
本身。 我想就可能导致问题的原因寻求帮助。 我在Windows 7上使用MinGW编译器。谢谢！
好消息是你的dftreal实现没有问题。
如果可以这样说的话，问题是您使用的测试波形包含相对于您的采样率linecount数非常低的频率linecount 。 相应地，DFT表明能量集中在前几个分箱中，由于波形频率分量不是采样速率的精确倍数，所以一些频谱泄漏到较高分箱中。
如果通过使频率相对于采样频率增加测试波形频率，例如：
//the test wavefunction
double theoretical(double t)
{
double f = 0.1*44100;
double a = sin(2 * M_PI * f * t) + 2 * sin(4 * M_PI * f * t) + 4 * sin(8 * M_PI * f * t);
return a;
}
你应该得到一个情节，例如：
警告：我对此有点生疏
正如我记得的那样，cos / real部分产生频谱，sin / imag部分产生功率谱。 所以，我认为你想要打印的是频谱（这只是outreal[i] ），而不是你所做的（即增加outreal和outimag是错误的）。 你可以用不同的颜色绘制两者，但我会简单地开始。
另外，我会从简单的数据输入开始。
我改变了theoretical做一个单一的余弦波[以及你的原始]。 这是可预测的，你应该得到一个单一的秒杀，你这样做，所以我认为dftreal是正确的。
我还添加了一些命令行选项，以便您可以尝试。
我还发现，划分i / linecount有时会被截断，所以我创建了一个FRAG宏来说明/纠正这个问题。 如果没有这个， cos(2 * pi * t)输入的尖峰以outreal[0] （DC部分？）而不是outreal[1]
此外，出于测试目的，您可以获得更少的样本（例如1000）。 这也可能有助于分散你的情节，以便更好地看到事物。
无论如何，这是代码[请原谅无偿风格的清理]：
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
//#define M_PI 3.14159265358979323846
#define M_2PI (M_PI * 2.0)
int opt_A; // algorithm to use
int linecount; // number of samples
int opt_a; // use absolute value on output
int opt_s; // use squared value on output
int opt_i; // use integer output index
int opt_j; // output imaginary part
int opt_c; // .csv compatibility
time_t todlast;
// the first was your original and would truncate
#if 0
#define FRAG(_i) ((_i) / linecount)
#else
#define FRAG(_i) ((double) (_i) / linecount)
#endif
void
pgr(int k,int n,int count)
{
time_t todnow = time(NULL);
if ((todnow - todlast) >= 1) {
todlast = todnow;
printf("r%d %d ",count,k);
fflush(stdout);
}
}
// the test wavefunction
double
theoretical(double t)
{
double a;
switch (opt_A) {
case 0:
a = 0.0;
a += sin(M_PI * t);
a += 2.0 * sin(2.0 * M_PI * t);
a += 4.0 * sin(4.0 * M_PI * t);
break;
default:
a = cos(M_2PI * t * opt_A);
break;
}
return a;
}
//-------------------------------------------------------------------------
void
dftreal(double inreal[], double outreal[], double outimag[], int linecount)
{
int n;
int k;
for (k = 0; k < linecount; k++) {
double sumreal = 0;
double sumimag = 0;
double omega_k = (M_2PI * k) / (double) linecount;
for (n = 0; n < linecount; n++) {
// omega k:
#if 0
double angle = M_2PI * n * FRAG(k);
#else
double angle = omega_k * n;
#endif
sumreal += inreal[n] * cos(angle);
sumimag += inreal[n] * sin(angle);
pgr(k,n,linecount);
}
outreal[k] = sumreal;
outimag[k] = sumimag;
}
}
//=========================================================================
int
main(int argc,char **argv)
{
char *cp;
--argc;
++argv;
for (; argc > 0; --argc, ++argv) {
cp = *argv;
if (*cp != '-')
break;
switch (cp[1]) {
case 'a': // absolute value
opt_a = ! opt_a;
break;
case 's': // square the output value
opt_s = ! opt_s;
break;
case 'c': // .csv output
opt_c = ! opt_c;
break;
case 'i': // integer output
opt_i = ! opt_i;
break;
case 'j': // imaginary output
opt_j = ! opt_j;
break;
case 'A':
opt_A = atoi(cp + 2);
break;
case 'N':
linecount = atoi(cp + 2);
break;
}
}
if (linecount <= 0)
linecount = 44100 / 2;
// creates all necessary arrays
double inreal[linecount];
double outreal[linecount];
double outimag[linecount];
#if 0
double p[linecount];
#endif
FILE *fout = fopen("Output.txt", "w");
for (int i = 0; i < linecount; ++i)
inreal[i] = theoretical(FRAG(i));
todlast = time(NULL);
// actually computes the transform
dftreal(inreal, outreal, outimag, linecount);
for (int i = 0; i < linecount; ++i) {
#if 0
p[i] = 2 * (outreal[i] * outreal[i] + outimag[i] * outimag[i]);
fprintf(fout, "%f %fn", (i / (double) linecount), p[i]);
#endif
#if 0
fprintf(fout, "%f %f %fn", (i / (double) linecount),
outreal[i] * outreal[i],
outimag[i] * outimag[i]);
#endif
#if 0
fprintf(fout, "%f %fn",
(i / (double) linecount),outreal[i] * outreal[i]);
#endif
if (opt_a) {
outreal[i] = fabs(outreal[i]);
outimag[i] = fabs(outimag[i]);
}
if (opt_s) {
outreal[i] *= outreal[i];
outimag[i] *= outimag[i];
}
if (! opt_c) {
if (opt_i)
fprintf(fout, "%d ",i);
else
fprintf(fout, "%f ",FRAG(i));
}
fprintf(fout, "%f",outreal[i]);
if (opt_j)
fprintf(fout, "%f",outimag[i]);
fprintf(fout, "n");
}
fclose(fout);
printf("nEnd of programn");
//getchar();
return 0;
}
你提到“......变换看起来不对，......”
您是否将结果与其他程序或软件包进行比较以确认结果确实是错误的？ 我最近在C ++和JavaScript中编写了DFT，并将输出与MATLAB，Maple和MathCAD的结果进行了比较。 有时候，差异只是缩放或格式化的问题。
你想输入的原始振幅数有多大？ 如果你在这里发布样本数据，我愿意通过我自己的程序运行它，看看我的程序是否输出与你的结果相同的结果。
链接地址: http://www.djcxy.com/p/33003.html