Forums

FFTW Accuracy vs Matlab Accuracy

Started by Clyde April 12, 2007
Hi,
i'm using FFTW to computer determine the SNR of a recorded sinusoidal 
signal. (1khz sinus for 5seconds @ 48khz)

I'm experiencing a strange behaviour:
if i compute the fft with FFTW i obtain really different values from 
what matlab compute.
To be more specific:

i compute a real 2 real fft with:
ApplyWindow(BLACKMANHARRIS);// <- some windowing to my signal
plan = fftw_plan_r2r_1d(dim, in[0], out[0], FFTW_R2HC, FFTW_ESTIMATE);
fftw_execute(plan);

If i try to do the same with matlab or any other software (like wavelab 
or audacity) i can see that with FFTW i obtain higher magnitudinae at 
low frequencies ( < 200hz ~).
I tried 6 different windowing  methods, but it didn't worked.
I tried to compute a
fftw_plan_dft_1d(dim, in, out, FFTW_FORWARD, FFTW_ESTIMATE) but it 
didn't worked as well.
Is there something i'm missing?

--

Clyde
Clyde wrote:
> Hi, > i'm using FFTW to computer determine the SNR of a recorded sinusoidal > signal. (1khz sinus for 5seconds @ 48khz) > > I'm experiencing a strange behaviour: > if i compute the fft with FFTW i obtain really different values from > what matlab compute.
Matlab also uses FFTW, so you must be doing something wrong.
"Andor" <andor.bariska@gmail.com> wrote in news:1176377377.655832.61120
@q75g2000hsh.googlegroups.com:

> Clyde wrote: >> Hi, >> i'm using FFTW to computer determine the SNR of a recorded sinusoidal >> signal. (1khz sinus for 5seconds @ 48khz) >> >> I'm experiencing a strange behaviour: >> if i compute the fft with FFTW i obtain really different values from >> what matlab compute. > > Matlab also uses FFTW, so you must be doing something wrong. >
Ditto. I haven't used FFTW directly, so I'm grasping, but I do recall some people posting about having problems getting the argument order right, or possibly the shape of the inputs, in FFTW. Double and triple check your function call. Also, Matlab scales the FFT a little funny, so you might be off by a constant 1/N or something, but both answers should have the same shape. -- Scott Reverse name to reply
Clyde wrote:

> Hi, > i'm using FFTW to computer determine the SNR of a recorded sinusoidal > signal. (1khz sinus for 5seconds @ 48khz)
> If i try to do the same with matlab or any other software (like wavelab > or audacity) i can see that with FFTW i obtain higher magnitudinae at > low frequencies ( < 200hz ~).
Can you show us the shape of the outpute data? Maybe it'll help.
> I tried 6 different windowing methods, but it didn't worked.
How it look without windowing?
> I tried to compute a > fftw_plan_dft_1d(dim, in, out, FFTW_FORWARD, FFTW_ESTIMATE) but it > didn't worked as well. > Is there something i'm missing?
You can try to generate one harmonic signal and do FFTW (without windowing) on it and then show us the result. -- pgw
pgw wrote:

> How it look without windowing?
They are different, and the problem still remains.
> You can try to generate one harmonic signal and do FFTW (without windowing) > on it and then show us the result.
Here is the result: http://i45.photobucket.com/albums/f77/yesitookmypills/FFT.png (blue is FFTW, red is Matlab FFT) I used fftw_plan_r2r_1d(dim, in[0], out[0], FFTW_R2HC, FFTW_ESTIMATE) and i compute the magnitude in this way: magnitude[i] = sqrt(out[i]*out[i]); Here's the data with a small script in matlab: http://rapidshare.com/files/25643904/FFTWvsFFTMatlab.zip.html -- Clyde
Scott Seidman wrote:

> Ditto. I haven't used FFTW directly, so I'm grasping, but I do recall some > people posting about having problems getting the argument order right, or > possibly the shape of the inputs, in FFTW. Double and triple check your > function call. Also, Matlab scales the FFT a little funny, so you might be > off by a constant 1/N or something, but both answers should have the same > shape.
I'm not a FFT guru, and FFTW has so many option that it can be trivial.. I just used my complex analysis book, ant a little of intuition. With FFTW you can compute 11 kind of real to real fft, but i have big troubles understanding the differences between them. The scale should be the same. The same "feature" is "documented" in fftw (at least they warn you). -- Clyde
Andor wrote:
> Matlab also uses FFTW, so you must be doing something wrong.
Yeah, i agree. But what i'm trying to underlying is not "who's the best". I just want to understand what differs, and where i do something wrong. -- Clyde
Clyde <iloveuppercaseEREKOSE@EMAIL.IT> wrote in news:461e738c$0$4786
$4fafbaef@reader4.news.tin.it:

> fftw_plan_r2r_1d
Don't you need to deal with the "half complex" format that fftw is asked to produce before you deal with the magnitude? -- Scott Reverse name to reply
Clyde <iloveuppercaseEREKOSE@EMAIL.IT> writes:

> pgw wrote: > >> How it look without windowing? > > They are different, and the problem still remains. > >> You can try to generate one harmonic signal and do FFTW (without >> windowing) on it and then show us the result. > > Here is the result: > http://i45.photobucket.com/albums/f77/yesitookmypills/FFT.png > (blue is FFTW, red is Matlab FFT) > > I used fftw_plan_r2r_1d(dim, in[0], out[0], FFTW_R2HC, FFTW_ESTIMATE) > and i compute the magnitude in this way: > magnitude[i] = sqrt(out[i]*out[i]); > > Here's the data with a small script in matlab: > > http://rapidshare.com/files/25643904/FFTWvsFFTMatlab.zip.html
One difference I can see is that you use abs() in Matlab but your own magnitude function in your other code (is it C?). Are you using doubles in the C code? Can we see your C code? -- % Randy Yates % "Watching all the days go by... %% Fuquay-Varina, NC % Who are you and who am I?" %%% 919-577-9882 % 'Mission (A World Record)', %%%% <yates@ieee.org> % *A New World Record*, ELO http://home.earthlink.net/~yatescr
Scott Seidman <namdiesttocs@mindspring.com> wrote in
news:Xns99109030D964Dscottseidmanmindspri@130.133.1.4: 

> Clyde <iloveuppercaseEREKOSE@EMAIL.IT> wrote in news:461e738c$0$4786 > $4fafbaef@reader4.news.tin.it: > >> fftw_plan_r2r_1d > > > Don't you need to deal with the "half complex" format that fftw is > asked to produce before you deal with the magnitude? >
See http://www.fftw.org/fftw3_doc/Real_002dto_002dReal-Transform- Kinds.html#Real_002dto_002dReal-Transform-Kinds -- Scott Reverse name to reply