Introduction to Discrete Fourier Transforms

APAM 1601y
Columbia University

Introduction

The Fourier transform expands functions in terms of sine and cosine functions of different wavelengths or frequencies. The Fourier transform is like "spectrum analyzer." It takes a function and returns a spectrum: a number proportional to the fraction of the function that is represented by each wavelength or frequency.

The "discrete Fourier transform" became in use in the 1960's as the world entered the  "computer age," Today, information is primarily represented (and transmitted) in digital form, and the field of "digital signal processing" is at the core of all forms of digitally-recorded music, video, and information.

This notebook introduces the discrete Fourier transform (DFT) and illustrates its use. We also describe the recursive "fast Fourier transform" (FFT) that is built into Mathematica.

John W. Tukey (1915-2000)

John Wilder Tukey was born in New Bedford, Massachusetts, on June 16, 1915.  He earned bachelor's and master's degrees in chemistry at Brown University in 1936 and 1937, respectively. He then went to Princeton, where, in 1938 he was a Jacobus Fellow, and upon receiving his doctorate in mathematics from Princeton in 1939, he was appointed Henry B. Fine Instructor in Mathematics.  A decade later, at age 35, he was became a full professorship.  He directed Princeton's Statistical Techniques Research Group from its founding in 1956.   At the same time, he was a Member of Technical Staff at AT&T Bell Laboratories in 1945, advancing to Assistant Director of Research, Communications Principles in 1958 and, in 1961, to Associate Executive Director-Research Information Sciences, which position he held until his retirement in 1985.

John Tukey has attracted international attention for his studies in mathematical and theoretical statistics and their applications to a wide variety of scientific and engineering disciplines.  He has led the way in the now-burgeoning fields of Exploratory Data Analysis and Robust Estimation, and his contributions to the Spectrum Analysis of Time Series and other aspects of Digital Signal Processing have been widely used in engineering and science.  He has been credited with coining the word "bit", a contraction of "binary digit", which refers to a unit of information, often as processed by a computer.

In addition to strong continuing interests in a wide variety of areas of statistical philosophy, techniques and application, Tukey has been active in improving the access of the scientist to the scientific literature, particularly through the development of citation and permutation indices to the literature of statistics and probability. In another area collaboration with a fellow mathematician resulted in the formation of the Cooley-Tukey Fast Fourier Transform (FFT) algorithm, a mathematical technique that greatly simplifies computation for Fourier series and integrals.  Other interests range through applications to such fields asbehavioral sciences, geophysics and pharmaceutical research.

In 1965, Tukey was named the first recipient of the Samuel S. Wilks Award of the American Statistical Association.  He received the National Medal of Science in 1973 "for his studies in mathematical and theoretical statistics…and for his outstanding contributions to the applications of statistics to the physical, social, and engineering sciences".  Tukey received the Shewhart Medal of the American Society for Quality Control in 1977; the Institute of Electronic and Electrical Engineers’1982 Medal of Honor , for the Cooley-Tukey Fast Fourier Transform (FFT) Algorithm; and the American Society for Quality Control’s Deming Medal in 1983.  In 1989 he was elected a Foreign Member of The Royal Society (London) .  Princeton University honored him in 1984 with the James Madison Medal , given annually to an alumnus of the Graduate School who has had a distinguished career, who has advanced the cause of graduate education, or who has achieved a record of outstanding public service. In 1989 he received the Monie A. Ferst Award of Sigma Xi, and in 1990 the Educational Testing Service Award for Distinguished Service to Measurement.

Tukey has taught on both the undergraduate and graduate levels and is widely sought as a seminar leader and lecturer.  He holds honorary degrees from Case Institute of Technology, the University of Chicago and Brown, Temple and Yale Universities.

President of the Institute for Mathematical Statistics in 1960,he has also served as Vice President of the American Statistical Association and as Chairman of both its Biometrics Section and its Section on Physical and Engineering Sciences.He has been a Member of the Council of the Biometric Society,and served on the NationalResearch Council for more than a decade.He is a member of the National Academy of Sciences,and is a member of the American Philosophical Society (vice president,1975-77) and the American Academy of Arts and Sciences,Honorary Fellow of the Royal Statistical Society (London),and a Foreign Member of The Royal Society (London).Among his other affiliations are the International Statistical Institute,the American Association for the Advancement of Science (vice president,Section A,1972 and Section U,1974) and Sigma Xi.

