fft-ish notes, see https://gist.github.com/pervognsen/55a999540ead46b73732f24d75295216

fft-ish.txt

note: j is used instead of i in these notes as it's much easier to differentiate from 1 in my editor's font :)

So, let's consider a 4-point (N=4) DFT of the following signal:

x = [0, 1, 0, 1]

for regular DFT, we have 4 basis vectors:

 - k=0: [ 1,  1,  1,  1] (exp(-j(2pi/N)kn) for k=0 yields exp(0) or 1 for all n)
 - k=1: [ 1, -j, -1,  j] (exp(-j(2pi/N)kn) for k=1 yields exp(-j(pi/2)n))
 - k=2: [ 1, -1,  1, -1] (exp(-j(2pi/N)kn) for k=2 yields exp(-j(pi)n))
 - k=3: [ j, -1, -j,  1] (exp(-j(2pi/N)kn) for k=3 yields exp(-j(3pi/2)n))

Working that out, our expected result is:

 - X[0] = sum([ 0,  1,  0,  1]) =  2
 - X[1] = sum([ 0, -j,  0,  j]) =  0
 - X[2] = sum([ 0, -1,  0, -1]) = -2
 - X[3] = sum([ 0, -1,  0,  1]) =  0

so,

X = [ 2,  0, -2,  0]

Assuming a scaling factor of N, this intuitively makes sense.
By inspection, our original signal can be described as:

x[n] = 0.5 - 0.5 * cos(pi * n)

and this is _exactly_ what X describes; two partials, one at DC and one at nyquist, with each 0.5 term multiplied by
N (4) to have magnitude 2; thus this (non-normalized) DFT result appears to be correct.

Now let's consider computing the DFT by means of divide-and-conquer, by considering separately the even/odd samples,
and combining their spectra, as in Per's notes. First, note that the bin corresponding to the nyquist in the original
spectrum (N/2=2) is nonzero, which means that we expect there to be aliasing in the even/odd spectra. Hopefully, this
cancels once we phase-shift the odd spectrum - does it?

First, let's calculate the even/odd spectra, denoted as E and O, respectively.

E is trivial. It corresponds to e, the even samples of x, which are all 0 in this example. Thus, E is also all 0, or:

E = [0, 0]

We also have that o <-> O, and o is defined as the odd samples of x, or o = [1, 1]. The basis vectors for a 2-point DFT
are simply [1, 1] and [1, -1] (corresponding to DC and nyquist, respectively), so this yields:

O = [2, 0]

Or, pure DC, with a scaling factor of N/2=2, as expected... erm, mostly as expected - where's the alias? We know that
with frequency content at nyquist in the original signal, that when decimating by a factor of 2, our new nyquist is
half that of the original, so we'd expect the frequency content at the larger nyquist to be reflected about the new
nyquist, which should yield additional frequency content at DC. So, it's possible that the DC component that we see
here actually includes _both_ the original signal's DC component, and an alias of the original signal's nyquist
component. Can we verify/prove this?

To combine the two spectra into X, we know that:

X[k] = E[k/2] + O[k/2] * w[k]

where w[k] represents the phase rotation ("twiddle factor") corresponding to a 1-sample delay for bin k that we've
implicitly factored out of the odd frequency component computations when separating out even/odd spectra like this.

W represents (negative) powers of the Nth root of unity, which we can denote as U_N. Thus:

w[k] = U_N^-k for 0 <= k < N.

For N=4, this yields:

W = [ 1, -j, -1,  j]

Now, we can determine X by combining E, O, and W as follows:

X[0] = E[0] + O[0] * W[0] = 0 + 2 *  1 =  2
X[1] = E[1] + O[1] * W[1] = 0 + 0 * -j =  0
X[2] = E[0] + O[0] * W[2] = 0 + 2 * -1 = -2
X[3] = E[1] + O[1] * W[3] = 0 + 0 *  j =  0

so,

X = [ 2,  0, -2,  0]

which we, as know from before, is correct.

So, what about the alias from the frequency content at nyquist in x? From above, we can see that O[0] ends up
distributed over both the DC and nyquist bins in X (that is, X[0] and X[2]), so indeed, that bin contained DC as
well as the aliased nyquist component that we wanted to see.

With proper normalization, we would see that:

E_norm = [  0, 0] (no change)
O_norm = [  1, 0] (components divided by N/2)
X_norm = [1/2, 0, -1/2, 0] (components divided by N)

This suggests that the first component of O_norm is 0.5 DC and 0.5 the alias of the cosine wave, which again,
is what we expect.