LECTURE 15: FOURIER METHODS • We discussed different bases for regression in lecture 13: polynomial, rational, spline/gaussian… • One of the most important basis expansions is Fourier basis: Fourier (or spectral) transform • Several reasons for its importance: the basis is complete (any function can be Fourier expanded) • ability to compute convolutions and spectral densities (power spectrum) • Ability to convert dome partial differential equations (PDE) into ordinary differential equations (ODE) • Main reason for these advantages: fast Fourier transform (FFT)

Definition and properties of Fourier transforms

• Complete basise

Properties of Fourier transforms

• Convolution theorem

Correlation function and power spectrum

Power spectrum in higher dimensions

Discrete sampling: sampling theorem • We sample interval in points of length D: hn=h(nD) • Nyquist frequency: fc=1/2D • Sampling theorem: if the function h(t) does not have frequencies above fc (h(f)=0 for f>fc) it is bandwidth limited. Then h(t) is completely determined by hn: • This says that the information content is limited. If we know the maximum bandwidth frequency then we know how to sample the function using fc=1/2D to get the full information content

If we are not bandwidth limited we get aliasing

• All the power outside –fc
Discrete Fourier transform (DFT) • We measure a function on an interval ND. Maybe the function is 0 outside the interval, otherwise we assume periodicity (periodic boundary conditions), because sin and cos are periodic • Frequency range from –fc to fc • DFT: • H(-fc)=H(fc) • Inverse DFT • Parseval’s theorem

Fast Fourier Transform (FFT) • How expensive is to do FT? Naively it appears to be a matrix multiplication, hence O(N2) • FFT: O(Nlog2N) • The difference is enormous: these days we can have N>1010. FFT is one of the most important algorithms of numerical analysis • FFT existence became widely known in 1960s (Cooley & Tukey), but known (and forgotten) since Gauss (1805)

How FFT works? • Assume for now N=2M: one can show that FT of length N can be rewritten as the sum of 2 FTs of length N/2 (Danielson & Lanczos 1942): even (e) and odd (o)

• This can be recursively repeated, until we reach the length of 1, at which point we have for some n

Bit reversal • To determine n we take the sequence eoeeoo…oee, reverse the pattern of e and o, assign e=0 and o=1, and we get n in binary form.

Extensions • So far this was FFT of a complex function. FFT of real function can be done as FFT of complex of length N/2. • FFT of sine and cosine: it is common to use FFT to solve differential equations. The classic problem in physics are oscillations on a string. If the string ends are fixed the values are 0 there, hence we use sine. If they are open we often use derivative=0, for which we use cosine expansion.

FFT in higher dimensions • One can write it as a sequence of 1d FFTs

• 2d FFT of real function used in image processing: high pass, low pass filters, convolutions, deconvolutions… • 3d FFT: solving Poisson’s equation…

Image processing

Low pass filter: smoothly set high f/w to 0 High frequency sharpening: increase high f/w Derivative operator in FT: multiply Fourier modes with iw

(De)convolutions with FFT

• We have data that have been convolved with some response. For example, we observe the sky through the telescope which does not resolve angles below l/R, where R is its diameter. Here response is Airy function. To simulate the process we convolve data s(t) with r(t), by multiplying their FTs. • If the data are perfect (no noise) we can deconvolve the process: we FT the data r*s(t) to get r(f)s(f) and divide by the convolution term r(f) to get s(f), then inverse FT to get s(t).

Discrete version

• Discrete convolution assumes periodic signal on the intervale

Zero padding • Because discrete FFTs assume periodicity we need to zero pad the data. If r has non-zero support over M and is symmetric then zero pad M/2 pixels beyond the data.

Power spectrum and correlation function • Cross correlation function Corr(g,h)=FT-1[FT(g)FT(h)*] • Cross-power spectrum FT(g)FT(h)* • Auto-power spectrum (power spectrum density PSD, or periodogram) Pg(fk)=FT(g)FT(g)*

Window effects

• Because we observe on a finite interval the data are multiplied with the top-hat window, and the power spectrum is convolved with FT of top hat window • Top hat window has large side-lobes leaking from one f to another: mode mixing

