Chronology Current Month Current Thread Current Date
[Year List] [Month List (current year)] [Date Index] [Thread Index] [Thread Prev] [Thread Next] [Date Prev] [Date Next]

Re: Fourier transform normalization



Corrected version. This morning's version was off by a factor of 1000 in a
couple of places. Publish in haste, repent at leisure....


Suppose we have a vector Vin, representing a voltage waveform sampled at
successive points in time. We can take the Fourier transform
Vf = fft(Vin)
and we can take the inverse transform:
Vout = ifft(Vf)

In a typical implementation, such as Scilab or Matlab, or in typical
textbooks, the normalization is as follows: Suppose that Vin is a delta
function, that is:
Vin = [1, 0, 0, ..... , 0]
then Vf will be all ones:
Vf = [1, 1, 1, ..... , 1]

But on the other hand, if Vf were a delta function in the frequency domain,
Vout would not be all ones; it would be smaller by a factor of N, where N
is the number of samples.

So.......... This leads to some questions:

1) Why are the forward transform and the inverse transform so
different? Conceptually they ought to be very, very symmetric. They need
to differ by the sign of (i) somewhere, but for real-valued inputs that
shouldn't matter very much.

2) What are the units on Vin?

3) What are the units on Vf?

4) We can call Vin and Vf the ordinates (in the time domain and in the
frequency domain, respectively). What are the corresponding abscissas?

5) Whatever happened to Parseval's theorem? I would hope that our signals
would uphold the law
norm(Vin) = norm(Vf) = norm(Vout)
so I guess we need to ask how to take the norm of such vectors.

===============================================

Here are some answers.

A) In this example, Vin has units of volts. Vf has units of volt-seconds.

B) The FFT and IFFT assume that the input abscissas are evenly spaced,
starting at zero -- but they doesn't keep track of the spacings. We have
to keep track of the spacings ourself. Suppose we have two seconds of
waveform data, sampled at intervals of
dt = 1 millisecond
for a total of
N = 2000 samples.

Then the abscissas in the frequency domain will have a spacing of
df = 1 / (N dt) (as a rule)
= 0.5 Hz (in this case)

C) I claim that when doing a forward FFT, the output should be scaled by a
factor of dt. In particular, if Vin is a delta function
Vin = [1, 0, 0, 0, .... , 0]
then the transform should be
myVf = myFFT(Vin)
= [.001, .001, .001, .... , .001]

D) I claim that when doing an inverse FFT, the output should be scaled by a
factor of df. In particular, if Vf is a delta function,
Vf = [1, 0, 0, 0, .... , 0]
then the inverse transform should be
myVout = myIFFT(Vf)
= [.5, .5, .5, .... , .5]

Note that when using the Scilab or Matlab inverse FFT, the software stuffs
in a factor of 1/N automagically, so you need to undo that (multiply by a
factor of N) before applying the factor of df. By the way, note that
multiplying by (N df) is the same as dividing by (dt).

E) Note that myIFFT is the inverse of myFFT. Their definitions are
completely symmetric; one has a factor of dt and the other has a factor of
df. Neither of them needs to have a factor of 1/N. This all works out
cleanly because
df * dt = 1/N

F) Parseval's theorem: The relevant norm-squared is not just
sum_i Vin_i * Vin_i i.e. Vin dot Vin
nor
sum_f Vf_f * Vf_f i.e. Vf dot Vf

Think about the Riemann integral: You don't just sum up a bunch of
ordinates; you need to weight each term in the sum by the _measure_ of the
corresponding abscissas.

Therefore the right answer is:
normsquared(Vin) = Vin dot Vin * dt ## units of volt^2 seconds
normsquared(Vf) = Vf dot Vf * df ## units of volt^2 seconds

and (big surprise :-) if you use myFFT and myIFFT these two quantities turn
out to be equal, which means we are upholding Parseval's theorem.

G) The same arrangements generalize cleanly to the case of Fourier
integrals; they are not limited to the discrete Fourier sums considered
above. The only additional trick is that you need to think of frequencies
(f) in circular measure (such as Hz == cycles per second) instead of
frequencies (omega) in angular measure (radians per second).

That is, in the inverse transform, be sure to write
Vout(t) = integral Vf(f) exp(2 pi i f t) df
rather than
Vout(t) = integral Vf(omega) exp(i omega t) d(omega)

... because the latter is off by a nasty factor of 2pi.