Amplitude and Phase of a discrete Fourier Spectrum

A. Dieckmann                  ELSA, Physikalisches Institut der Universität Bonn

DFT_1.gif

DFT_2.gif

    This tutorial describes the calculation of the amplitude and the phase from DFT spectra with finite sampling. The DFT matrix is intuitively
understood as a set of probes, each sensitive for a certain frequency and corresponding phase information present in a sequential sample of data,
distributed in time or space.

Principle

    The multiplication of the complex n × n matrix (there exist several conventions regarding phase and normalization, all working from similar principles...)

DFT_3.gif

with an array of n real numbers DFT_4.gif yields a complex valued vector DFT_5.gif

DFT_6.gif

that can be interpreted as "Fourier Spectrum". It delivers information about the amplitudes, frequencies and phases of the basic cosines
spanning the function of which an interval having a length of integer number of periods is being sampled into the vector DFT_7.gif. The vector DFT_8.gif consists
of n equally spaced samples, where r is the ‘running’ sample counter. The index r may count samples taken progressing in time or space.
The first sample DFT_9.gif corresponds to the value of the sampled function at the left edge DFT_10.gif of the interval DFT_11.gif and the last sample DFT_12.gif is
taken at DFT_13.gif In this way we get n samples (not n + 1), separated by a distance tt/n.
    The complex exponential in the matrix above can be imagined as unit vector that rotates in the complex plane with ‘frequency’ (s-1) / tt  and
progressing in time or space DFT_14.gif  according to the nature of the interval tt. During the matrix multiplication the vector of angular orientation
corresponding to r is multiplied with the sample of value DFT_15.gif and the vector sum is calculated over all r, giving a complex number DFT_16.gif. This is repeated
for all s from 1 to n, each DFT_17.gif signaling the presence of a contribution in the samples oscillating with a frequency related to s.
The set of all DFT_18.gif represents then the spectrum. Note the cancellation of the interval tt inside the argument of the exponential. Depending on the
application ‘spatial’ variables k (wave number) and x (location) or ‘time’ variables f (frequency) and t (time) may be used.
    If the samples contain parts changing with the same period as the rotating unit vector this vector sum will exhibit a large magnitude,
because the biggest product vectors all point roughly in the same direction. If the changes of the samples are of a different period with respect to
the rotating unit vector the sums magnitude will be smaller. This will be shown in the next four pictures, their production code is also given.

DFT_19.gif

DFT_20.gif

    Above we see the rotating unit vector multiplied with a sampled sine function, so the magnitude of the vector is proportional to the
value of the sine. Both are ‘in phase’ (rotation and sine have the same period). The y-t plane of the data and the complex plane of the ‘rotator’
were overlayed. The sine and the vectors appear in the graphic to lie in the same coordinate system and show their relation in a clear way.
The vectors were drawn starting from the position of the sampled points. Negative samples let the rotation phase jump by π.
    All rotating vectors (complex numbers) are summed and the result (see the red line below, beginning in the origin (0,0) and ending
in the point pe (upper black point)) is large. It corresponds to the Fourier amplitude of frequency 1Hz.

DFT_21.gif

Graphics:vector sum

    Next the same function is sampled, but the unit vector is rotating more frequently (not once but four times!) within the period of the sine function.
They are ‘out of phase’ and the resulting vector sum will be small.

DFT_23.gif

DFT_24.gif

    Above we see the rotating unit vector multiplied with a sampled sine function being out of phase (different frequency). The red line
representing the sum of rotating vectors has completely disappeared behind the overlapping black points, see below. The scale is enlarged
compared to the previous vector sum plot (note the length of the edges).

DFT_25.gif

Graphics:vector sum

    The length of the red line is too close to zero to be displayed :

DFT_27.gif

DFT_28.gif

Reconstruction of different Functions from the Fourier Amplitude and Phase

    Now we consider more complicated functions. For example let's define a periodic function g[t]

DFT_29.gif

    With parameter values given in this example (m = 3) as three triples of amplitude (a), integer frequency (f), phase (φ) and constant DFT_30.gif

DFT_31.gif

and a sampling time interval, that covers a full period interval tt, calculated by 1 / (the greatest common divisor of all contributing frequencies),

DFT_32.gif

DFT_33.gif

the function and its cosines look like

DFT_34.gif

Graphics:Function composed of 3 cosines

    Nyquist's theorem tells us we need at least two samples per period of the highest component frequency ( n > DFT_36.gif=  2 tt * DFT_37.gif ).
The lower limit for n lies in this case at

DFT_38.gif

DFT_39.gif

