Download Notes on Numerical Analysis and Computing - Lecture notes 15 | MATH 541 and more Study notes Mathematics in PDF only on Docsity!
Trigonometric PolynomialsApplications of the FFT Numerical Analysis and Computing^ Lecture Notes #15— Approximation Theory —The Fast Fourier Transform, with Applications^ Peter Blomgren,^ 〈[email protected]
Department of Mathematics and Statistics^ Dynamical Systems GroupComputational Sciences Research Center^ San Diego State UniversitySan Diego, CA 92182-7720^ http://terminus.sdsu.edu/^ Fall 2009 Peter Blomgren, 〈[email protected]〉^ The Fast Fourier Transform, w/Applications
— (1/38)
Trigonometric PolynomialsApplications of the FFT Outline 1 Trigonometric Polynomials Trigonometric Interpolation: Introduction Historical Perspective^ Ã^ the FFT 2 Applications of the FFT Recap, Notes, & Historical Perspective 1080p Video, ..., Image Processing Peter Blomgren, 〈[email protected]〉^ The Fast Fourier Transform, w/Applications
— (2/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Trigonometric Polynomials: Least Squares
⇒^ Interpolation. Last Time:^ We used trigonometric polynomials,
i.e.^ linear combinations of the functions:^ ^ Φ(x)^0
(^1) = (^2) Φ(x) = cos(kx),^ k^ = 1,... ,k n Φ(x)^ =^ sin(kx),^ n+k^
k^ = 1,... ,^ n^ −^1 to find^ least squares^ approximations (where
n^ <^ m) to equally spaced data (2m^ points) in the interval [
−π, π], at the node points x=^ −π^ + (jπ/m),^ j^ = 0j^
,^1 ,... ,^ (2m^ −^ 1). This Time:^ We will find the
interpolatory^ (n^ =^ m) trigonometricpolynomials... and we will figure out how to do it fast! Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (3/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Why use Interpolatory Trigonometric Polynomials?^ Interpolation of large amounts of
equally spaced^ data by trigonometric polynomials produces very good results (close tooptimal,^ c.f.^ Chebyshev interpolation). Some Applications^ Digital Filters (Lowpass, Bandpass, Highpass)^ Signal processing/analysis^ Antenna design and analysis^ Quantum mechanics^ Optics^ Spectral methods numerical solutions of equations.^ Image processing/analysis^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (4/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Interpolatory Trigonometric Polynomials^ Let^ xbe 2m^ equally spaced node points in [i^
−π, π], and^ f=^ f^ (x) thei^ i^ function values at these nodes. We can find a trigonometric polynomialof degree^ m:^ P(x)^ ∈ Twhich interpolates the data:m^ aa^0 m^ S(x) =^ +^ m^2
m−^1 ∑cos(mx) +^ [acos(kxk^ k= ) +^ bsin(kx)]^ ,k^ where^2 m−^1 ∑^1 a= k^ m^ j=
2 m−^1 ∑ (^1) fcos(kx) b= j j k m^ j= fsin(kx).j^ j^ The only difference in this formula compared with the one correspondingto the least squares approximation,
S(x),^ n^ <^ m^ is the division by twon of the^ acoefficient.m^ Where does the factor of 2 come from???^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (5/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Finding the Factor of 2: Breakdown of the Lemma^ The factor of two comes from the failure of the second part of the lemma(which we showed last time):^ Lemma^ If the integer r is not a multiple of
2 m, then 2 m− 12 m−^1 ∑∑ cos(rx) =^ sin(rx) = 0.j j^ j=0^ j= Moreover, if r is^ not^ a multiple of m, then^2 m−^1 ∑^ [cos(^ j=
2 m−^1 ∑ (^2 2) rx)]=^ [sin(rx)]=j j^ j= m. Now, since^ n^ =^ m, we end up with one instance (the cos(
mx)-sum)j^ where^ r^ =^ m, and the lemma fails.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (6/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Finding the Factor of 2: Computing the Sum^ Now, since we are interpolating, the basis functionΦ(x) = cos(mx) is part of the set. When we computem
〈Φ,^ Φ〉mm we need^2 m−^1 ∑^ [cos(mxj^ j=
2 m−^1 ∑ (^2) )]=^ [cos(−πm^ + j= jπ^2 m )]^ m^2 m− (^1) ∑^2 = [cos((j − m)π)] j=0 2 m− (^1) ∑2(j−m) = (−1) j=0 2 m− (^1) ∑ = 1 = 2 m. j= Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (7/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Historical Perspective^ Much of the analysis was done by
Jean Baptiste Joseph Fourier in the early 1800s, but the use of the Fourier series representationwas not practical until 1965.Why? The straight-forward implementation requires about
(^2) 4m operations in order to compute the coefficients
˜ ˜a, and b. In 1965^ Cooley and Tukey
published a 4-page paper titled
“An algorithm for the machine calculation of complex Fourier series”
in the journal^ Mathematics of Computation.
The paper describes an algorithm which computes the coefficients using only O(m^ logm) operations.^2 It is hard to overstate the importance of this paper!!! The algorithm is now known as the
“Fast Fourier Transform”
or just the^ “FFT”.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (8/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Comparing Operation Counts
1 of 2 (^2) m 4m 3m^ +^ m^ logm^ Speedup^2 16 1,^
112 9 64 16,^
576 28 256 262,^
2,816^93 1,024^ 4,194,
13,312^315 4,096^ 67,108,
61,440^ 1, 16,384^ 1,073,741,
278,528^ 3, 65,536^ 17,179,869,
1,245,184^ 13, 262,144^ 274,877,906,
5,505,024^ 49, 1,048,576^ 4,398,046,511,
24,117,248^ 182, 4,194,304^ 70,368,744,177,
104,857,600^ 671, 8,388,608^ 281,474,976,710,
218,103,808^ 1,290, 16,777,216^ 1,125,899,906,842,
452,984,832^ 2,485, Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (13/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Comparing Operation Counts
2 of 2 m^ 4m
2 3m^ +^ m^ logm^ Speedup^2 1,048,576^ 4,398,046,511,
24,117,^
8,388,608^ 281,474,976,710,
218,103,808^ 1,290,
m^ =^8 ,^388 ,^608 roughly corresponds to an 8-Megapixel cameraIf a 3.8 GHz Pentium chip could perform one addition ormultiplication per clock-cycle (which it can’t), we could computethe Fourier coefficients for the 8-Megapixel image in
FFT^ Slow FT 0.057 seconds^ 20.576 hours Each “FFT second” translates roughly to 15 “Slow-FT days.”^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (14/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Implementing the FFT^ Assigning implementation of the FFT falls in the category of
cruel and unusual punishment
Matlab implements fft^ the Fourier transform (data
→^ coefficients) the inverse Fourier transform (coefficients ifft
→^ data) and the helper functions fftshift^ Shifting (and unshifting the coefficient (mostly forifftshiftdisplay purposes)The 2-dimensional version
[i]fft2, and^ n-dimensional
[i]fftn version of the FFT are also implemented.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (15/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT The Fastest Fourier Transform in the West (FFTW)
1 of 2 [http://www.fftw.org/
]
FFTW is a C / Fortran subroutine library for computing theDiscrete Fourier Transform (DFT) in one or more dimensions, ofboth real and complex data, and of arbitrary input size.FFTW is free software, distributed under the GNU General PublicLicense.Benchmarks, performed on on a variety of platforms, show thatFFTW’s performance is typically superior to that of other publiclyavailable FFT software. Moreover, FFTW’s performance isportable: the program will perform well on most architectureswithout modification.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (16/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT The Fastest Fourier Transform in the West (FFTW)
2 of 2 It is difficult to summarize in a few words all the complexities thatarise when testing many programs, and there is no ”best” or”fastest” program.However, FFTW appears to be the fastest program most of thetime for in-order transforms, especially in the multi-dimensionaland real-complex cases. Hence the name, ”FFTW,” which standsfor the somewhat whimsical title of ”Fastest Fourier Transform inthe West.” Please visit the benchFFT[http://www.fftw.org/benchfft/
]^ home page for a more extensive survey of the results.The FFTW package was developed at MIT by Matteo Frigo andSteven G. Johnson Matlab uses FFTW these days... see
help fftw. Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (17/38)
Trigonometric Polynomials^ Trigonometric Interpolation: IntroductionApplications of the FFTHistorical Perspective
Ã^ the FFT Homework #9, Not Due Fall 2009
HW-Extra-For-Fun Read the matlab^ help^ for^ fft
,^ ifft,^ fftshift, and^ ifftshift
. [1] Let^ x=(-pi:(pi/8):(pi-0.1))’
. Let^ f1=cos(x),^ f2=cos(2*x),
f3=cos(7x), f4=cos(8x),^ f4=cos(9*x),^ g1=sin(x)
,^ g2=sin(2x). Compute the^ fft^ of these functions.[2] Let^ f=1+cos(2x). Compute^
fftshift(fft(f)). Let^ g=1+sin(3*x). Compute^
fftshift(fft(g)). [3] Use your^ observations^ from
[1]^ and [2]^ to^ construct^ a^
low-pass^ filter:^ Given a function^ f, compute the^ fft(f)
.^ For a given N, keep only the coefficients corresponding to the N lowest frequencies (set all the others to zero).
ifft^ the result.Let^ f=1+cos(x)+sin(2x)+5cos(7*x)
and apply the above with N=4. Plot both the initial and the filtered f.[4] Use the FFT to determine the trigonometric interpolating polynomial of degree 8^2 for^ f^ (x) =^ xcos^ x^ on [−π, π]. Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (18/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing^ FFT Applications^ Figure:^ Some additional reading? Peter Blomgren, 〈[email protected]〉^ The Fast Fourier Transform, w/Applications
— (19/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing The Fast Fourier Transform, Recap The Fast Fourier Transform (FFT) is an
O(m^ log(m)) algorithm^2 for computing the 2m^ complex coefficients
cin the discretek^ Fourier transform:^2 m−^1 ∑^1 S(x) = m^ m^ k=
(^2) ikx (^) ce, where c=k k m−^1 ∑ikπj/m^ fej^ j= whereas the straight-forward implementation of the sum would^2 require^ O(m) operations.We noted [last time] that for a problem of size 2
23 = 8,^388 ,^608
this reduced the computation time by a factor of a 2 weeks —
i.e. a computation can be done in
t^ seconds using the FFTs would require^ ∼^2 t^ weeks to complete using the “Slow” FT.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (20/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Communications From:^ www.spd.eee.strath.ac.uk
2 of 4 If we take simple digital pulse that is to be sent down a telephoneline, it will ideally look like this: If we take the Fourier transform of this to show what frequenciesmake up this signal we get something like:^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (25/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Communications From:^ www.spd.eee.strath.ac.uk
3 of 4 This means that the square pulse is a sum of infinite frequencies.However if the telephone line only has a bandwidth of 10 MHz thenonly the frequencies below 10 MHz will get through the channel.This will cause the digital pulse to be distorted
e.g. This fact has to be considered when trying to send large amountsof data down a channel, because if too much is sent then the datawill be corrupted by the channel and will be unusable.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (26/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Communications From:^ www.spd.eee.strath.ac.uk
4 of 4 Extending the example of the telephone line, whenever you dial anumber on a ” touch-tone” phone you hear a series of differenttones. Each of these tones is composed of two different frequenciesthat add together to produce the sound you hear.The Fourier transform is an ideal method to illustrate this, as itshows these two frequencies
e.g. Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (27/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Side Scan Sonar From:^ www.spd.eee.strath.ac.uk
1 of 3 Another use of DSP techniques (including DFT/FFT) is in sonar.The example given here is side-scan sonar which is a little differentfrom the normal idea of sonar.With this method a 6.5-kHz sound pulse is transmitted into theocean toward the sea floor at an oblique angle. Because the signalis transmitted at an oblique angle rather than straight down, thereflected signal provides information about the inclination of thesea floor, surface roughness, and acoustic impedance of the seafloor. It also highlights any small structures on the sea floor, folds,and any fault lines that may be present.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (28/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Side Scan Sonar From:^ www.spd.eee.strath.ac.uk
2 of 3 Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (29/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Side Scan Sonar From:^ www.spd.eee.strath.ac.uk
3 of 3 The image on the previous slide came from the United StatesGeological Survey (USGS) which is using geophysical Long RangeASDIC (GLORIA) sidescan sonar system to obtain a plan view ofthe sea floor of the Exclusive Economic Zone (EEZ). The pictureelement (pixel) resolution is approximately 50 meters. The dataare digitally mosaiced into image maps which are at a scale of1:500,000. These mosaics provide the equivalent of ”aerialphotographs” that reveal many physiographic and geologic featuresof the sea floor.To date the project has covered approximately 2 million squarenautical miles of sea floor seaward of the shelf edge around the 30coastal States. Mapping is continuing around the American FlagIslands of the central and western Pacific Ocean.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (30/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Astronomy From:^ www.spd.eee.strath.ac.uk
1 of 3 Venus is Earth’s closest companion, and is comparable in it’s sizeand diameter and yet scientists knew virtually nothing about it, asit is perpetually covered with a cloud layer which normal opticaltelescopes can’t penetrate, so Magellan (a satellite) had radar andadvanced Digital Signal Processing that was designed to ”see”through this cloud layer. It’s mission was map 70% of the planetwith radar and to reveal surface features as small as 250 metersacross.The black and white pictures that it sent back were strips of theplanets surface, about 20km wide, from the north pole to thesouth pole.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (31/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Astronomy From:^ www.spd.eee.strath.ac.uk
2 of 3 Figure:^ One example image of a surface feature called ”Pandora Corona.”If you lookat the image, you will see a two black lines through the picture. This is just amismatch between the strips sent back by Magellan. It also gives you an idea of thescale of the image as each strip is 20km wide.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (32/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Example: Image Processing — Highpass Filtering
4 of 5 −150 −100 −50 0 50 100 150 −100 −50 0 50 100 Figure:^ Example of high-pass filtering the image. We have filtered out the Fouriercoefficients corresponding to the lowest frequencies.^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (37/38)
Trigonometric Polynomials^ Recap, Notes, & Historical PerspectiveApplications of the FFT1080p Video, ..., Image Processing Example: Image Processing — Lowpass Filtering
5 of 5 −150 −100 −50 0 50 100 150 −100 −50^0 50 Figure:^ Example of low-pass filtering the image. We have filtered out the Fouriercoefficients corresponding to the highest frequencies. [See movies]^ Peter Blomgren,^ 〈[email protected]
〉^ The Fast Fourier Transform, w/Applications
— (38/38)
