DSPRelated.com
Forums

FFT output interpretation (FFTW v. Ooura)

Started by rela...@hotmail.co.uk September 30, 2008
I'm new(ish) to FFT and to broaden my knowledge/experience I've been comparing FFTW and the Ooura FFT library. I'm getting a result with the latter that I don't quite understand, however. I have been analysing 1024 samples from a WAV file (16-bit, sr 44100) using FFTW's real-to-complex transform and the equivalent rdft from the Ooura library. I get the same output apart from the DC bin. FFTW gives me 22203 0, whereas Ooura gives me 22203 328721. I understand that the DC is just the average of the samples (22203 in this case), but shouldn't the phase necessarily be 0, as in the FFTW result? Why does Ooura give me 328721 for the imaginary part of the DC?

Many thanks in advance,
Cortez
On Tue, 30 Sep 2008, r...@hotmail.co.uk wrote:

> FFTW gives me 22203 0, whereas Ooura gives me 22203 328721.

DC has no phase, so the phase of it is meaningless.

--
Grzegorz Kraszewski
http://teleinfo.pb.edu.pl/~krashan
Cortez-

> I'm new(ish) to FFT and to broaden my knowledge/experience I've been
> comparing FFTW and the Ooura FFT library. I'm getting a result with the
> latter that I don't quite understand, however. I have been analysing 1024
> samples from a WAV file (16-bit, sr 44100) using FFTW's real-to-complex
> transform and the equivalent rdft from the Ooura library. I get the same
> output apart from the DC bin. FFTW gives me 22203 0, whereas Ooura gives
> me 22203 328721. I understand that the DC is just the average of the
> samples (22203 in this case), but shouldn't the phase necessarily be 0,
> as in the FFTW result? Why does Ooura give me 328721 for the imaginary
> part of the DC?

If you calculate magnitude for bin 0 you should get identical results for both
algorithms, and obviously that's not going to happen, so something is strange. My
guesses would be that you fed something less than the FFT length to the test
programs, so their default value for zero filling was different, or possibly you are
running under Linux and the Gcc tools + optimization options produced executables
that produce some anomaly.

-Jeff
> If you calculate magnitude for bin 0 you should get identical
results for both
> algorithms, and obviously that's not going to happen, so something
is strange. My
> guesses would be that you fed something less than the FFT length to
the test
> programs, so their default value for zero filling was different, or
possibly you are
> running under Linux and the Gcc tools + optimization options
produced executables
> that produce some anomaly.
>
> -Jeff
>

Thanks Jeff. I should perhaps come clean and admit that I'm actually
calling rdft from within Common Lisp via an FFI (I do the same with
FFTW). I'm on WinXP, but I made both a unix dynamic library (.so) and
a .dll using Dev-C++ (it uses gcc). So possibly there was some anomaly
in compiling, as you suggest. I've tried both the .so and the .dll,
and I get the same result. My data is just the first 1024 samples from
the wav file.

Cortez
Cortez-

> > If you calculate magnitude for bin 0 you should get identical
> results for both
> > algorithms, and obviously that's not going to happen, so something
> is strange. My
> > guesses would be that you fed something less than the FFT length to
> the test
> > programs, so their default value for zero filling was different, or
> possibly you are
> > running under Linux and the Gcc tools + optimization options
> produced executables
> > that produce some anomaly.
> >
> > -Jeff
> > Thanks Jeff. I should perhaps come clean and admit that I'm actually
> calling rdft from within Common Lisp via an FFI (I do the same with
> FFTW). I'm on WinXP, but I made both a unix dynamic library (.so) and
> a .dll using Dev-C++ (it uses gcc). So possibly there was some anomaly
> in compiling, as you suggest. I've tried both the .so and the .dll,
> and I get the same result. My data is just the first 1024 samples from
> the wav file.

What if you use CoolEdit or MATLAB or other 'independent' tool? What do they give
for the dc bin?

Also, how do you know you got the first 1024 samples from the .wav file? What are
you using for a .wav reader? If you send me your .wav file I can run it with
Hypersignal and give you back real/imag and mag plots.

If I had to guess I'd say the Ooura test program somehow doesn't "see" quite the same
input data as FFTW.

-Jeff

PS. Please post to the group and cc me.
> What if you use CoolEdit or MATLAB or other 'independent' tool?
What do they give
> for the dc bin?

I've got Scilab, but I'll have to go through the fft tutorial first...

> Also, how do you know you got the first 1024 samples from the .wav
file? What are
> you using for a .wav reader? If you send me your .wav file I can
run it with
> Hypersignal and give you back real/imag and mag plots.

I'm using my own Lisp code for reading wav files. I read the samples
into a vector, and for this test I'm just selecting the first 1024
values. Just to be on the safe side I then read these into a C array
via my FFI code, and pass this as the data to rdft (same with FFTW).

The file I'm using can be found at
http://www.captain.at/data/theevent-zoom.wav (it's a sonification of a
cosmic ray event). I'll send it to you privately.

