Class 26. Fourier Transforms, Part 1Introduction• Cf. NRiC §12.0.• Fourier theorem: a well-behaved function can be represented by a series of sines andcosines of different frequencies and amplitudes.• Often useful to know what these frequencies and amplitudes are. Can do this with aFourier transf orm:H(f) =Z∞−∞h(t)e2πiftdt,where −∞ < f < ∞ is the frequency and H(f) is the amplitude (H is often complex,i.e., contains phase info).• Inverse Fourier transform:h(t) =Z∞−∞H(f)e−2πiftdf.• Units: if t is in seconds, f is in Hertz. If have h(x), x in m, then get H(n), n =wavenumber ( m−1).• FTs are linear ops:FT(g + h) = FT(g) + FT(h),FT(αh) = αFT(h).• h(t) may have special symmetries, e.g., pure real or pure imaginary, even (h(t) = h(−t))or odd (h(t) = −h(−t)) =⇒ can increase computational efficiency:h(t) pure real =⇒ H(−f) = H?(f)h(t) pure imaginary =⇒ H(−f ) = −H?(f)h(t) real & even =⇒ H(f) real & evenh(t) real & odd =⇒ H(f) imaginary & oddetc.Other properties, and combinations• If h(t) ⇐⇒ H(f) are a FT pair, thenh(at) ⇐⇒1|a|H(fa) “time scaling”1|b|h(tb) ⇐⇒ H(bf) “frequency scaling”h(t − t0) ⇐⇒ H(f)e2πift0“time shifting”h(t)e−2πif0t⇐⇒ H(f − f0) “frequency shifting”1• Combinations: if h(t) ⇐⇒ H(f ) and g(t) ⇐⇒ G(f), then1. Convolution:g ? h ≡Z∞−∞g(τ)h(t − τ) dτ.– Function of time. Note g ? h = h ? g.2. Convolution theorem:g ? h ⇐⇒ G(f)H(f)– E.g., instrumental pro file (point spread function): observe star, get PSF (con-volution of instrumental profile with delta function), now observe target, takeFT, divide by FT of PSF, take inverse FT to get deconvolved image.3. Correlation:corr(g, h) =Z∞−∞g(τ + t)h(τ) dτ.– Function of time, called “lag.”– Note corr(g, h) ⇐⇒ G(f)H(−f) = G(f)H?(f) if h(t) real.– Correlation used to compare dat a sets: it’s large at some t if functions areclose copies of each other but lead or lag in time by t. E.g., Doppler shift!4. Wiener-Khinchin theorem(autocorrelation) :corr(g, g) ⇐⇒ |G(f)|2.5. Parseval’s theorem:total power =Z∞−∞|h(t)|2dt =Z∞−∞|H(f)|2df.• Often interested in power between f and f + df. Usually regard f as varying f rom 0(D.C.) t o +∞ =⇒ one-sided power spectral density (PSD):Ph(f) ≡ |H(f )|2+ | H( −f)|2, 0 ≤ f < ∞.If h(t) real, Ph(f) = 2|H(f)|2.• If h(t) goes endlessly from −∞ < t < ∞, total power and PSD will generally b einfinite. Instead compute PSD per unit time, i.e. PSD/sample length. Area thencorresponds to mean square amplitude. As sample length → ∞, PSD per unit time →delta functions for pure sines and cosines.Discretely Sampled Data• Cf. NRiC §12.1.• For real data, often have hk≡ h(tk), tk= k∆, k = 0, 1, ..., N − 1. Here ∆ is thesampling interval; 1/∆ is the sampling rate.2• Define Nyquist critical frequency fc≡12∆. Critical sampling of a sine wave of frequencyfcis two points per cycle.• Sampling theorem: if signal is bandwidth limited such that H(f) = 0 for all |f| ≥ fc,then entire information content of signal can be recorded by sampling at ∆−1= 2fc.• If h(t) has power in frequencies outside −fc< f < fc, sampling h(t) causes power tospuriously move inside this range =⇒ aliasing:− fcfcH(f)H(f)ht f– Solution: filter signal and sample at least 2 points/cycle for highest fr equency.• If h(t) finite in time, N points should sample entire interval. If h(t) infinite, userepresentative portion.• N inputs =⇒ N outputs:fn≡nN∆, n = −N2, ...,N2.(For simplicity, assume N is even.) Extreme values of n ⇐⇒ Nyquist frequency range.• Now approximate:H(fn) =Z∞−∞h(t)e2πifntdt 'N−1Xk=0hke2πifntk∆ = ∆N−1Xk=0hke2πikn/N|{z }≡ Hn(DFT).• Note H−n= HN−nif n = 1, 2, ... (period N). Convention: let n = 0, 1, ..., N − 1 so nand k vary over same range. ∴ n = 0 ⇐⇒ zero frequency, n = N/2 ⇐⇒ f = fcandf = −fc. Hence:1 ≤ n ≤ N/2 − 1 ⇐⇒ 0 < f < fc,N/2 + 1 ≤ n ≤ N − 1 ⇐⇒ −fc< f < 0.Also note H(−f) ⇐⇒ Hn−N.• Discrete inverse Fourier transform:hk=1NN−1Xn=0Hne−2πikn/N.Very similar to Hn=⇒ can use same code...3Application: Solving Poisson’s Equation• Cf. NRiC §19.4.• Recall in 2-D the pro totypical elliptic equation is given by∂2u∂x2+∂2u∂y2= ρ(x, y).• The FD version is (assuming ∆x = ∆y ≡ ∆)uj−1,k− 2uj,k+ uj+1,k∆2+uj,k−1− 2uj,k+ uj,k+1∆2= ρj,k. (1)• Consider letting uj,kbe the 2-D inverse DFT of the Fourier-domain equivalent of u:uj,k=1JKJ−1Xm=0K−1Xn=0ˆum,ne−2πimj/Je−2πink/K. (2)(In multi-D, FTs can be computed independently in each dimension.)• Similarly,ρj,k=1JKJ−1Xm=0K−1Xn=0ˆρm,ne−2πimj/Je−2πink/K. (3)• Substituting (2) and (3) into (1), we getˆum,ne2πim/J+ e−2πim/J+ e2πin/K+ e−2πin/K− 4= ˆρm,n∆2,orˆum,n=ˆρm,n∆22cos2πmJ+ cos2πnK− 2. (4)• Strategy:1. Compute ˆρm,nas the FTˆρj,k=J−1Xj=0K−1Xk=0ρj,ke2πimj/Je2πink/K.2. Compute ˆum,nfrom (4).3. Compute uj,kby inverse FT (2).• Procedure valid only for periodic boundary conditions, i.e., for uj,k= uj+J,k= uj,k+K.• All we need now is a f ast way to compute the