• Top hat has sharp edges: if we smooth it then the side-lobes are reduced and the leakage is reduced (see project 2)

Apodization or windowing We multiply the signal with a window function, reducing sidelobes Note that we are reducing signal at the edges: suboptimal

w • w

w • w

w • w

w • w

Optimal (Wiener) filtering • Suppose we have noisy and convolved time stream of gaussian data and we wish to determine its best reconstruction in absence of noise and smearing • In absence of noise we would deconvolve, but in presence of noise this would amplify noise • This corresponds to a gaussian process: how do we choose the kernel KN? KN=Corr. • We also have kernel-basis function duality: a gaussian process corresponds to an infinite basis. Fourier basis is complete, so we can rewrite the kernel in terms of Fourier basis: Corr becomes power spectrum.

Wiener filter derivation • • • • •

Response convolution Noise Ansatz: determine F Least squares: minimize c2 No cross-term between N and S

• Take derivative wrt F: • Typically: F(f)=1 for low f, F(f)=0 for high f.

Wiener filter

Matched filtering • We have stationary but not diagonal noise and a template for signal, but we do not know the signal origin. • Example from project 2: LIGO chirp template on top of a correlated noise

Derivation of matched filter: 2d image • Stationary noise • Template t(x) • Linear ansatz • Bias • Variance • We minimize variance subject to zero bias: Lagrange multiplier

Matched filter

• Minimizing L wrt y

• MF solution filters out frequencies where signal to noise t/P is low • Signal to noise estimate • Procedure: we want to evaluate S/N for every point, so we want to convolve data with matched filter: we FT D(x) and template t, evaluate y(k), multiply with D(k), FT back, divide by s to get S/N, look for S/N peaks.

Project 2: LIGO matched filter analysis • W

Summary • Fourier methods revolutionized fields like image processing, Poisson solvers, PDEs, convolutions… • At the core is FFT algorithm scaling as Nlog2N • Many physical sciences use these techniques: see project 2, LIGO analysis of black hole black hole merger event

Literature • NR, Press etal. 12, 13

lecture 15: fourier methods - GitHub

LECTURE 15: FOURIER METHODS. â¢ We discussed different bases for regression in lecture. 13: polynomial, rational, spline/gaussianâ¦ â¢ One of the most important basis expansions is ... dome partial differential equations. (PDE) into ordinary differential equations (ODE). â¢ Main reason for these advantages: fast Fourier.

Recommend Documents

Lecture 1 - GitHub
Jan 9, 2018 - We will put special emphasis on learning to use certain tools common to companies which actually do data ... Class time will consist of a combination of lecture, discussion, questions and answers, and problem solving, .... After this da

Transcriptomics Lecture - GitHub
Jan 17, 2018 - Transcriptomics Lecture Overview. â¢ Overview of RNA-Seq. â¢ Transcript reconstrucån methods. â¢ Trinity de novo assembly. â¢ Transcriptome quality assessment. (coffee break). â¢ Expression quanæ¿aån. â¢ Differené¶¯l express