(From the Bell Labs Memorial site: <http://cm.bell-labs.com/cm/ms/departments/sia/tukey/>.)

Euler's Identity for (e^(Iθ))

ExpToTrig[Exp[I θ]]

TrigToExp[Cos[θ]]

TrigToExp[Sin[θ]]

Plot[{Re[(E^(I θ))],Im[(E^(I θ))]},{θ,0,2 π},PlotStyle->{Blue,Red}]

Defining the Discrete Fourier Transform

Mathematica has a powerful command the discrete Fourier transform (DFT) that can be used for digital signal processing.

??Fourier

??InverseFourier

Example (64 Samples)

nSize = 2^6

x=Table[N[Sin[((2 π (i-1) 2)/nSize)]],{i,1,nSize}];
ListPlot[x]

h[n_] := Table[Exp[2 π I (i - 1)n/nSize], {i, 1, nSize}];

h[0]

h[1]

h[63]

h[-2] - h[nSize - 2]

h[3] . Conjugate[h[3]]// N  // Chop

h[3] . Conjugate[h[4]]// N  // Chop

Table[h[i] . Conjugate[h[i]] // N // Chop, {i, 0, nSize-1}]

Table[h[i] . Conjugate[h[10]] // N // Chop, {i, 0, nSize-1}]

xf = (x.#)& /@ Table[h[i], {i, 0, nSize-1}]// Chop (* DFT *)

DFT Orthogonal Modes

ListPlotComplex[x_List,options___] := ListPlot[{Re[x], Im[x]},PlotStyle-> {Blue,Red},options]

ListPlotComplex[h[3]]

modeTable=Table[ListPlotComplex[h[i 2+j],PlotLabel->"n = "<>ToString[i 2+j], Joined -> True],{i,0,3},{j,0,1}];
Show[GraphicsGrid[modeTable]]

modeTable=Table[ListPlotComplex[h[i 2+j],PlotLabel->"n = "<>ToString[i 2+j], Joined -> True],{i,13,16},{j,0,1}];
Show[GraphicsGrid[modeTable]]

modeTable=Table[ListPlotComplex[h[i 2+j],PlotLabel->"n = "<>ToString[i 2+j], Joined -> True],{i,28,31},{j,0,1}];
Show[GraphicsGrid[modeTable]]

DFT (Long Way)

dft[x_List] := Module[{n,e0, eVec},
n = Length[x];
e0 = Exp[2 π I/n];
eVec = Table[e0^(i-1), {i, 1, n}];
Table[(x .eVec^(i-1)), {i, 1, n}]/Sqrt[n]];

idft[xf_List] := Module[{n,e0, eVec},
n = Length[x];
e0 = Exp[-2 π I/n];
eVec = Table[e0^(i-1), {i, 1, n}];
Table[(xf .eVec^(i-1)), {i, 1, n}]/Sqrt[n]];

xf = dft[x];

xf // Chop

ListPlotComplex[xf, PlotRange-> All, PlotStyle -> PointSize[0.02]]

xt = idft[xf];

ListPlot[Chop[xt]]

Print["Max transform error = ",Max[Abs[xt-x]]];ListPlot[Abs[xt-x],PlotRange->All,PlotStyle->PointSize[0.02]]

Example ("in-between" mode)

nSize = 2^6

x=Table[N[Sin[((2 π (i-1) 3.4)/nSize)]],{i,1,nSize}];
ListPlot[x,PlotLabel->"3.4 Periods Long"]

xf=dft[x];
ListPlot[Abs[xf],PlotRange->All,PlotStyle->PointSize[0.02`],PlotLabel->"3.4 Periods Long"]

When the sampled signal, x[n], have oscillations that do not "fit" precisely into the sample period, then the Fourier spectrum is not simply related to the primary oscillations.

Example (Shift)

nSize=(2^6);
x1=Table[N[Sin[((2 π (i-1) 3)/nSize)]],{i,1,nSize}];
x2=RotateRight[x1,2];
p1=ListPlot[x1];
p2=ListPlot[x2,PlotStyle->Red];
Show[p1,p2]

x1f = dft[x1];
x2f = dft[x2];

ListPlotComplex[x1f, PlotRange -> All, PlotLabel -> "No Shift"]
ListPlotComplex[x2f, PlotRange -> All, PlotLabel -> "Shift"]

x1f // Chop

x2f // Chop

A shift in the time sequence results in a phase-shift in the frequency representation.

dx=idft[x1f-x2f];ListPlot[Re[dx],PlotLabel->"Shifted Difference"]

The difference of a shifted signal results in a derivative. (Of course!)

Example (1k Samples)

nSize=(2^10);
s1=Table[Sin[(((i-1) 2 π 6.3)/nSize)]+Sin[(((i-1) 2 π 15.7)/nSize)],{i,1,nSize}];
ListPlot[s1]

(sf1 = dft[s1];) // Timing

(sf2 = Fourier[s1];) // Timing

Max[Abs[sf1 - sf2]]

Our "home-made" DFT is numerically equivalent to Mathematica's built-in Fourier[…]. However, our "home-made" DFT is very, very slow. This is because Fourier uses the "fast-Fourier transform" algorithm (FFT) discussed next.

Recursive FFT (J. W.  Cooley and J. W. Tukey, 1965)

The following expressions illustrate the recursive algorithm of the FFT.

Clear["s"];
nSize = 2^3;
sData = Table[s[i], {i, 1, nSize}]

(Sqrt[2] dft[sData] // Simplify)// TableForm

Transpose[Partition[sData, 2]]

{fts1, fts2} = (dft /@ Transpose[Partition[sData, 2]]) // Simplify

(sfExample1 = ( (fts1 ~ Join ~ fts1) + Table[Exp[2 π I/nSize]^(m-1), {m ,1, nSize}] (fts2 ~ Join ~ fts2) ) // Simplify)//TableForm

(sfExample2 = dft[sData]  Sqrt[2]// Simplify) // TableForm

sfExample1 -  sfExample2// Simplify

Recursive FFT

fftRecursive[x_List] := Module[{n,xeven, xodd, g, w},
If[(n = Length[x]) == 1, Return[x]];
{xeven, xodd} = fftRecursive /@ Transpose[Partition[x,2]];
g = Exp[2 π I/n];
w = Table[g^(m-1), {m, 1, n}];
(xeven ~ Join~ xeven) + w (xodd ~ Join ~ xodd)]

fft[x_List] := (* defaullt normalization *)fftRecursive[x]/Sqrt[Length[x]]

(sf3 = fft[s1]; ) // Timing

ListPlotComplex[sf3, PlotRange-> All];

Max[Abs[sf3 - sf2]]

A Larger Example

nSize = 2^18 (* 256k Samples *)

rate = 41000 (* samples per second *)

nSize/rate  // N(* total time, Second *)

frequency = 2000.0; (* cycles per second, Hz *)
rate/frequency (* samples per cycle *)

g[i_] := Sin[2π (i -1) (frequency/rate)] Sin[2 π (i - 1) 2/nSize]^2

x=Table[N[g[i]],{i,1,nSize}];
ListPlot[Take[x,{1,Length[x],100}]]

ListPlay[x, SampleRate -> rate]

(xf = Fourier[x]; )//Timing

ListPlot[Log[Abs[Take[xf,{1,Length[xf],100}]]],PlotRange->All]

frequencies = Table[i rate/nSize, {i, 0, nSize/2 }] ~Join~ Reverse[Table[-i rate/nSize, {i, 1, nSize/2-1}]] // N;

xff = Transpose[{frequencies, xf}];

xff // Shallow

ListPlot[Apply[{#1,Log[Abs[#2]]}&,Take[xff,{1,Length[xff],100}],{1}],PlotRange->All,PlotLabel->"Frequency Spectrum Hz"]

Summary

The discrete Fourier transform (DFT) is fundamental to modern digital signal processing. The transform takes as an input a sequence of equally-spaced samples of a signal and produces on output the complex frequency spectrum of the signal.

The DFT of a N samples requires N dot-products of the sample-list with a Fourier modal vector. The number of computations grow as (N^2)!

With the "discovery" of the recursive FFT algorithm, the DFT can be performed with many fewer operations. The number of operations scales like N Log[N]. The FFT of one million samples can be performed in about 2 seconds. Our direct algorithm would take about 2 days!

The FFT (and similar transforms) have enabled the recent revolution in digital communications, media, and storage.