and the vector containing the n samples DFT_40.giftaken within the interval [tt[ (“closed” to the left, “open” to the right), is calculated. The Drop command
is used to get the number of samples right :

DFT_41.gif

DFT_42.gif

DFT_43.gif

    The next graph shows the function g[t], the samples DFT_44.gif and the 'cosine components' of the function

DFT_45.gif

Graphics:Function composed of 3 cosines

    Recall that the first sample sits on the left edge of the interval, whereas the last sample is taken inside the interval near the right edge.
THIS IS IMPORTANT ! Only iff the point on the right edge of the sampling interval of the period is dropped, the correct reconstruction and
periodic continuation is ensured...The reconstructed function will repeat itself outside the interval tt  (this is the definition of periodic).
From the samples we get the complex Fourier array

DFT_47.gif

DFT_48.gif

and plot the absolute value  as spectrum :

DFT_49.gif

DFT_50.gif

Graphics:Fourier spectrum with dominant peaks

    The most relevant bins as well as their mirrors are readily identified (offset in s by 1, because the first term in the spectrum is the one
describing the zero frequency part of the function), all others being zero. Their corresponding frequencies are given by f = (s - 1) / tt  in Hertz,
where tt is the full period interval that has been sampled. Apart from the first point the spectrum shows mirror symmetry with respect to
the index DFT_52.gif = Ceiling[ (n+1)/2 ], where Ceiling is a function returning the next integer number ≥ its argument.

Why is the Fourier Array showing a Mirror Symmetry ?

    Whenever the phase angle φ of the exponential function in the Fourier matrix would be larger than π, 2π is subtracted and φ stays
in the range -π ≤ φ ≤ π. Instead of counting the phase angle by the left turn it is now (if φ > π) measured by the right turn (in the negative sense of rotation).
This operation leaves the matrix element, which is just a complex number, unchanged. Within the Fourier matrix the phase angle φ increases
with row number reaching the π border at the middle. There will then be pairs of rows with the same (row number) distance to the middle.
The entries of these paired rows will be complex conjugates of each other. So the upper half of the Fourier array, which is a complex vector,  
mirrors the lower half with negative phase. In the spectrum plot using absolute numbers we get with increasing s a sequence of amplitude values
up to the highest frequency DFT_53.gif, that is visible by the Fourier transformation and that lies in the middle of the frequency indices DFT_54.gif = 11 in the above plot).
Going further up in s we get the same sequence, but in reverse order.

Finding Amplitude and Phase for select frequencies

    The next statement collects all bin numbers of the spectrum with bin values over threshold into a list s. The threshold is used to suppress insignificant bins
(not necessary in this example which is very clean, but add some noise ...) :

DFT_55.gif

DFT_56.gif

    Now using the Fourier spectrum the offset and amplitudes are found (see next section) to be

DFT_57.gif

DFT_58.gif

DFT_59.gif

DFT_60.gif

and the phases are

DFT_61.gif

DFT_62.gif

from which we can perfectly reconstruct the original function :

DFT_63.gif

Graphics:Function and reconstructed function

    If the original amplitude is negative (in this example a[[3]]),  the corresponding phase is offset by π...

DFT_65.gif

DFT_66.gif

    If a phase lies outside  the interval (- π ≤ φ ≤ π) just add or subtract 2 π.

The Connection between Amplitude and Fourier Amplitude

    The proportionality factor of the amplitude depends in this normalization convention on the number of samples n.
It can be determined by applying the Fourier transformation explicitly to samples of one simple cosine DFT_67.gif + A cos[ 2π f t - φ]:

DFT_68.gif

    A look at every frequency s in the spectrum reveals only three non zero entries: The peak in the spectrum lies at s = f + 1 (f ∈ Integers),
its mirror at s = n - f +1 and the zero frequency term at s = 1 :

DFT_69.gif

DFT_70.gif

    The complex number at f + 1 (== Fourier bin) has magnitude ADFT_71.gif and phase φ. The original amplitude A is therefore obtained
by multiplication of the discrete Fourier amplitude with 2 / DFT_72.gif. All other bins in the lower half (s ≠ f + 1) are zero except the
zero frequency term (offset) which comes out as DFT_73.gif.

Conclusion

    Each entry DFT_74.gif (s ≠ 1) in the lower half of the Fourier spectrum represents a cosine of frequency corresponding to s with
amplitude DFT_75.gif and phase DFT_76.gif being present in the original samples and DFT_77.gifis DFT_78.gif. If there are not enough
sample points (Nyquist) or if the interval tt is too short, the analysis gets less precise...The sampled interval is expected to cover exactly
the period of the periodic function. If the function is meant to be non-periodic, ignore the reconstructed part outside the interval.

Additional Example

    As an additional example using many more data points consider a function displaying a periodic sequence of pulses defined by edge points:

DFT_79.gif

    This function 'pulse' is sampled within the periodic interval at distances d of the argument x into an array ps with n data points. The increment d
must be small enough to 'grasp the structure' of the pulse (Nyquist theorem...) :

DFT_80.gif

DFT_81.gif

DFT_82.gif

    Find the Fourier transformed array :

DFT_83.gif

    The spectrum of the above pulse looks in a logarithmic presentation like

DFT_84.gif

Graphics:Logarithmic Amplitude of Pulse Spectrum

    We then get a good reconstruction of the original rectangular pulse function from a superposition of cosine functions with
amplitudes and phases calculated as outlined above using the first half of the (symmetric) Fourier transform array. The first term
represents the 'zero frequency' part, the offset. The frequencies of the contributing cosines are again calculated from their position
in the Fourier array divided by the span of the periodic interval.  Also notice the 'Gibbs Phenomenon', the overshoot appearing
at the kinks present in the original function. Shown are different approximations (all terms (red), 9 terms (pink) and 2 terms (orange) of the
Fourier expansion, i.e. different frequency cuts in the spectrum, the black dashed line uses an amplitude cut – only important terms
larger than ampCut are taken into account, the actual position of ampCut in the spectrum is shown in the graphic above with another dashed line.

DFT_86.gif

DFT_87.gif

Graphics:Rectangular Pulses and Fourier Approximations

    Now with a lengthened sampling interval the period gets longer and the transformation looks differently :

DFT_89.gif

Graphics:Rectangular Pulses and Fourier Approximations

Other Aspects of a Fourier Transformation, Effect of a High Pass

    The FourierTransform of a Gaussian pulse is :

DFT_91.gif

DFT_92.gif

    As σ is a positive quantity  FourierTransform gives the same result :

DFT_93.gif

DFT_94.gif

    Plot a pulse:

DFT_95.gif

DFT_96.gif

    Get the spectrum of amplitudes as the absolute value of the Fourier Transform (as the distribution is again Gaussian
with width 1/σ –> we plot 3 times the width) :

DFT_97.gif

DFT_98.gif

    In the plot above also the real and imaginary part of the spectrum is shown, both are formed by a mixture of
frequency dependent amplitude and phase (real part=a(ω) cos(φ(ω)),imaginary part=a(ω) sin(φ(ω))). The value of μ
for example can be found from the period Δω (within frequency space) of real or imaginary part as μ=2π/Δω.
    The phase of the Fourier Transform is given by the imaginary part of the argument of the complex exponential divided by the
imaginary unit, it contains the information about the position µ of the pulse given as the slope of the line describing the phase as function of ω :

DFT_99.gif

DFT_100.gif

    Sample the Gaussian pulse. This function is assumed non periodic, so the sampling interval is chosen such that it covers
the part of the function, that should be analyzed for frequency content or be reconstructed.  But the drop of the last sample is still
necessary to get μ and σ  - the horizontal scale - right :

DFT_101.gif

DFT_102.gif

DFT is

DFT_103.gif

    The normalization DFT_104.gif makes FourierTransform and Fourier 'equal', jmax  corresponds to the 3 sigma width in the Gaussian spectrum :

DFT_105.gif

DFT_106.gif

DFT_107.gif

    The next phase plot shows the same line as in the phase plot above, because - π is equivalent to π (each time the value rises above π,
2π is subtracted...). The Chop on fgs eliminates numerical noise that lets phase values fluctuate from π to -π and backwards.

DFT_108.gif

DFT_109.gif

    Plot both reconstructed function and original function, and we see that the DFT continues the original interval as the period of the sampled function :

DFT_110.gif

DFT_111.gif

    The red points in the next plot represent  the spectrum multiplied with a high pass characteristic, which damps just the amplitude linearly from DFT_112.gif to 0
and leaves the phase untouched :

DFT_113.gif

DFT_114.gif

    Damp the spectrum at low frequencies:

DFT_115.gif

    The reconstructed pulse (red) shows undershoots :

DFT_116.gif

DFT_117.gif

similar for rectangular Pulse, Edge Detection, no low Frequencies

    Use same function as above:

DFT_118.gif

DFT_119.gif

DFT_120.gif

DFT_121.gif

    Find the transformed array :

DFT_122.gif

    Two basic (complex) transfer functions f(ω) are DFT_123.gif (special phase convention by FT in Mathematica...).
This time lets choose the high pass characteristic f(ω) = ω / (ω + i DFT_124.gif) with DFT_125.gif = 1 / (R*C) = 333 Hz :

DFT_126.gif

Graphics:Typical High Pass Filter,  ω  = 333 Hz                                   g

    Blue and red show amplitude and phase of the transfer function, green and yellow the real and imaginary part.
Now we calculate damped spectrum by binwise multiplication of the complex transfer function with the complex Fourier transformed array:

DFT_128.gif

DFT_129.gif

    The reconstructed pulse shows mainly the edges of the original, the low frequency part ( ' DC  part ' ) is suppressed and change is
emphasized in a differentiated signal:

DFT_130.gif

Graphics:Rectangular Pulses and Fourier Approximation (red) after high-pass filter

Afterthought

    The point of this tutorial is to demonstrate the working of the Fourier Transform to the newcomer and to show how to avoid certain pitfalls
by giving detailed examples, that may serve as orientation. Good luck!

Spikey Created with Wolfram Mathematica 8.0       Go to Download Page