lecture 12: distributional approximations - GitHub
We have data X, parameters Î¸ and latent variables Z (which often are of the ... Suppose we want to determine MLE/MAP of p(X|Î¸) or p(Î¸|X) over q: .... We limit q(Î¸) to tractable distributions. â¢ Entropies are hard to compute except for tractable

lecture 4: linear algebra - GitHub
Inverse and determinant. â¢ AX=I and solve with LU (use inv in linalg). â¢ det A=L00. L11. L22 â¦ (note that Uii. =1) times number of row permutations. â¢ Better to compute ln detA=lnL00. +lnL11. +â¦

lecture 16: ordinary differential equations - GitHub
Bulirsch-Stoer method. â¢ Uses Richardson's extrapolation again (we used it for. Romberg integration): we estimate the error as a function of interval size h, then we try to extrapolate it to h=0. â¢ As in Romberg we need to have the error to be in

Old Dominion University Lecture 2 - GitHub
Old Dominion University. Department of ... Our Hello World! [user@host ~]\$ python .... maxnum = num print("The biggest number is: {}".format(maxnum)) ...

EE 396: Lecture 14-15
Calculating the posterior probability, p({Ri,ui}N i=1|I) using Bayes' Rule, and then calculating the MAP estimate is equivalent to minimizing the energy. E({Ri,ui}N.

lecture 5: matrix diagonalization, singular value ... - GitHub
We can decorrelate the variables using spectral principal axis rotation (diagonalization) Î±=XT. L. DX. L. â¢ One way to do so is to use eigenvectors, columns of X. L corresponding to the non-zero eigenvalues Î» i of precision matrix, to define new

C2M - Team 101 lecture handout.pdf - GitHub
Create good decision criteria in advance of having to make difficult decision with imperfect information. Talk to your. GSIs about their experience re: making ...

LECTURE 8: OPTIMIZATION in HIGHER DIMENSIONS - GitHub
We want to descend down a function J(a) (if minimizing) using iterative sequence of steps at a t . For this we need to choose a direction p t and move in that direction:J(a t. +Î·p t. ) â¢ A few options: fix Î·. â¢ line search: vary Î· ... loss (co

lecture 2: intro to statistics - GitHub
Continuous Variables. - Cumulative probability function. PDF has dimensions of x-1. Expectation value. Moments. Characteristic function generates moments: .... from realized sample, parameters are unknown and described probabilistically. Parameters a

lecture 10: advanced bayesian concepts - GitHub
some probability distribution function (PDF) of perfect data x, but what we measure is d, a noisy version of x, and noise is ... We can also try to marginalize over x analytically: convolve true PDF with noise PDF and do this for each ... We would li

GPU Multiple Sequence Alignment Fourier-Space Cross ... - GitHub
May 3, 2013 - consists of three FFTs and a sliding dot-product, both of these types of ... length n, h and g, and lets call this sum f. f indicates the similarity/correlation between ... transformed back out of Fourier space by way of an inverse-FFT.

Fast Fourier Transform Based Numerical Methods for ...
Analyzing contact stress is of a significant importance to the design of mechanical ...... A virtual ground rough surface can be formed through periodically extend-.

Brooklyn Community District 15 - GitHub
Public/Institutional. Open Space. Parking. Vacant. Other. 17,753. 2,161. 285. 1,150. 712. 35. 53. 245. 54. 194. 677. 87. BeltPkwy. N o stra n d. A v. Kings Hwy ... Highway, Manhattan Beach, Plumb Beach, Sheepshead Bay. Top 3 pressing issues identifie

Fourier series
Fourier series in an interval of length i2. Even Function. Odd Function. Convergence of Fourier Series: â¢ At a continuous point x = a, Fourier series converges to ...

lecture 6: information theory, entropy, experiment design - GitHub
LECTURE 6: INFORMATION THEORY,. ENTROPY, EXPERIMENT DESIGN. â¢ The concept of information theory and entropy appears in many statistical problems. â¢ Here we will develop some basic theory and show how it can be applied to questions such as how to

lecture 13: from interpolations to regressions to gaussian ... - GitHub
LECTURE 13: FROM. INTERPOLATIONS TO REGRESSIONS. TO GAUSSIAN PROCESSES. â¢ So far we were mostly doing linear or nonlinear regression of data points with simple small basis (for example, linear function y=ax+b). â¢ The basis can be arbitrarily larg

Lecture 11 â November 22, 2016 Volkswagen Emissions ... - GitHub
Emissions Workshop, an academic conference, in May 2014. Regulators and ... âThis VW Diesel Scandal is Much Worse Than a Re- call.â 21 September 2015.

LECTURE 7: NONLINEAR EQUATIONS and 1-d ... - GitHub
solving nonlinear equations and 1-d optimization here. â¢ f(x)=0 (either in 1-d or many ... dimensions. â¢ Requires a gradient: f(x+Î´)=f(x)+Î´f'(x)+â¦ â¢ We want f(x+Î´)=0, hence Î´=-f(x)/f'(x). â¢ Rate of convergence is quadratic (NR 9.4) Îµi+