1How to Multiplyintegers, matrices, and polynomialsCOS 423Spring 2007slides by Kevin WayneConvolution and FFTChapter 303Fourier AnalysisFourier theorem. [Fourier, Dirichlet, Riemann] Any periodic functioncan be expressed as the sum of a series of sinusoids.sufficiently smoothtN = 1N = 5N = 10N = 100! y(t) = 2"sin ktkk=1N#4Euler's IdentitySinusoids. Sum of sines and cosines.Sinusoids. Sum of complex exponentials.eix = cos x + i sin xEuler's identity5Time Domain vs. Frequency DomainSignal. [touch tone button 1]Time domain.Frequency domain.! y(t) = 12sin(2"# 697 t) + 12sin(2"# 1209 t)Reference: Cleve Moler, Numerical Computing with MATLABfrequency (Hz)amplitude0.5time (seconds)soundpressure6Time Domain vs. Frequency DomainSignal. [recording, 8192 samples per second]Magnitude of discrete Fourier transform.Reference: Cleve Moler, Numerical Computing with MATLAB7Fast Fourier TransformFFT. Fast way to convert between time-domain and frequency-domain.Alternate viewpoint. Fast way to multiply and evaluate polynomials.“ If you speed up any nontrivial algorithm by a factor of a million or so the world will beat a path towards finding useful applications for it. ” — Numerical Recipeswe take this approach8Fast Fourier Transform: ApplicationsApplications.! Optics, acoustics, quantum physics, telecommunications, radar,control systems, signal processing, speech recognition, datacompression, image processing, seismology, mass spectrometry…! Digital media. [DVD, JPEG, MP3, H.264]! Medical diagnostics. [MRI, CT, PET scans, ultrasound]! Numerical solutions to Poisson's equation.! Shor's quantum factoring algorithm.! …“ The FFT is one of the truly great computational developments of [the 20th] century. It has changed the face of science and engineering so much that it is not an exaggeration to say that life as we know it would be very different without the FFT. ” — Charles van Loan9Fast Fourier Transform: Brief HistoryGauss (1805, 1866). Analyzed periodic motion of asteroid Ceres.Runge-König (1924). Laid theoretical groundwork.Danielson-Lanczos (1942). Efficient algorithm, x-ray crystallography.Cooley-Tukey (1965). Monitoring nuclear tests in Soviet Union andtracking submarines. Rediscovered and popularized FFT.Importance not fully realized until advent of digital computers.10Polynomials: Coefficient RepresentationPolynomial. [coefficient representation]Add. O(n) arithmetic operations.Evaluate. O(n) using Horner's method.Multiply (convolve). O(n2) using brute force. ! A(x) = a0+ a1x + a2x2+L + an"1xn"1 ! B(x) = b0+ b1x + b2x2+L + bn"1xn"1 ! A(x) + B(x) = (a0+ b0) + (a1+ b1)x + L + (an"1+ bn"1)xn"1 ! A(x) = a0 + (x (a1 + x (a2 + L + x (an"2 + x (an"1)) L ) )! A(x) " B(x) = cixii = 02n#2$, where ci= ajbi# jj = 0i$11A Modest PhD Dissertation Title“ New Proof of the Theorem That Every Algebraic Rational Integral Function In One Variable can be Resolved into Real Factors of the First or the Second Degree. ” — PhD dissertation, 1799 the University of Helmstedt12Polynomials: Point-Value RepresentationFundamental theorem of algebra. [Gauss, PhD thesis] A degree npolynomial with complex coefficients has exactly n complex roots.Corollary. A degree n-1 polynomial A(x) is uniquely specified by itsevaluation at n distinct values of x.xyxjyj = A(xj )13Polynomials: Point-Value RepresentationPolynomial. [point-value representation]Add. O(n) arithmetic operations.Multiply (convolve). O(n), but need 2n-1 points.Evaluate. O(n2) using Lagrange's formula. ! A(x) : (x0, y0), K, (xn-1, yn"1) B(x) : (x0, z0), K, (xn-1, zn"1) ! A(x) + B(x) : (x0, y0+ z0), K, (xn-1, yn"1+ zn"1)! A(x) = yk(x " xj)j#k$(xk" xj)j#k$k=0n"1% ! A(x) " B(x) : (x0, y0" z0), K, (x2n-1, y2n#1" z2n#1)14Converting Between Two RepresentationsTradeoff. Fast evaluation or fast multiplication. We want both!Goal. Efficient conversion between two representations ! all ops fast.coefficientrepresentationO(n2)multiplyO(n)evaluatepoint-valueO(n) O(n2)! a0, a1, ..., an -1 ! (x0, y0), K, (xn"1, yn"1)coefficient representationpoint-value representation15Converting Between Two Representations: Brute ForceCoefficient ! point-value. Given a polynomial a0 + a1 x + ... + an-1 xn-1,evaluate it at n distinct points x0 , ..., xn-1.Running time. O(n2) for matrix-vector multiply (or n Horner's). ! y0y1y2Myn"1 # $ % % % % % % & ' ( ( ( ( ( ( = 1 x0x02L x0n"11 x1x12L x1n"11 x2x22L x2n"1M M M O M1 xn"1xn"12L xn"1n"1 # $ % % % % % % & ' ( ( ( ( ( ( a0a1a2M an"1# $ % % % % % % & ' ( ( ( ( ( ( 16Converting Between Two Representations: Brute ForcePoint-value ! coefficient. Given n distinct points x0, ... , xn-1 and valuesy0, ... , yn-1, find unique polynomial a0 + a1x + ... + an-1 xn-1, that has givenvalues at given points.Running time. O(n3) for Gaussian elimination. ! y0y1y2Myn"1 # $ % % % % % % & ' ( ( ( ( ( ( = 1 x0x02L x0n"11 x1x12L x1n"11 x2x22L x2n"1M M M O M1 xn"1xn"12L xn"1n"1 # $ % % % % % % & ' ( ( ( ( ( ( a0a1a2M an"1# $ % % % % % % & ' ( ( ( ( ( ( Vandermonde matrix is invertible iff xi distinctor O(n2.376) via fast matrix multiplication17Divide-and-ConquerDecimation in frequency. Break up polynomial into low and high powers.! A(x) = a0 + a1x + a2x2 + a3x3 + a4x4 + a5x5 + a6x6 + a7x7.! Alow(x) = a0 + a1x + a2x2 + a3x3.! Ahigh (x) = a4 + a5x + a6x2 + a7x3.! A(x) = Alow(x) + x4 Ahigh(x).Decimation in time. Break polynomial up into even and odd powers.! A(x) = a0 + a1x + a2x2 + a3x3 + a4x4 + a5x5 + a6x6 + a7x7.! Aeven(x) = a0 + a2x + a4x2 + a6x3.! Aodd (x) = a1 + a3x + a5x2 + a7x3.! A(x) = Aeven(x2) + x Aodd(x2).18Coefficient to Point-Value Representation: IntuitionCoefficient ! point-value. Given a polynomial a0 + a1x + ... + an-1 xn-1,evaluate it at n distinct points x0 , ..., xn-1.Divide. Break polynomial up into even and odd powers.! A(x) = a0 + a1x + a2x2 + a3x3 + a4x4 + a5x5 + a6x6 + a7x7.! Aeven(x) = a0 + a2x + a4x2 + a6x3.! Aodd (x) = a1 + a3x + a5x2 + a7x3.! A(x) = Aeven(x2) + x Aodd(x2).! A(-x) = Aeven(x2) - x Aodd(x2).Intuition. Choose two points to be ±1.! A( 1) = Aeven(1) + 1 Aodd(1).! A(-1) = Aeven(1) - 1 Aodd(1).Can evaluate polynomial of degree " nat 2 points by evaluating two polynomialsof degree " !n at 1 point.we get to choose which ones!19Coefficient to Point-Value Representation: IntuitionCoefficient ! point-value. Given a polynomial a0 + a1x + ... + an-1 xn-1,evaluate it at n