Something astray in my FFT understanding

Hi again all

As you might know from my other thread Im implementing an FFT from scratch
to get a better understanding of the whole time/frequency mapping. However
Ive run into a bit of a roadblock. I have implemented small bit of code
that takes N samples and performs a N point FFT. My code is running fine,
it gives the same result as when I do it out by hand, so whats at fault is
that my understanding of some stage of the transform is incorrect. I was
hoping that someone here might be able to spot where Im going wrong.

Rather that post up code or type out a long explanation of how I currently
understand the FFT to work, I thought that posting some examples of the
output Im getting might be a better idea. Hopefully theyll show symptoms
of a common trap or something similar.

The FFT is real, DIT and in-place.

Impulse	
-------
Time	Frequency
10	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10
0	10

Not sure if its working for an Impulse. Ive read the FFT of an impulse is
the impulse magnitude in every bin (in which case the code does work for
this input), but shouldn't the imaginary part be all zeros?

Impulse Shifted x[1]=10
-----------------------	
Time	Frequency
0	10
10	13.06455197
0	17.06825181
0	22.29890628
0	29.13252199
0	38.06033475
0	49.72412212
0	64.96233775
0	-10
0	-13.06455197
0	-17.06825181
0	-22.29890628
0	-29.13252199
0	-38.06033475
0	-49.72412212
0	-64.96233775

This is wrong as the FFT of an impulse shifted by one sample should be, in
the realpart, a cosine of frequency 1, and in the imaginary part, a sine of
frequency 1. Instead, Im getting exponentials...

Impulse Shifted x[2]=10	
-----------------------
Time	Frequency
0	10
0	14.1421345
10	19.99999683
0	28.28426452
0	-10
0	-14.1421345
0	-19.99999683
0	-28.28426452
0	10
0	14.1421345
0	19.99999683
0	28.28426452
0	-10
0	-14.1421345
0	-19.99999683
0	-28.28426452

Again weird exponential results.

One way I look at the results is if you consider the entire output array
as being the real component of the FFT result, I seem to be getting the
correct frequency but the 'wave' is the wrong 'shape'. This means that the
real part is twice as long as it should be and the 'wave shape' which
should be a cosine (for the real part) is instead concave shaped. The
imaginary part seems to have been chewed up and digested somewhere inside
the routine.

However Im not sure how to map these observations back to my understanding
of the FFT to find out where Im going wrong so any ideas or suggestion
would be really helpful.

Thanks in advance

Hi,

I am not sure what you mean by frequency here. The output of the DFT
should be complex valued. Can you post the complex values for the
routine so we can take a look ?

SeanOg skrev:
> Hi again all
>
> As you might know from my other thread Im implementing an FFT from scratch
> to get a better understanding of the whole time/frequency mapping. However
> Ive run into a bit of a roadblock. I have implemented small bit of code
> that takes N samples and performs a N point FFT. My code is running fine,
> it gives the same result as when I do it out by hand, so whats at fault is
> that my understanding of some stage of the transform is incorrect. I was
> hoping that someone here might be able to spot where Im going wrong.
>
> Rather that post up code or type out a long explanation of how I currently
> understand the FFT to work, I thought that posting some examples of the
> output Im getting might be a better idea. Hopefully theyll show symptoms
> of a common trap or something similar.
>
> The FFT is real,

What does that mean? That the FFT accepts only real-valued inputs?
That it returns only the real part of the DFT?

> DIT and in-place.
>
> Impulse
> -------
> Time	Frequency
> 10	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
> 0	10
>
> Not sure if its working for an Impulse. Ive read the FFT of an impulse is
> the impulse magnitude in every bin (in which case the code does work for
> this input), but shouldn't the imaginary part be all zeros?

Yes, it should. With the caveat that the impulse is in time bin 0 and
allowing for numerical inaccuracies.

> Impulse Shifted x[1]=10
> -----------------------
> Time	Frequency
> 0	10
> 10	13.06455197
> 0	17.06825181
> 0	22.29890628
> 0	29.13252199
> 0	38.06033475
> 0	49.72412212
> 0	64.96233775
> 0	-10
> 0	-13.06455197
> 0	-17.06825181
> 0	-22.29890628
> 0	-29.13252199
> 0	-38.06033475
> 0	-49.72412212
> 0	-64.96233775
>
> This is wrong as the FFT of an impulse shifted by one sample should be, in
> the realpart, a cosine of frequency 1, and in the imaginary part, a sine of
> frequency 1. Instead, Im getting exponentials...

