PDA

View Full Version : Strange FFT Behavior in MATLAB


patrickjennings
09-20-2007, 01:01 PM
I think we can all agree that FFT { x*(n) } = X*(N-k), or the FFT of th
conj of x(n) is the conj of the reversed version of the FFT of x(n).

But in MATLAB if a= [1+2j 3+4j 5+6j 7+8j] then

fft(conj(a)) = [16+20j -8+0j -4-4j 0-8j]

and

conj(fliplr(fft(a))) = [-8+0j -4-4j 0-8j 16+20j]

Any ideas? Another engineer and I spent most of a day looking at th
model before finding the fundamental problem.


Cheers

/Patrick

fatnbafan
09-20-2007, 03:04 PM
This is due to the fact that MATLAB index starts at 1 while when you say

FFT{x*(n)} = X*(N-k),

you are assuming that k = 0,...,N-1.

Note that when you do X(N-k) when k = 0, the first sample is X(N), whic
is the same as X(0). However, if you use fliplr, this is no longer true.

HTH,

>I think we can all agree that FFT { x*(n) } = X*(N-k), or the FFT of the
>conj of x(n) is the conj of the reversed version of the FFT of x(n).
>
>But in MATLAB if a= [1+2j 3+4j 5+6j 7+8j] then
>
>fft(conj(a)) = [16+20j -8+0j -4-4j 0-8j]
>
>and
>
>conj(fliplr(fft(a))) = [-8+0j -4-4j 0-8j 16+20j]
>
>Any ideas? Another engineer and I spent most of a day looking at the
>model before finding the fundamental problem.
>
>
>Cheers
>
>/Patrick
>
>
>
>

Andor
09-20-2007, 03:19 PM
Patrick wrote:
> I think we can all agree that FFT { x*(n) } = X*(N-k), or the FFT of the
> conj of x(n) is the conj of the reversed version of the FFT of x(n).

It's not exactly the reversed version. Your formula is correct, you
just have to look at it closer. Consider

X[k] := DFT{ x[n] }
X1[k] := DFT{ x*[n] }

So for k=0 you get

X1[0] = X*[N-0] = X*[N]. However, this index is beyond the limit
(indices go from k=0,1,...N-1), so X*[N] is in fact X*[0] (all indices
are modulo N). If you write it out you get

(X1[0], X1[1], ... , X1[N-1]) = (X*[0], X*[N-1], ...., X*[1]).

Presumably Matlab's 1-based indexing doesn't make this matter any
easier.

> But in MATLAB if a= [1+2j 3+4j 5+6j 7+8j] then
>
> fft(conj(a)) = [16+20j -8+0j -4-4j 0-8j]
>
> and
>
> conj(fliplr(fft(a))) = [-8+0j -4-4j 0-8j 16+20j]

For that a, I get

X = fft(a) = [16+20j -8 -4-4j -8j]
X1 = fft(conj(a)) = [16-20j 8j -4+4j -8 ]

which is just what the formula says.

> Any ideas? Another engineer and I spent most of a day looking at the
> model before finding the fundamental problem.

No wonder. If I were to look all day at a model*, I would never find
any problems!

Regards,
Andor

*Favourite model to look at: Claudia.

patrickjennings
09-20-2007, 05:14 PM
Ah ha! X(N-k) isn't a "reversal" of X, it's a circular shift... that wa
a silly mistake, thanks for the quick answers all. If anyone is ever i
Dallas, TX beers are owed.

Cheers.

/ptj

