Discrete Fourier transform giving incorrect results
I am trying to do discrete Fourier transforms in C.
Initially just the brute-force method. First I had the program open a data file (of amplitudes) and put the data into an array (just one, since I'm limiting myself to real-valued input).
But the transform looked wrong, so instead I tried to generate a simple wavefunction and check whether it transforms properly.
Here's my code, stripped of the bells and whistles:
#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;
}
The program compiles, completes, but instead of several sharp peaks on a power(frequency) plot, I get this: 
.
A single frequency or different frequencies give the exact same inverted-bathtub-curve.
I checked several sources about the DFT and I still don't know what is going wrong, there do not seem to be any glaring errors with the function:
dftreal
itself. I'd like to ask for help on what could be causing the issue. I'm using the MinGW compiler on Windows 7. Thanks!
The good news is that there is nothing wrong with your dftreal implementation.
The problem, if one could call it that, is that the test waveform you are using includes frequency components which are of very low frequencies relative to your sampling rate linecount . Correspondingly, the DFT shows that the energy is concentrated in the first few bins, with some spectral leakage into higher bins since the waveform frequency components aren't exact multiples of the sampling rate.
If you increase the test waveform frequencies by making the frequency relative to the sampling frequency with for example:
//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;
}
You should get a plot such as:
Caveat: I'm a bit rusty on this
As I remember it, the cos/real part yields the frequency spectrum and the sin/imag part yields the power spectrum. So, what I think you want to print is the frequency spectrum (which is just outreal[i] ) rather than what you did (ie adding outreal and outimag is wrong). You could plot both in different colors but I'd start simply.
Also, I'd start with simpler data input.
I changed the theoretical to do a single cosine wave [as well as your original]. This is predictable and you should get a single spike, which you do, so I think dftreal is correct.
I also added some command line options so you can experiment.
I also found that dividing i / linecount got truncated sometimes, so I created a FRAG macro to illustrate/correct this. Without this, the spike for an input of cos(2 * pi * t) ended up in outreal[0] (DC part?) instead of outreal[1]
Also, for testing purposes, you can get by with far fewer samples (eg 1000). This may also help spread out your plot horizontally to be better able to see things.
Anyway, here's the code [please pardon the gratuitous style cleanup]:
#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;
}
You mention that "... the transform looked wrong, ..."
Did you compare the results with another program or software package to confirm that the results are, indeed, wrong? I recently wrote a DFT in C++ and JavaScript and compared the outputs to results from MATLAB, Maple, and MathCAD. Sometimes the difference is just a matter of scaling or formatting.
How big was the original array of amplitudes you wanted to input? If you post sample data here, I am willing to run it through my own program and see if my program outputs the same results as yours.
链接地址: http://www.djcxy.com/p/33004.html上一篇: Youtube Data API v3 PlaylistItems更新不适用于“稍后观看”播放列表
下一篇: 离散傅里叶变换给出不正确的结果