> If I had to guess I'd say the Ooura test program somehow doesn't
"see" quite the same
> input data as FFTW.
>
> -Jeff
>
> PS. Please post to the group and cc me.

Apologies

Best,
Cortez
Ah, it seems there was a bug in my wav file reading program. The data
I was processing contained errors. Apologies, I shall look at it again
once I've fixed the (minor) bug.

Cortez
Christopher-

> > What if you use CoolEdit or MATLAB or other 'independent' tool?
> What do they give
> > for the dc bin?
>
> I've got Scilab, but I'll have to go through the fft tutorial first...
>
> > Also, how do you know you got the first 1024 samples from the .wav
> file? What are
> > you using for a .wav reader? If you send me your .wav file I can
> run it with
> > Hypersignal and give you back real/imag and mag plots.
>
> I'm using my own Lisp code for reading wav files. I read the samples
> into a vector, and for this test I'm just selecting the first 1024
> values. Just to be on the safe side I then read these into a C array
> via my FFI code, and pass this as the data to rdft (same with FFTW).
>
> The file I'm using can be found at
> http://www.captain.at/data/theevent-zoom.wav (it's a sonification of a
> cosmic ray event). I'll send it to you privately.

Got your .wav file. Here are some Hypersignal captures:

1) .wav file in time domain display:

http://www.signalogic.com/images/evt_zoom.jpg

2) Real/imag data:

http://www.signalogic.com/images/evt_zoom_RI.jpg

The above display used 1024 pt FFT, no overlap, no window. Hypersignal is showing
values of -433.3 for the dc bin real and zero for dc bin imag. I think you would
multiply by 1024 to get the "raw" dc bin value.

3) 2-D Spectrograph:

http://www.signalogic.com/images/evt_zoom_2Dspec.jpg

The above display used 1024 pt frame, 4096 pt FFT, 75% overlap, Hamming window, and a
log-log format to expand the area around your "event".

What do you get for dc bin real after your .wav reader bug is fixed?

-Jeff
Much obliged, Jeff.

> What do you get for dc bin real after your .wav reader bug is fixed?
>
> -Jeff
>

I'm scaling my input between -1 and 1 now, like Matlab/Scilab, and I'm
getting precisely the same values for real and imaginary parts as
Scilab and FFTW - apart, still, from the imaginary part of DC. So,
scaling between -1 and 1 Scilab and FFTW give me an absolute value of
1.3224182 for the real part of the DC bin, and I get the same with the
Ooura rdft (I'm calling both rdft and FFTW's real2complex from Lisp).
But DC has 0.0317687 for its imaginary part with rdft.

I'm beginning to suspect that the parameters I'm passing just aren't
correct. The documentation for the Ooura library is quite terse and to
be honest I'm not entirely sure I've interpreted it correctly. The
exact parameters are -

void rdft(int n, int isgn, double *a, int *ip, double *w)

n is obviously just the length of the data (a), isgn is the direction
of the transform, and ip is the work area for bit reversal (with
length >= 2+sqrt(n/2)). Is w a window parameter? Right now I'm passing
an array filled with zeros (doubles), of length n/2, assuming that
this is equivalent to an unwindowed transform.

Grateful for any clarification,
Chris
Chris-

> Much obliged, Jeff.
>
> > What do you get for dc bin real after your .wav reader bug is fixed?
> >
> > -Jeff
> > I'm scaling my input between -1 and 1 now, like Matlab/Scilab, and I'm
> getting precisely the same values for real and imaginary parts as
> Scilab and FFTW - apart, still, from the imaginary part of DC. So,
> scaling between -1 and 1 Scilab and FFTW give me an absolute value of
> 1.3224182 for the real part of the DC bin, and I get the same with the
> Ooura rdft (I'm calling both rdft and FFTW's real2complex from Lisp).
> But DC has 0.0317687 for its imaginary part with rdft.
>
> I'm beginning to suspect that the parameters I'm passing just aren't
> correct. The documentation for the Ooura library is quite terse and to
> be honest I'm not entirely sure I've interpreted it correctly. The
> exact parameters are -
>
> void rdft(int n, int isgn, double *a, int *ip, double *w)
>
> n is obviously just the length of the data (a), isgn is the direction
> of the transform, and ip is the work area for bit reversal (with
> length >= 2+sqrt(n/2)). Is w a window parameter? Right now I'm passing
> an array filled with zeros (doubles), of length n/2, assuming that
> this is equivalent to an unwindowed transform.

Did you try passing all 1's for *w? If that's actually window data, then "no window"
would be a rectangular window which would be all 1's. Also I wouldn't assume that w
is symmetric. There could be reasons for not doing that.

> Grateful for any clarification,

Given your Scilab and FFTW results -- well-documented routines -- I still question
your .wav reader. If I just take the first 1024 points of event_zoom.wav,
Hypersignal is showing a value of -0.113 for the real dc bin without the 1024 scaling
and about -116 with it. Are you sure the reader is handling the .wav file header
correctly? Hypersignal is showing a 44 byte header in this case.

I can post the 1024 pt plots if you want.

-Jeff