On this page, I want to reorganise my FFT routines again. Why so?
Till now, my
FFT's correlate the input with functions over the interval [0, 1] of
the fundamental cycle, or [0, 2pi] in radians. For an impression of
that, here is a
plot of the first 5 harmonic cosines over that interval:
Phase angles run from zero till 2pi here. This can be inconvenient.
Spectrum coefficients are complex, and frequency bins are neighboring
branches of the complex plane. Theoretically, the branch cut can be
defined anywhere, but for practical purposes it is often defined
between pi and pi radians. That is on the negative real axis, from the
origin to minus infinity.
The complex arctangent function in the standard C library has that
range, to
mention one important consideration. At least for the forward
transform, it
would be better to
correlate over the interval [0.5, 0.5] of the fundamental cycle, and
thus have phase angles of the spectrum output between pi and pi
radians:
Since we have no negativesigned indexes in computer
implementations, it
seems that we need to right shift the functions over half an interval
or pi radians. That is what I thought at first sight. But ooopps, small
problem:
A timeshift has different phaseeffects for each frequency. The odd
harmonics must be shifted by an odd multiple of pi radians and the even
harmonics by an even multiple of pi radians. That means the even
harmonics should remain like they are, while the odd harmonics can be
multiplied with 1. Then we get this:
That looks better. For the sine parts a similar modification holds.
In the forward transform we normally have negative sine functions like
these:
Shifting the virtual zero point to halfway the interval requires the
even harmonic sine functions to be negative, and the odd harmonics
positive:
I have copied a Fourier matrix from an earlier page, and
cut&pasted a phasecentered version from it, which is on the right.
Recall that the W^{kn} factors are complex exponentials e^{i2*pi*k*n/N}.
It is essentially a matter of swapping the left and right halves of
the matrix. As you can see, the beautiful symmetry of the
original Fourier matrix is lost in the phasezerocentered version. The
rows
are no longer equal to the columns. Therefore, row permutations are
needed for the transpose in the inverse transform. That is illustrated
below. On the left is the normal transpose and on the right is the
permutated version.
It is also possible to center both time input, and frequency output,
in the forward transform. In that case, the matrix is it's own
transpose again, apart from the rotation of the complex exponentials in
the inverse transform, as is always the case. You have to get used to
the rotated spectrum then. The positive frequencies appear in the upper
half of the output.

Since we are specially concerned with the behaviour of the odd
harmonics, I want to sketch an impression of the impact on these. On
the left below, two harmonics are sketched in a ordered Fourier
situation, and on the right I centered the quarters according to input
and output. Remarkably, the cosines flip but the sines do not. Still,
taking the shifted spectrum into account, the relevant sine direction
is flipped as well.
Fortunately, we do not need elaborate mathematical formula's to
achieve these effects. There is only one mechanism at work. Multiplying
the input with (1)^{n} will center the output. And centering
the input will multiply the output with (1)^{n}. No matter
what input: time signal or spectrum. So at both ends of the transform
the same things happen. You can do both: center the input and multiply
it with (1)^{n}. You then get an output which is centered and
alternated as well. That is the effect of the Fourier matrix in the
last figure. I went through quite some mishaps with code typo's and
other errations before getting fully convinced of this simple truth.
The centering and alternations can be done outside the FFT routine,
but it is more efficient to integrate some signflippings in the first
FFT stage.
Below is an illustration of that in my familiar square paper style so you can compare with my earlier FFT pages. Multiplication of the input array with (1)^{n} can be implemented as an adaptation of the first stage add/subtract factor matrix. The original is on the left and the adaptation on the right.
Centering the signal input, a virtual shift of time zero to the
frame centre, is slightly less straigthforward. We have to signflip
the odd correlation vectors in the FFT matrix. Basically, all odd
harmonics in the FFT are 'inheritants' of the
fundamental rotation, which is also odd. One option, of many, is to
signflip the fundamental twiddle factors, both sine and cosine. All
other twiddle factors are even harmonics and must be left untouched.
Swapping the add/subtract factors is equal to signflipping the
twiddles, so this again can be integrated in the same matrix, as shown
on the left below. On the right is the combination of the two signflip
actions.
I found it even more convenient to do the modifications in an FFT
scheme for decimated input. The first stage of that scheme has no
twiddles, therefore it is advantageous to peel it anyway and do
normalisation there. In a decimated sequence, the odd numbers and the
highest half of the numbers in the array have their positions
exchanged.
On the page with bitreversal
tables you can check which numbers go where in case of bitreversed
permutation. It happens that similar signflipping actions can be done
for this FFT scheme. On the left is the normal add/subtract factor
matrix, and on the
right the version that effectuates centering of the input.
To compensate for the phase/timezero centering, the modification
showed on
the left below could be used in the inverse transform. Or both
modifications can be combined in forward and inverse transform alike,
as shown on the right.
Notice that an FFT for decimated input has a reversed sequence of
operations, as compared to an FFT for decimated output.
Four different FFT radix 2 schemes are illustrated on my page Radix 2 FFT Forms.
Let us now check the output of a forward FFT with time zero
centered. I want to see if I get a
familiar pattern when I enter a pulse at x[(N/2)+1], now being time
one:
Indeed it is the same output that I had from a pulse at x[1] with my
noncentered FFT. Apparently, the simple signflip does it's job well.
I am curious how a pulse at x[33] would transform with a
noncentered FFT, for comparison. Here it is:
The above plot illustrates the alternating spectrum, as it results
from timeshifted input.
Now I return to the time zero centered input, and I want to combine
with a frequency zero centered output. Harmonic zero (DC) is at index
N/2 in this case. The spectrum is rotated. The benefit of this
transform is a better symmetry between forward and inverse. But you
have to remind this index offset, which may be inconvenient. I prefer
not to make regular use of this representation.
Now that I know how to implement a virtual timezero and/or
frequencyzero centering, one
dangling question remains. I have learned or at least concluded
sometimes that multiplying a signal with (1)^{n} effectuates a
'rotation over pi radian' of the spectrum. On this page, I have
carefully avoided that word, or replaced it by 'centered spectrum'.
What does it actually mean, a rotated spectrum? In the plot above, the
spectrum can be considered as being shifted rightward, while the part
with negative frequencies is shifted in from the left. If you would
make a circle of the xaxis over the range of the framesize, the shift
would be a rotation over half the framesize. A circle... can we draw
the spectrum on the complex plane? That is hard. To draw both phases
and frequencies would require different branches of the complex plane
in one figure.
Yet I want to produce some visualisation. Let me draw real and
imaginary parts of a frequency as points with different colours. Then
rotate the whole thing over half a circle. A low frequency becomes a
high frequency. You get the alias(es). In addition to that, the
negative frequencies become the positive frequencies and vice versa.
If you would flip the signs of the imaginary parts to conjugate the
spectrum, after rotation, you would get exactly the inverse of the
original spectrum. Why would one create an inverse of a spectrum? It is
an important aspect of Wavelet theory to create orthogonal halfband
filters. I was just curious how I could visualize that, though it is
quite offtopic here. In Wavelets technique, they do alternateflip the
filter impuls response: reverse the signal's direction and multiply
with (1)^{n}. This is almost the equivalent to rotation and
conjugation of the filter kernel's frequency response. Except for that
onesample discrepancy you get when timereversing a discrete signal.
(See the page Invertible FFT).