Hi. I am trying to finally get a grip on what a FFT really does by the learning by doing approach. I ripped some code from qjacktuner (thanks to Karsten Wiese) and rewrote it such that it uses fftw3 (which I'd prefer as a backend for now). However, the code executes, but returns very strange results, and so I am asking the DSP coders of you if you perhaps can see the probably quite obvious error I am making? When I pluck the A string on my guitar, I receive output like this: Sample rate: 48000 0 = {freq=143.807441, amp=8.776460} 1 = {freq=162.811147, amp=9.086173} 2 = {freq=181.814852, amp=9.604586} 3 = {freq=200.818557, amp=9.631306} 0 = {freq=148.405443, amp=9.330607} 1 = {freq=271.871957, amp=nan} 2 = {freq=148.405443, amp=49.469994} 3 = {freq=148.405443, amp=45.180032} 0 = {freq=19.003705, amp=47.543385} 2 = {freq=19.003705, amp=66.636058} 3 = {freq=394.003705, amp=63.710202} 0 = {freq=19.003705, amp=19.350767} 2 = {freq=19.003705, amp=28.120214} 3 = {freq=394.003705, amp=56.406310} 0 = {freq=19.003705, amp=42.607500} 2 = {freq=19.003705, amp=78.669607} 3 = {freq=241.974062, amp=38.201785} 0 = {freq=19.003705, amp=2.835922} 2 = {freq=241.974062, amp=93.692646} 3 = {freq=394.003705, amp=48.011389} 0 = {freq=19.003705, amp=45.488988} 2 = {freq=19.003705, amp=46.009611} 3 = {freq=394.003705, amp=60.820928} 2 = {freq=19.003705, amp=54.438861} 3 = {freq=394.003705, amp=47.402751} In fact, it seems no matter what sound input I give it, it ends up reporting 19.003705 as the loudest frequency. This is obviously wrong. Insight appreciated :-) /* foolowing code is only partial, so it doesn't compile, I can provide * you with the full source of the program if you want. * Variable "rate" is the current sampling rate (48000 in my case). */ #include <fftw3.h> #define M_PI 3.14159265358979323846 #define MAX_FFT_LENGTH 16384 static float sampleFIFO[MAX_FFT_LENGTH]; static float *sample = NULL; static float lastPhase[MAX_FFT_LENGTH/2+1]; static int fftSize; static double *fftIn; static fftw_complex *fftOut; static fftw_plan fftPlan; void fftInit (int size) { fftIn = fftw_malloc(sizeof(double) * size); fftOut = fftw_malloc(sizeof(fftw_complex) * (size/2+1)); fftSize = size; fftPlan = fftw_plan_dft_r2c_1d(fftSize, fftIn, fftOut, FFTW_MEASURE); memset(sampleFIFO, 0, MAX_FFT_LENGTH*sizeof(float)); memset(lastPhase, 0, (MAX_FFT_LENGTH*sizeof(float))/2); } void fftMeasure (int nframes, int osamp, float *indata) { double freqPerBin, expct; long i,k, qpd, inFifoLatency, stepSize; stepSize = fftSize/osamp; freqPerBin = rate/(double)fftSize; expct = 2.*M_PI*(double)stepSize/(double)fftSize; inFifoLatency = fftSize-stepSize; if (!sample) sample = sampleFIFO + inFifoLatency; for (i=0; i<nframes; i++) { *sample++ = indata[i]; if (sample-sampleFIFO >= fftSize) { sample = sampleFIFO + inFifoLatency; for (k=0; k<fftSize; k++) { double window = -.5*cos(2.*M_PI*(double)k/(double)fftSize)+.5; fftIn[k] = sampleFIFO[k] * window; } fftw_execute(fftPlan); for (k=0; k<=fftSize/2; k++) { double magn, phase, tmp, real, imag; real = fftOut[k][0]; imag = fftOut[k][1]; /* compute magnitude and phase */ magn = 2.*sqrt(real*real + imag*imag); phase = atan2(imag,real); /* compute phase difference */ tmp = phase - lastPhase[k]; lastPhase[k] = phase; /* subtract expected phase difference */ tmp -= (double)k*expct; /* map delta phase into +/- Pi interval */ qpd = (long)(tmp / M_PI); if (qpd >= 0) qpd += qpd&1; else qpd -= qpd&1; tmp -= M_PI*(double)qpd; /* get deviation from bin frequency from the +/- Pi interval */ tmp = osamp*tmp/(2.*M_PI); /* compute the k-th partials' true frequency */ tmp = (double)k*freqPerBin + tmp*freqPerBin; if (tmp > 0.0) { printf("%d = {freq=%f, amp=%f}\n", k, tmp, magn); } } /* move input FIFO */ for (k=0; k<inFifoLatency; k++) sampleFIFO[k] = sampleFIFO[k+stepSize]; } } } void fftFree () { fftw_destroy_plan(fftPlan); fftw_free(fftIn); fftw_free(fftOut); } -- Thanks, Mario