help-octave
[Top][All Lists]

RE: Octave vs. Matlab vs. IDL

 From: Steve Goncalo Subject: RE: Octave vs. Matlab vs. IDL Date: Fri, 28 Sep 2001 08:39:45 -0400

Zero fill is how most fft algorithms pad out an odd sized array
to a power of two before performing the transform. That can do
some strange stuff, but your problem is a lot simpler than that.

Even for real valued inputs, the FFT can produce complex
data at its output. With a sine wave input like you used
the largest amplitude value is really imaginary.
In theory, the real part of that fft output should
have been zero, and all you were seeing was round-off
error in the finite length arithmetic. Order of operations
can significantly affect round off errors, so it would
be surprising if Matlab and Octave did give exactly the
same round off errors.

If you did the same thing with a cosine, the peak is
in the real part of the output.

In general, you probably want to be looking at the
absolute value of the complex data when trying to
find a peak response.

octave:28> % compute time series, sine and cosine values
octave:28> t=[0:(2*pi)/288:50*2*pi-(2*pi)/288]';
octave:29> sint=sin(t);
octave:30> cost=cos(t);
octave:31>
octave:31> fs=fft(sint); % complex fft out
octave:32> fc=fft(cost);
octave:33>
octave:33> [val,index]=max(abs(real(fs))) %peak in the real part of
fft(sint)
val =  8.1311e-11
index = 51
octave:34> [val,index]=max(abs(imag(fs))) %peak in the imag part
val = 7200.0
index = 51
octave:35> [val,index]=max(abs(fs)) %peak in the amplitude
val = 7200.0
index = 51
octave:36>
octave:36> [val,index]=max(abs(real(fc))) %peak in the real part of
fft(cost)
val = 7200.0
index = 51
octave:37> [val,index]=max(abs(imag(fc))) %peak in the imag part
val =  8.1387e-11
index = 51
octave:38> [val,index]=max(abs(fc)) %peak in the amplitude
val = 7200.0
index = 51
octave:39>
octave:39> wavet=cos(t+pi/4); % give sinusoid a 45 degree phase shift
octave:40> fw=fft(wavet);
octave:41> [val,index]=max(abs(real(fw))) %peak in the real part of
fft(wavet)
val = 5091.2
index = 51
octave:42> [val,index]=max(abs(imag(fw))) %peak in the imag part
val = 5091.2
index = 14351   %% If you graph this you will see two peaks of nearly the
%% same amplitude 50 cells in from the start and end of
%% the array of 14400 cells
octave:43> [val,index]=max(abs(fw)) %peak in the amplitude
val = 7200.0
index = 51
octave:44>

-----Original Message-----
Sent: Thursday, September 27, 2001 5:28 PM
To: Help-Octave
Subject: fft: Octave vs. Matlab vs. IDL

I've taken the FFT of a sine function defined by:

>> t = [0:(2*pi)/288:50*2*pi - (2*pi)/288]';
>> sint = sin (t);
>> f = real(fft(sint));

in MATLAB, Octave, and IDL.  All three give me different answers,
although they do give a peak in power at 50, as expected.

Matlab -> peak at 50 of ~  2e-11
Octave -> peak at 50 of ~ -8e-11 (yes that "-" is real)
IDL    -> peak at 50 of ~  2e-06 (PI here was single precision)

I am trying to use the FFT to calculate an auto-correlation.  I know
this sounds strange, but the very slight differences I get between
calculating auto-correlations this way, and by using something like
xcov.m, significantly affects the stability, and predictive accuracy of
the linear prediction filters I am trying to generate...at least when
testing it on a simple sine function.  As much as I hate to say it, I
suspect that the FFT from IDL is giving the most reasonable answer,
because it spits out a relatively smooth function, with only a single
peak at 50, whereas both Octave and Matlab spit out very noisy data,
with local peaks almost as large as the main one at 50!

Can somebody who knows more about FFTs than me suggest a reason why
these would all return different values (and more importantly, how might
I get them to return more similar answers)?  I might guess this has
something to do with padding my vectors with zeros in a certain way, but
that is only because I've heard people say this...I don't really
understand what that means, or how it affects the FFT algorithm.

-EJR

-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

reply via email to