Could this have to do with the FFT being real? If not, check out if the

exponentials in th ereal and imaginary parts might have been switched.
I am pretty sure this has to do with the phase term in the twiddle
factors.

> Impulse Shifted x[2]=10
> -----------------------
> Time	Frequency
> 0	10
> 0	14.1421345
> 10	19.99999683
> 0	28.28426452
> 0	-10
> 0	-14.1421345
> 0	-19.99999683
> 0	-28.28426452
> 0	10
> 0	14.1421345
> 0	19.99999683
> 0	28.28426452
> 0	-10
> 0	-14.1421345
> 0	-19.99999683
> 0	-28.28426452
>
> Again weird exponential results.
>
> One way I look at the results is if you consider the entire output array
> as being the real component of the FFT result, I seem to be getting the
> correct frequency but the 'wave' is the wrong 'shape'. This means that the
> real part is twice as long as it should be and the 'wave shape' which
> should be a cosine (for the real part) is instead concave shaped. The
> imaginary part seems to have been chewed up and digested somewhere inside
> the routine.

A bad sign... could it be that the FFT algorithm is based on a
real-valued
input, which is guaranteed to produce a symmetric spectrum, and
pops out half the number of real and imaginary coefficients?

> However Im not sure how to map these observations back to my understanding
> of the FFT to find out where Im going wrong so any ideas or suggestion
> would be really helpful.

Rune

Many thanks for the quick replies. Sorry for the confusion, Im probally
using the terminology a bit to quickly without thinking exactly what it
means.

What my routine does is it accepts an array of N points. This array
represents the signal to be transformed. This is what I meant about the
FFT being real, it accepts one array of real, not complex, numbers, each
point in the array representing one sample.

The output of the routine is an array of N points which represents (at
least what I believed up to now) the entire DFT; the real part of the DFT,
stored in the first half, and the imaginary part, stored in the second
half.

So using an impulse function of 8 points as an example (input is a single
8 point array);

Input                     Output
1 0 0 0 0 0 0 0 -> FFT -> 1 1 1 1 0 0 0 0
                          |real | |img  |

