Notes 3, Computer Graphics 2, 15-463Fourier Transforms and theFast Fourier Transform (FFT) AlgorithmPaul HeckbertFeb. 1995Revised 27 Jan. 1998We start in the continuous world; then we get discrete.Definition of the Fourier TransformThe Fourier transform (FT) of the function f (x) is the function F(ω), where:F(ω) =Z∞−∞f (x)e−iωxdxand the inverse Fourier transform isf (x) =12πZ∞−∞F(ω)eiωxdωRecall that i =√−1andeiθ= cosθ +i sinθ.Think of it as a transformation into a different set of basis functions. The Fourier trans-form uses complex exponentials (sinusoids) of various frequencies as its basis functions.(Other transforms, such as Z, Laplace, Cosine, Wavelet, and Hartley, use different basisfunctions).A Fourier transform pair is often written f (x) ↔ F(ω),orF( f (x)) = F(ω) where Fis the Fourier transform operator.If f (x) is thought of as a signal (i.e. input data) then we call F(ω) the signal’s spectrum.If f is thought of as the impulse response of a filter (which operates on input data to produceoutput data) then we call F the filter’s frequency response. (Occasionally the line betweenwhat’s signal and what’s filter becomes blurry).1Example of a Fourier TransformSuppose we want to create a filter that eliminates high frequencies but retains low frequen-cies (this is very useful in antialiasing). In signal processing terminology, this is called anideal low pass filter. So we’ll specify a box-shaped frequency response with cutoff fre-quency ωc:F(ω) =1 |ω|≤ωc0 |ω| >ωcWhat is its impulse response?We know that the impulse response is the inverse Fourier transform of the frequencyresponse, so taking off our signal processing hat and putting on our mathematics hat, all weneed to do is evaluate:f (x) =12πZ∞−∞F(ω)eiωxdωfor this particular F(ω):f (x) =12πZωc−ωceiωxdω=12πeiωxixωcω=−ωc=1πxeiωcx− e−iωcx2i=sinωcxπxsince sinθ =eiθ−e−iθ2i=ωcπsinc(ωcπx)where sinc(x) = sin(πx)/(πx). For antialiasing with unit-spaced samples, you want thecutoff frequency to equal the Nyquist frequency, so ωc= π.Fourier Transform PropertiesRather than write “the Fourier transform of an X function is a Y function”, we write theshorthand: X ↔ Y.Ifz is a complex number and z = x +iy where x and y are its real andimaginary parts, then the complex conjugate of z is z∗= x − iy. A function f (u) is even iff (u) = f (−u),itisodd if f (u) =−f (−u), it is conjugate symmetric if f (u) = f∗(−u),and it is conjugate antisymmetric if f (u) =−f∗(−u).2discrete ↔ periodicperiodic ↔ discretediscrete, periodic ↔ discrete, periodicreal ↔ conjugate symmetricimaginary ↔ conjugate antisymmetricbox ↔ sincsinc ↔ boxGaussian ↔ Gaussianimpulse ↔ constantimpulse train ↔ impulse train(can you prove the above?)When a signal is scaled up spatially, its spectrum is scaled down in frequency, and viceversa: f(ax) ↔ F(ω/a) for any real, nonzero a.Convolution TheoremThe Fourier transform of a convolution of two signals is the product of their Fourier trans-forms: f ∗ g ↔ FG. The convolution of two continuous signals f and g is( f ∗ g)(x) =Z+∞−∞f (t)g(x −t) dtSoR+∞−∞f (t)g(x −t) dt ↔ F(ω)G(ω).The Fourier transform of a product of two signals is the convolution of their Fouriertransforms: fg↔ F ∗ G/2π.Delta FunctionsThe (Dirac) delta function δ(x) is defined such that δ(x) = 0 for all x 6= 0,R+∞−∞δ(t) dt = 1,and for any f (x):( f ∗δ)(x) =Z+∞−∞f (t)δ(x − t) dt = f (x)The latter is called the sifting property of delta functions. Because convolution with a deltais linear shift-invariant filtering, translating the delta by a will translate the output by a:f (x) ∗δ(x −a)(x) = f (x −a)3Discrete Fourier Transform (DFT)When a signal is discrete and periodic, we don’t need the continuous Fourier transform.Instead we use the discrete Fourier transform, or DFT. Suppose our signal is anfor n =0...N −1, and an= an+jNfor all n and j. The discrete Fourier transform of a, also knownas the spectrum of a,is:Ak=N−1Xn=0e−i2πNknanThis is more commonly written:Ak=N−1Xn=0WknNan(1)whereWN= e−i2πNand WkNfor k = 0...N −1arecalledtheNth roots of unity. They’re called this because, incomplex arithmetic, (WkN)N= 1 for all k. They’re vertices of a regular polygon inscribedin the unit circle of the complex plane, with one vertex at (1, 0). Below are roots of unityfor N = 2, N = 4, and N = 8, graphed in the complex plane.W42ReImN=2W20W21N=4W40W43W411−1 −11i-iW84N=8W80W86W82−11i-iW87W85W83W81Powers of roots of unity are periodic with period N, since the Nth roots of unity arepoints on the complex unit circleevery 2π/ N radiansapart, and multiplyingby WNis equiv-alent to rotation clockwise by this angle. Multiplication by WNNis rotation by 2π radians,that is, no rotation at all. In general, WkN= Wk+jNNfor all integer j. Thus, when raising WNto a power, the exponent can be taken modulo N.The sequence Akis the discrete Fouriertransform of the sequence an. Each is a sequenceof N complex numbers.The sequence anis the inverse discrete Fourier transform of the sequence Ak.Thefor-mula for the inverse DFT isan=1NN−1Xk=0W−knNAk4The formula is identical except that a and A have exchanged roles, as have k and n.Also,the exponent of W is negated, and there is a 1/N normalization in front.Two-point DFT (N=2)W2= e−iπ=−1, andAk=1Xn=0(−1)knan= (−1)k·0a0+ (−1)k·1a1= a0+ (−1)ka1soA0= a0+a1A1= a0−a1Four-point DFT (N=4)W4= e−iπ/2=−i,andAk=3Xn=0(−i)knan= a0+(−i)ka1+(−i)2ka2+(−i)3ka3= a0+(−i)ka1+(−1)ka2+ika3soA0= a0+a1+a2+a3A1= a0−ia1−a2+ia3A2= a0−a1+a2−a3A3= a0+ia1−a2−ia3This can also be written as a matrix multiply:A0A1A2A3=11111 −i −1 i1 −11−11 i −1 −ia0a1a2a3More on this later.To compute A quickly, we can pre-compute common subexpressions:A0= (a0+a2) + (a1+a3)A1= (a0−a2) − i(a1−a3)A2= (a0+a2) − (a1+a3)A3= (a0−a2) + i(a1−a3)5This saves a lot of adds. (Note that each add and multiply here is a complex (not real) op-eration.)If we use the following diagram for a complex multiply and
View Full Document
Unlocking...