8 Sep
2004
8 Sep
'04
9:51 p.m.
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