(the above output is what it should be, I get all 1's)


>Hi,
>I am not sure what you mean by frequency here. 

Again bad use of a term on my part. What I was getting at was that in the
example of the impulse shifted by two, the output looked like the output
of the 1 point shifted impulse but with 'twice the frequency'. 

>I am pretty sure this has to do with the phase term in the twiddle
>factors.

Im checking that out at the moment. I have a niggling in the back of my
head that something about how I consider how the phase factors should be
calculated isnt 100% correct.


>A bad sign... could it be that the FFT algorithm is based on a
>real-valued
>input, which is guaranteed to produce a symmetric spectrum, and
>pops out half the number of real and imaginary coefficients?

Ah, thats news to me. Looks like Ive got a very mixed up view of the whole
real/complex FFT. I always thought that; 

x[N] -> FFT -> X[N]

where x[N] is an array of real numbers (ie no complex component)

      first half X[N] is real part of resulant DFT
      second half X[N] is imaginary part of resultant DFT


Seems I need to go back over this topic. I think then the problem is the
FFT Ive implemented is specifcally for complex inputs or something...Gotta
hate when you think youve finally understood something you regard as
difficult and then realise that youre still as confused as ever :P

thanks again, its great to have a soundboard that I can use to figure out
these things rather muttering incoherent abuse at my books. Ill let you
know how I get on.

SeanOg skrev:
> Many thanks for the quick replies. Sorry for the confusion, Im probally
> using the terminology a bit to quickly without thinking exactly what it
> means.
>
> What my routine does is it accepts an array of N points. This array
> represents the signal to be transformed. This is what I meant about the
> FFT being real, it accepts one array of real, not complex, numbers, each
> point in the array representing one sample.

OK...

> The output of the routine is an array of N points which represents (at
> least what I believed up to now) the entire DFT; the real part of the DFT,
> stored in the first half, and the imaginary part, stored in the second
> half.

Maybe I'm just splitting a few hairs here, but we need to agree on the
details.
An N-point DFT produces N complex-valued coefficients, that contain
one real and one imaginary part, 2N in all. Due to a symmetry
when the input signal is real, namely

   X(w) =  re{X(w)} + j*im{X(w)}
   X(-w) = conj(X(w)) = re{X(w)} - j*im{X(w)}

the full spectrum is known from N/2 complex-valued coefficients.
So if your FFT routine accepts only real-valued signals and produces
N real numbers, it is almost certain that N/2 of those numbers
represent
the real part and the other N/2 the imaginary parts...

> So using an impulse function of 8 points as an example (input is a single
> 8 point array);
>
> Input                     Output
> 1 0 0 0 0 0 0 0 -> FFT -> 1 1 1 1 0 0 0 0
>                           |real | |img  |
>
> (the above output is what it should be, I get all 1's)

...like what you indicate here.

>
> >Hi,
> >I am not sure what you mean by frequency here.
>
> Again bad use of a term on my part. What I was getting at was that in the
> example of the impulse shifted by two, the output looked like the output
> of the 1 point shifted impulse but with 'twice the frequency'.
>
> >I am pretty sure this has to do with the phase term in the twiddle
> >factors.
>
> Im checking that out at the moment. I have a niggling in the back of my
> head that something about how I consider how the phase factors should be
> calculated isnt 100% correct.

Try restrict this to N = 2^M, M integer, first. That's how most FFTs
are
designed. Only after that one works, does one extend to general N.

> >A bad sign... could it be that the FFT algorithm is based on a
> >real-valued
> >input, which is guaranteed to produce a symmetric spectrum, and
> >pops out half the number of real and imaginary coefficients?
>
> Ah, thats news to me. Looks like Ive got a very mixed up view of the whole
> real/complex FFT. I always thought that;
>
> x[N] -> FFT -> X[N]
>
> where x[N] is an array of real numbers (ie no complex component)

Nope. See below.

>       first half X[N] is real part of resulant DFT
>       second half X[N] is imaginary part of resultant DFT

Only if the FFT routine is specifically designed to accept only
real-valued sequences. In general, you get N complex-valued
coefficients.

>
> Seems I need to go back over this topic. I think then the problem is the
> FFT Ive implemented is specifcally for complex inputs or something...Gotta
> hate when you think youve finally understood something you regard as
> difficult and then realise that youre still as confused as ever :P

Been there, done that. Not with the FFT (never attempted that one
from scratch, but you've got me thinking...), but still. Never mind,
all that is forgotten when you, in the end, know that you finally got
it right...

> thanks again, its great to have a soundboard that I can use to figure out
> these things rather muttering incoherent abuse at my books. Ill let you
> know how I get on.

Rune

SeanOg wrote:

...
> Input                     Output
> 1 0 0 0 0 0 0 0 -> FFT -> 1 1 1 1 0 0 0 0
>                           |real | |img  |

A real FFT of N points should return N/2+1 real values and N/2-1
imaginary values, ie.

1 0 0 0 0 0 0 0 -> FFT -> 1 1 1 1 1 + j[ (0) 0 0 0 (0) ]

The terms in brackets are the imaginary part at DC and Nyquist
frequency - they are always equal to zero (because of the Hermitian
symmetry for real inputs) and therefore not explicitely mentioned.

Sorry I didn't read the whole post, the above point was just an obvious
blunder that I saw.

Regards,
Andor

SeanOg wrote:
> As you might know from my other thread Im implementing an FFT from scratch
> to get a better understanding of the whole time/frequency mapping.

Sean,

Let me make a more basic suggestion.  If your goal is to understand the
time-frequency mapping, don't implement an FFT.  Implement the DFT
directly from the definition (i.e. the "slow" O(N^2) algorithm).

An FFT is just an optimized version of the DFT, and doesn't tell you
anything especially new about the properties of this transformation,
except in a computer-science sense.  An FFT just does the DFT faster.
The problem is that an FFT is more complicated to implement, and
obscures the underlying DFT definition, and it sounds like this is your
problem: you are not seeing the forest for the trees here.

For example, the DFT intrinsically takes N complex inputs and produces
N complex outputs.  If the inputs are real, then there is a redundancy
in the outputs that you can exploit to save a factor of 2 in
time/storage, but again this is an optimization...it comes at the price
of making things more complicated, so it is *not* a good place to
start.

Once you understand the DFT from its definition, then you should just
download a fast FFT code off the web and then it will compute exactly
the same thing much more quickly.

(FFT algorithms form a fascinating study in themselves, intersecting
lots of interesting mathematical areas from number theory to polynomial
algebras.  You should certainly dive in if you are interested, but be
warned that the study of FFT algorithms is actually fairly independent
of signal processing or of understanding the properties of time vs.
frequency representations.)

Cordially,
Steven G. Johnson

stevenj@alum.mit.edu wrote:
> SeanOg wrote:
> > As you might know from my other thread Im implementing an FFT from scratch
> > to get a better understanding of the whole time/frequency mapping.
>
> Sean,
>
> Let me make a more basic suggestion.  If your goal is to understand the
> time-frequency mapping, don't implement an FFT.  Implement the DFT
> directly from the definition (i.e. the "slow" O(N^2) algorithm).
>
> An FFT is just an optimized version of the DFT, and doesn't tell you
> anything especially new about the properties of this transformation,
> except in a computer-science sense.  An FFT just does the DFT faster.
> The problem is that an FFT is more complicated to implement, and
> obscures the underlying DFT definition, and it sounds like this is your
> problem: you are not seeing the forest for the trees here.
>
> For example, the DFT intrinsically takes N complex inputs and produces
> N complex outputs.  If the inputs are real, then there is a redundancy
> in the outputs that you can exploit to save a factor of 2 in
> time/storage, but again this is an optimization...it comes at the price
> of making things more complicated, so it is *not* a good place to
> start.

You can also construct a DFT operator with no complex math at all.
It's just the correlation against N real and othogonal vectors (1 DC
sum, N/2-2 cosines, and N/2-1 sines, the missing sine wave being all
zeros).  A complex DFT just moves the sine correlations into some
imaginary components, then duplicates these complex values to give
a redundant result (for an all real imput).

> Once you understand the DFT from its definition, then you should just
> download a fast FFT code off the web and then it will compute exactly
> the same thing much more quickly.

Or you can try factoring and unrolling a small DFT operator using
various trig identities until you get something that breaks into
O NLog(N) operations instead of N^2.  Then you will have some
understanding of what's behind an FFT.  The fast FFT code often
obscures the mathematical idea which lots of implementation
specific details (address arithmetic tricks, pipeline scheduling,
cache blocking, etc.) well beyond the DSP idea.


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

Many thanks to everyone for their comments, theyre all greatly
appreciated.

I went back over the DFT (as Steve suggested) and as was pointed out by
Rune and others, the problem was my understanding of the whole
real/complex issue. What I had actually implemented was some weird hybrid
of a complex/real FFT which included some of my own custom innovations, in
others words, a very confused routine! The symptom of using different types
of code to try and understand the process and not making sure I understood
the theory behind it first. Mea Culpa.

Im now much clearer on the whole issue (the reason why I always implement
things Im trying to understand, code doesnt let you away with
misconceptions).

Ive redone my routine and its now very much better behaved, giving nearly
the expected results, however there is a glitch which I havent been able
to fix or identify the source of, hopefully people wont mind me poking
around their grey matter yet again...

The routine is a DIT, in-place FFT accepting 2^K (K is an integer) complex
samples.
The size of the fft is equal to the number of input samples
No windowing

INPUT
-----
It accepts two arrays;
real[] representing the real part of the input
imaginary[] representing the imaginary part of the input

Im inputting signals that only have a real component (ie all imaginary[]
elements are assigned 0)

OUTPUT
------
real[] contains the real part of the transform
imaginary[] contains the imaginary part of the transform

The Glitch (going to use two examples)
----------
If I input a COSINE wave of F frequency, I expect (ignoring symmetry):

1. real[] - peak at element/bin F, zeros everywhere else
2. imaginary[] - zeros everywhere

My routine gives 1. but for the imaginary part, I get a small peak at F
and a number of other much smaller peaks


If I input a SINE wave of F frequency, I expect (ignoring symmetry):

1. real[] - zeros everywhere
2. imaginary[] - peak at element/bin F, zeros everywhere else

My routine gives 2. but for the real part, I get a small peak at F and a
number of other much smaller peaks.

Again Im hoping this is a 'common' mistake or suggests some mistakes Im
making, I understand its hard to make suggestion from little information.

Once I get this working properly, Ill send everyone the price of a beer
over Paypal!!

Thanks

PS thought i better mention that the unexpected peaks are more than just
internal computational inaccuracies, theyre of the magnitude of +/- 3E+00