ESCI 386 Scientific Programming, Analysis And .

ESCI 386 – Scientific Programming,Analysis and Visualization withPythonLesson 17 - Fourier Transforms1

Spectral Analysis Most any signal can be decomposed into asum of sine and cosine waves of variousamplitudes and wavelengths.2

Fourier Coefficients For each frequency of wave contained in thesignal there is a complex-valued Fouriercoefficient. The real part of the coefficient containsinformation about the amplitude of the cosinewaves The imaginary part of the coefficient containsinformation about the amplitude of the sinewaves3

Discrete Fourier Transform (DFT) The discrete Fourier transform pair N is number of data points1 N 1F (m) f (n) exp i 2 mn N Forward transformN n 0N 1f (n) F (m) exp i 2 mn N m 0Backward orinverse transform4

Discrete Fourier Transforms f (n) is value of function f at grid point n. F(m) is the spectral coefficient for the mthwave component.1 N 1F (m) f (n) exp i 2 mn N N n 0N 1f (n) F (m) exp i 2 mn N m 05

Discrete Fourier Transforms f (n) may be real or complex-valued. F(m) is always complex-valued.1 N 1F (m) f (n) exp i 2 mn N N n 0N 1f (n) F (m) exp i 2 mn N m 06

DFT There will be as many Fourier coefficients (waves)as there are grid points. The natural frequencies/wave numbersassociated with the Fourier coefficients arem m ; m 0,1, 2, sN mN m 1 ; sN, N 2for positive frequenciesm N 2 1, N 2 2, N 2 3,,N N 2for negative frequencies7

The Frequencies/Wave Numbers If the input function f(n) is entirely real thenthe negative frequencies are completelyredundant, giving no additional information. The highest natural frequency (wave number)that can be represented in a digital signal iscalled the Nyquist frequency, and is given as1/2 where is the sampling interval/period.8

The zeroth Coefficient The zeroth Fourier coefficient is the averagevalue of the signal over the interval.9

Show Visualizations10

Power Spectrum The Fourier coefficients, F(m), are complexnumbers, containing a real part and animaginary part. The real part corresponds to the cosine wavesthat make up the function (the even part ofthe original function), and the negative of theimaginary terms correspond to the sine waves(the odd part of the original function).11

Power Spectrum For a sinusoid that has an integer number ofcycles within the N data points, the amplitude ofthe sine or cosine wave is twice its Fouriercoefficient.– So, a pure cosine wave of amplitude one would have asingle real Fourier coefficient at its frequency, and thevalue of that coefficient would be 0.5. The magnitude squared of the Fouriercoefficients , F(m) 2 , is called the power.– A plot of the power vs. frequency is referred to as thepower spectrum.12

Power Spectra for Square Signal13

Power Spectra for Gaussian Signal14

The numpy.fft ModuleFunctionPurposefft(s)Computes the forward DFT and The returned array is a complex array.returns the coefficients Fifft(F)Computes the inverse DFT andreturns the signal sfftfreq(n,d) Computes the naturalfrequencies/wavenumbers. d isan optional sample spacing(default is 1).fftshift(F)RemarksThe zero frequency is in the firstposition of the array, followed by thepositive frequencies in ascendingorder, and then the negativefrequencies in descending orderShifts the zero frequency to the This can be used on either thecenter of the array.frequencies or the spectral coefficientsto put the zero frequency in the center.15

The numpy.fft Module (cont.) There are also functions for taking FFTs in twoor more dimensions, and for taking FFTs ofpurely real signals and returning only thepositive coefficients. See documentation fordetails. We normally import this asfrom numpy import fft16

The numpy.fft.fft() Function The fft.fft() function accepts either a realor a complex array as an input argument, andreturns a complex array of the same size thatcontains the Fourier coefficients. For the returned complex array:– The real part contains the coefficients for the cosineterms.– The imaginary part contains the negative of thecoefficients for the sine terms.17

The fft.fft() Function (cont.) IMPORTANT NOTICE!!! There are differentdefinitions for how a DFT should be taken.18

The fft.fft() Function (cont.) NumPy actually uses these equations! This means that NumPy’s Fourier coefficientswill be N times larger than expected!N 1F (m) f (n) exp i 2 mn N n 01f ( n) NN 1 F (m) exp i2 mn N m 019

The fft.fftfreq() Function The natural frequencies associated with thespectral coefficients are calculated using thefft.fft() function. The zeroth frequency is first, followed by thepositive frequencies in ascending order, andthen the negative frequencies in descendingorder.20

FFT Example Programfrom numpy import fftimport numpy as npimport matplotlib.pyplot as pltn 1000# Number of data pointsdx 5.0# Sampling period (in meters)x dx*np.arange(0,n)# x coordinatesw1 100.0 # wavelength (meters)w2 20.0# wavelength (meters)fx np.sin(2*np.pi*x/w1) 2*np.cos(2*np.pi*x/w2) # signalFk fft.fft(fx)/n# Fourier coefficients (divided by n)nu fft.fftfreq(n,dx) # Natural frequenciesFk fft.fftshift(Fk)# Shift zero freq to centernu fft.fftshift(nu)# Shift zero freq to centerf, ax plt.subplots(3,1,sharex True)ax[0].plot(nu, np.real(Fk))# Plot Cosine termsax[0].set ylabel(r' Re[F k] ', size 'x-large')ax[1].plot(nu, np.imag(Fk))# Plot Sine termsax[1].set ylabel(r' Im[F k] ', size 'x-large')ax[2].plot(nu, np.absolute(Fk)**2) # Plot spectral powerax[2].set ylabel(r' \vert F k \vert 2 ', size 'x-large')ax[2].set xlabel(r' \widetilde{\nu} ', size 'x-large')plt.show()File: fft-example.py21

FFT Result22

FFT Leakage There are no limits on the number of data points whentaking FFTs in NumPy. The FFT algorithm is much more efficient if the number ofdata points is a power of 2 (128, 512, 1024, etc.). The DFT assumes that the signal is periodic on the interval0 to N, where N is the total number of data points in thesignal. Care needs to be taken to ensure that all waves in thesignal are periodic within the interval 0 to N, or aphenomenon known as leakage will occur.23

Leakage Example In prior example if wavelength of the sine wave ischanged to 110 meters from 100 meters.Leakage24

The fft.rfft() Function Since for real input the negative frequenciesare redundant, there is an fft.rfft()function that only computes the positivecoefficients for a real input. The helper routine fftfreq() only returnsthe frequencies for the complex FFT.25

The fft.rfft() Function To return the frequencies associated withfft.rfft() function for an input of n datapoints, use:CaseUseFor n evenf fft.fftfreq(n)freq np.abs(f[0:n/2 1])For n oddf fft.fftfreq(n)freq f[0:n/2 1]26

The numpy.fft.fft() Function The fft.fft() function accepts either a real or a complex array as an input ar