I'm baffled by the results I'm getting from FFT and would appreciate any help.
I'm using FFTW 3.2.2 but have gotten similar results with other FFT implementations (in Java). When I take the FFT of a sine wave, the scaling of the result depends on the frequency (Hz) of the wave--specifically, whether it's close to a whole number or not. The resulting values are scaled really small when the frequency is near a whole number, and they're orders of magnitude larger when the frequency is in between whole numbers. This graph shows the magnitude of the spike in the FFT result corresponding to the wave's frequency, for different frequencies. Is this right??
I checked that the inverse FFT of the FFT is equal to the original sine wave times the number of samples, and it is. The shape of the FFT also seems to be correct.
It wouldn't be so bad if I were analyzing individual sine waves, because I could just look for the spike in the FFT regardless of its height. The problem is that I want to analyze sums of sine waves. If I'm analyzing a sum of sine waves at, say, 440 Hz and 523.25 Hz, then only the spike for the one at 523.25 Hz shows up. The spike for the other is so tiny that it just looks like noise. There must be some way to make this work because in Matlab it does work-- I get similar-sized spikes at both frequencies. How can I change the code below to equalize the scaling for different frequencies?
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <fftw3.h>
#include <cstdio>
using namespace std;
const double PI = 3.141592;
/* Samples from 1-second sine wave with given frequency (Hz) */
void sineWave(double a[], double frequency, int samplesPerSecond, double ampFactor);
int main(int argc, char** argv) {
/* Args: frequency (Hz), samplesPerSecond, ampFactor */
if (argc != 4) return -1;
double frequency = atof(argv[1]);
int samplesPerSecond = atoi(argv[2]);
double ampFactor = atof(argv[3]);
/* Init FFT input and output arrays. */
double * wave = new double[samplesPerSecond];
sineWave(wave, frequency, samplesPerSecond, ampFactor);
double * fftHalfComplex = new double[samplesPerSecond];
int fftLen = samplesPerSecond/2 + 1;
double * fft = new double[fftLen];
double * ifft = new double[samplesPerSecond];
/* Do the FFT. */
fftw_plan plan = fftw_plan_r2r_1d(samplesPerSecond, wave, fftHalfComplex, FFTW_R2HC, FFTW_ESTIMATE);
fftw_execute(plan);
memcpy(fft, fftHalfComplex, sizeof(double) * fftLen);
fftw_destroy_plan(plan);
/* Do the IFFT. */
fftw_plan iplan = fftw_plan_r2r_1d(samplesPerSecond, fftHalfComplex, ifft, FFTW_HC2R, FFTW_ESTIMATE);
fftw_execute(iplan);
fftw_destroy_plan(iplan);
printf("%s,%s,%s", argv[1], argv[2], argv[3]);
for (int i = 0; i < samplesPerSecond; i++) {
printf("\t%.6f", wave[i]);
}
printf("\n");
printf("%s,%s,%s", argv[1], argv[2], argv[3]);
for (int i = 0; i < fftLen; i++) {
printf("\t%.9f", fft[i]);
}
printf("\n");
printf("\n");
printf("%s,%s,%s", argv[1], argv[2], argv[3]);
for (int i = 0; i < samplesPerSecond; i++) {
printf("\t%.6f (%.6f)", ifft[i], samplesPerSecond * wave[i]); // actual and expected result
}
delete[] wave;
delete[] fftHalfComplex;
delete[] fft;
delete[] ifft;
}
void sineWave(double a[], double frequency, int samplesPerSecond, double ampFactor) {
for (int i = 0; i < samplesPerSecond; i++) {
double time = i / (double) samplesPerSecond;
a[i] = ampFactor * sin(2 * PI * frequency * time);
}
}