Re: fftshift help in octave 2.1.35, incorrect?

[Roberto Hernandez]

Mirek Kwasniak wrote:
> On Sat, Apr 13, 2002 at 03:40:58PM -0500, Hugo van der Merwe wrote:
> [...]
>
> > - Function File:  fftshift (V)
> >     Perform a shift of the vector V, for use with the `fft' and `ifft'
> >     functions, in order the move the frequency 0 to the center of the
> >     vector or matrix.
> >
> >     If V is a vector of E elements corresponding to E time samples
> >     spaced of Dt each, then `fftshift (fft (V))' corresponds to
> >     frequencies
> >
> >          f = linspace (-E/(4*Dt), (E/2-1)/(2*Dt), E)
>
>
> Yes, it's a bug :)

You're right. Fortunately it's only a bug in the help text, the function 
works well. Have you sent the patch to bug-octave?

> > The way I understand the fourier transform, I would believe that the
> > highest frequency you would get out of an fft would be half the sampling 
> > frequency. E.g. if you sample at 44.1kHz, from an fft you should get a 
> > spectrum from -22.05kHz to nearly 22.05kHz.
>
>
> Yes

True for even number of samples. If the number of samples is odd, the 
spectrum samples would go from -a to a, where a is a positive number 
such that a < samplerate/2 (notice that's a strict less than). In fact, 
for E odd, a = (E-1)/(2*E) * samplerate.

I guess the best way to look at it is the following. The distance 
between frequency samples (f0) is given by (for E odd or even) f0=fs/E, 
where fs is the sampling frequency and E is the total number of samples. 
So the samples will be located at frequencies: ..., -2*fs/E, -fs/E, 0, 
fs/E, 2*fs/E, ... If E is an odd number then there is no integer n such 
that (n*fs/E = fs/2).

> > I would suggest
> >
> >  samplerate = 1/Dt
> >  f = -samplerate/2 : samplerate/E : samplerate/2-samplerate/E
>
>
> Non-interger step in loop is a classical bug (see any beginners book for
> programming) :(
>
> I always use a formula:
>
>   f = (-E/2:E/2-1)/E/Dt
>
> I assume that E is even (i don't know how calculate it if E is odd).

If E is odd, the expression would be:

f = ( -(E-1)/2 : (E-1)/2 ) / E/Dt

The following will give you the frequency vector, independent of whether 
E is odd or even:

f = (0:E-1) / E/Dt;
a = (f >= 1/(2*Dt));
f(a) = f(a) - 1/Dt;
f = fftshift(f);

HTH,
Roberto



-------------------------------------------------------------
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
-------------------------------------------------------------