>Patrick wrote:
>> I think we can all agree that FFT { x*(n) } = X*(N-k), or the FFT o
the
>> conj of x(n) is the conj of the reversed version of the FFT of x(n).
>
>It's not exactly the reversed version. Your formula is correct, you
>just have to look at it closer. Consider
>
>X[k] := DFT{ x[n] }
>X1[k] := DFT{ x*[n] }
>
>So for k=0 you get
>
>X1[0] = X*[N-0] = X*[N]. However, this index is beyond the limit
>(indices go from k=0,1,...N-1), so X*[N] is in fact X*[0] (all indices
>are modulo N). If you write it out you get
>
>(X1[0], X1[1], ... , X1[N-1]) = (X*[0], X*[N-1], ...., X*[1]).
>
>Presumably Matlab's 1-based indexing doesn't make this matter any
>easier.
>
>> But in MATLAB if a= [1+2j 3+4j 5+6j 7+8j] then
>>
>> fft(conj(a)) = [16+20j -8+0j -4-4j 0-8j]
>>
>> and
>>
>> conj(fliplr(fft(a))) = [-8+0j -4-4j 0-8j 16+20j]
>
>For that a, I get
>
>X = fft(a) = [16+20j -8 -4-4j -8j]
>X1 = fft(conj(a)) = [16-20j 8j -4+4j -8 ]
>
>which is just what the formula says.
>
>> Any ideas? Another engineer and I spent most of a day looking at the
>> model before finding the fundamental problem.
>
>No wonder. If I were to look all day at a model*, I would never find
>any problems!
>
>Regards,
>Andor
>
>*Favourite model to look at: Claudia.
>
>

robert bristow-johnson
09-20-2007, 05:56 PM
On Sep 20, 8:01 am, "patrickjennings" <[email protected]>
wrote:
> I think we can all agree that FFT { x*(n) } = X*(N-k),

actually the folks at The Math Works do not agree. as fatnbafan said,
MATLAB is hard-wired or hard-coded so that the indices of all arrays
begin with 1, not 0 as it should for the DFT or FFT.

so in MATLAB, if N=length(x); y = conj(x); X = fft(x); Y = =
fft(y); then

Y(k+1) = conj( X(mod(N-k+1, N)) ); % for 0 <= k < N

or stated so elegantly that it's amazing we all don't just sing the
praises of MATLAB,

Y(k) = conj( X(mod(N-k+2, N) ); % for 1 <= k <= N

gee, isn't that elegant?

r b-j

> or the FFT of the
> conj of x(n) is the conj of the reversed version of the FFT of x(n).
>
> But in MATLAB if a= [1+2j 3+4j 5+6j 7+8j] then
>
> fft(conj(a)) = [16+20j -8+0j -4-4j 0-8j]
>
> and
>
> conj(fliplr(fft(a))) = [-8+0j -4-4j 0-8j 16+20j]
>
> Any ideas? Another engineer and I spent most of a day looking at the
> model before finding the fundamental problem.
>
> Cheers
>
> /Patrick

Scott Seidman
09-20-2007, 06:10 PM
robert bristow-johnson <[email protected]> wrote in
news:[email protected] ups.com:

> actually the folks at The Math Works do not agree. as fatnbafan said,
> MATLAB is hard-wired or hard-coded so that the indices of all arrays
> begin with 1, not 0 as it should for the DFT or FFT.


It's an index-- just a conceptual place holder, not a definition. X(1),
in a ones indexed language, merely points to the first memory location in
array X.

Most likely the ones indexing in Matlab has its roots in Fortran.

Aside from having to remember what language you're working in and getting
bitten every now and again when moving back and forth between languages,
it really isn't that problematic.


> Y(k) = conj( X(mod(N-k+2, N) ); % for 1 <= k <= N
>
> gee, isn't that elegant?

As compared to what? If you're comparing to finding and downloading the
c implemenation of fftw, compiling it, linking to it, then figuring out
which of the many arcane fftw routines you need to be calling, and then
figuring out what arcane pattern your returned arrays are in, maybe
creating complex variable types from two doubles, and then first doing
your elegant zero-indexed circular shift, its a veritable tango. If
you're talking about doing the same in some zero-indexed matlab clone,
probably not as elegant.


--
Scott
Reverse name to reply