Fourier Transforms and Spectral Analysis Chaiwoot Boonyasiriwat - - PowerPoint PPT Presentation

Chaiwoot Boonyasiriwat

September 17, 2020

Fourier Transforms and Spectral Analysis

Discrete-Time Fourier Transform

2

▪ The discrete-time Fourier transform (DTFT) of a discrete-time signal x(k) with sampling interval T is the Z transform X (z) evaluated along the unit circle on complex plane: ▪ “As a consequence, X ( f ) is well defined if and only if the region of convergence of X (z) includes the unit circle.”

Schilling and Harris (2012, p. 233)

Analysis Equation

Discrete-Time Fourier Transform

3

▪ “For a causal signal x(k), the DTFT X ( f ) exists if x(k) is absolutely summable or if the poles of X (z) all lie strictly inside the unit circle.” ▪ The DTFT X ( f ) is called the spectrum of x(k) and can be written in the polar form as where is the magnitude spectrum and is the phase spectrum of x(k).

Schilling and Harris (2012, p. 234)

Analysis and Synthesis Equations

4

▪ “The DTFT is called the analysis equation because it decomposes a signal into its spectral components.” ▪ “The signal x(k) can be recovered from its spectrum X ( f ) using the inverse transform” ▪ “The IDTFT is called the synthesis equation because it reconstructs or synthesizes the signal from its spectral components.”

See the derivation in Schilling and Harris (2012, p. 234, equation 4.2.3)

Properties of DTFT

5

Schilling and Harris (2012, p. 235)

▪ Since X ( f ) is periodic with period fs, the frequency is typically restricted to the interval [-fs/2, fs/2]. ▪ “For real signals, information about the spectrum is contained in the frequency range [0, fs/2] due to the symmetry of X ( f ).

Spectrum of Causal Exponential

6

▪ Consider the causal exponential ▪ The Z transform of this signal is ▪ Since the pole of X (z) is at z = c, DTFT converges only when |c| < 1. ▪ The spectrum of x(k) is ▪ The magnitude and phase spectra are

Schilling and Harris (2012, p. 235)

Spectrum of Causal Exponential

7

Schilling and Harris (2012, p. 236)

Normalized frequency

Frequency Ranges of Some Signals

8

Schilling and Harris (2012, p. 236)

9

Schilling and Harris (2012, p. 239)

2

( ) X f L =

IDTFT of Ideal Lowpass Filter

10

“An ideal lowpass filter with a cutoff frequency of 0 < Fc < fs/2 has a phase response of ( f ) = 0 and a magnitude response consisting of a rectangular window of radius Fc.” The impulse response hlow(k) can be computed by

Schilling and Harris (2012, p. 240)

IDTFT of Ideal Lowpass Filter

11

Schilling and Harris (2012, p. 240)

Basic DTFT Pairs

12

Schilling and Harris (2012, p. 241)

The DTFT of causal signal x(k) is This leads to two drawbacks: ▪ An infinite number of floating point operations (FLOPs) is required to evaluate X ( f ). ▪ The transform must be evaluated at an infinite number

• f frequencies f.

These drawbacks can be overcome by ▪ Focusing on finite signals ▪ Evaluating the transform at N equally spaced values of frequency f

Drawbacks of DTFT

13

Schilling and Harris (2012, p. 241)

Discrete Fourier Transform (DFT)

14

If x(k) is absolutely summable, then For a sufficiently large value of N, X( f ) can be approximated by the finite sum Then X ( f ) is evaluated at N discrete frequencies equally spaced over one period of X ( f ): Let zi be the point in the complex plane corresponding to frequency fi: . Since |zi| = 1 , the N evaluation points are equally around the unit circle.

Schilling and Harris (2012, p. 241)

Evaluation Points of DFT

15

Schilling and Harris (2012, p. 242)

N = 8

DFT of causal N-point signal x(k) is then defined as IDFT:

DFT and IDFT

16

Schilling and Harris (2012, p. 242-243)

DFT is a mapping from an input vector x to an output vector X: DFT is a linear transformation and can be represented by an NN matrix whose elements are . Example: N = 5,

Matrix Formulation of DFT

17

Schilling and Harris (2012, p. 243)

DFT can be represented in the matrix form as Multiplying both sides by W-1 yields Comparing the equations defining DFT and IDFT, we know that So the matrix form of IDFT is

Matrix Formulation of DFT and IDFT

18

Schilling and Harris (2012, p. 243-244)

Suppose the input are Thus, N = 4 and DFT of x is

Matrix Formulation of DFT: Example

19

Schilling and Harris (2012, p. 244)

From the previous example, we have IDFT of X is

Matrix Formulation of IDFT: Example

20

Schilling and Harris (2012, p. 245)

▪ “Periodic signals have a discrete spectrum.” ▪ “Suppose xa(t) is a periodic continuous-time signal with period T0 and fundamental frequency F0 = 1/T0.” ▪ Then, xa(t) can be expanded into a complex Fourier series as ▪ The Fourier coefficient is ▪ Let x(k) = xa(kT) be the kth sample of using the sampling interval of T. ▪ We want to compute the spectrum of x(k).

Fourier Series and Discrete Spectra

21

Schilling and Harris (2012, p. 245)

The IDTFT of X ( f ) = a(f – F0) is

• r

So, the spectrum of the periodic signal x(k) is which is zero everywhere except at the harmonic frequencies iF0 – discrete-frequency spectrum.

Fourier Series and Discrete Spectra

22

Schilling and Harris (2012, p. 245)

Periodic continuous-time signal xa(t) can be approximated by truncating the Fourier series to M harmonics: Suppose that xa(t) is sampled at N = 2M points using a sampling rate of fs = NF0. So, T0 = 1/F0 = N/fs = NT. That is the N samples cover one period of xa(t).

Fourier Coefficients

23

Schilling and Harris (2012, p. 246)

The Fourier coefficient can be approximated as and can be obtained from the DFT of the samples of xa(t).

Fourier Coefficients

24

Schilling and Harris (2012, p. 246)

▪ “For a real signal xa(t), c-i = ci

*. Thus, the complete set

• f Fourier coefficients is

▪ The Fourier series can also be expressed as The coefficients di and phase angle i can be obtained using the DFT by

Fourier Coefficients

25

Schilling and Harris (2012, p. 246)

DFT Properties

26

Schilling and Harris (2012, p. 251, 255)

Parseval’s Identity

27

Schilling and Harris (2012, p. 254)

▪ Let x and y be two N-point signals with DFT X and Y. ▪ It can be shown that ▪ When x = y, we then have ▪ The power density spectrum is ▪ The average power of N-point signal is

▪ The DFT of N-point signal x(k) is ▪ “The N2 values of do not depend on x(k), so they can be precomputed and stored in an NN matrix.” ▪ “Each point X(i) requires N complex multiplications.” ▪ “So, the total number of complex FLOPs required to compute the entire DFT is .” ▪ The computational cost of DFT is of order O(N2). ▪ “A computational algorithm is of order O(Np) if and

• nly if the number of FLOPs n satisfies

Fast Fourier Transform (FFT)

28

Schilling and Harris (2012, p. 256)

▪ Suppose the number of data points is N = 2r for some integer r. ▪ Decimate the N-point signal x(k) into two N/2-point signals xe(k) and xo(k) corresponding to the even and

• dd indices of x, respectively:

▪ “The DFT of x(k) then can be recast as two sums, one corresponding to the even values of k, and the other corresponding to the odd values of k.”

Decimation in Time FFT

29

Schilling and Harris (2012, p. 256)

The DFT of x(k) then becomes where Xe(i) = DFT{xe(k)} and Xo(i) = DFT{xo(k)} are N/2- point transforms of the even and odd parts of x(k).

Decimation in Time FFT

30

Schilling and Harris (2012, p. 257)

▪ The total FLOPs for the even-odd decomposition is ▪ Breaking the merging formula into two cases where yields ▪ The periodic property of the N/2-point transforms gives ▪ WN also has the symmetry property

Decimation in Time FFT

31

Schilling and Harris (2012, p. 257)

We then have, for , This computation is called the ith-order butterfly which reduces the number of multiplication from 2 to 1 but increases a storage by one scalar.

Decimation in Time FFT

32

Schilling and Harris (2012, p. 257)

Signal flow graph of the ith-order butterfly computation

FFT when N = 8

33

Schilling and Harris (2012, p. 258)

▪ The even-odd decomposition can be repeatedly applied. ▪ The even samples xe(k) and the odd samples xo(k) can each be further decimated so the N/2-point transform is turned into a pair of N/4-point transforms. ▪ Since N = 2r, this process can be continued for r times. ▪ “In the end, we are left with a collection of elementary 2-point DFTs which are the zeroth-order butterfly.” ▪ “The algorithms is called the radix-two fast Fourier transform or FFT.”

FFT

34

Schilling and Harris (2012, p. 258)

Signal Flow Graph of FFT: N = 8

35

Schilling and Harris (2012, p. 259)

The number of FLOPs for the FFT of N-point signal is So, the FFT algorithm is of order O(N log2 N).

Computational Cost of FFT

36

Schilling and Harris (2012, p. 260)

The IDFT formula is and the complex conjugate of WN is We then have

IDFT using FFT

37

Schilling and Harris (2012, p. 261-262)

Given causal signal x(k) of length M and h(k) of length N. The linear convolution of h(k) with x(k) is which can be computed using circular convolution of zero-padded signals hz(k) and xz(k) of length N = L+M+p where : Using the DFT property, we have So, an efficient way to perform circular convolution is

Fast Convolution

38

Schilling and Harris (2012, p. 263-264)

▪ To perform a fast linear convolution using FFT, the length of zero-padded signals must be a power of 2: where nextpow2 is a MATLAB function that finds the smallest power of two that is greater than or equal to its argument. ▪ The linear convolution can then be performed using the fast convolution

Fast Convolution

39

Schilling and Harris (2012, p. 264)

Block Diagram of Fast Convolution

40

Schilling and Harris (2012, p. 264)

Given two signals of length N (power of two). Cost of fast convolution is Cost of direct convolution is

Fast Convolution Cost

41

Schilling and Harris (2012, p. 265)

The linear cross-correlation of L-point signal y(k) with M-point signal x(k), M  L, is given by where is the circular cross-correlation of N-point zero-padded signal yz(k) by N-point zero-padded signal xz(k), N = L + M + p, p  -1.

Fast Correlation

42

Schilling and Harris (2012, p. 270)

Using the DFT property the circular cross-correlation can be performed by To use FFT, N = L + M + p must be a power of two: The linear cross-correlation can then be efficiently performed using FFT by The cost of fast correlation is also O(N log 2N).

Fast Correlation

43

Schilling and Harris (2012, p. 270-271)

▪ The probability that a random variable x lies in an interval [c, d] is ▪ The expected value of f (x) is denoted as E[ f (x)] and is defined as ▪ The kth moment of x is defined as for k  1. ▪ The mean of x is defined as ▪ The kth central moment is defined as ▪ Variance of x (second central moment)

▪ Standard deviation  is the square root of variance.

Review of Statistics

44

Schilling and Harris (2012, p. 275-276)

▪ If x and y are independent random variables, then

Review of Statistics

45

Schilling and Harris (2012, p. 275-276)

▪ White noise is a type of random signal useful for modeling physical signals corrupted with noise. ▪ White noise contains power at all frequencies and, therefore, is suitable as a test signal for system identification. ▪ Here, we are introducing the two important types of random signals: uniform white noise and Gaussian white noise.

White Noise

46

Schilling and Harris (2012, p. 274)

▪ Uniform white noise is a random signal characterized by a random variable x with the uniform probability density function ▪ White noise contains power at all frequencies similar to white light that contains all colors. ▪ “To create an array v of N random numbers uniformly distributed over [a, b], the following code can be used:” v = a + (b-a)*rand(N,1);

Uniform White Noise

47

Schilling and Harris (2012, p. 274, 276)

▪ The average power of x is the second moment E[x2]. ▪ The average power of the uniformly distributed random variable x is ▪ When [a, b] = [-c, c], Px = c2/3. ▪ In summary, the average power of uniform white noise is a constant.

Uniform White Noise

48

Schilling and Harris (2012, p. 276-277)

▪ Gaussian white noise is modeled using a normal or Gaussian probability density function where  is the mean and  the standard deviation. ▪ To create an array v of N random numbers with a Gaussian distribution of mean mu and standard deviation sigma, the following MATLAB code can be used: v = mu + sigma*randn(N,1); ▪ The average power of a Gaussian random variable x is

Gaussian White Noise

49

Schilling and Harris (2012, p. 278-280)

▪ Consider the continuous-time periodic signal produced by an AM mixer: ▪ “Suppose F1 = 300 Hz, F2 = 100 Hz, and xa(t) is corrupted with Gaussian noise v(k) with zero mean and standard deviation  = 0.8. The result is then sampled at fs = 1 kHz using N = 1024 samples, the discrete-time signal is

Gaussian White Noise: Example

50

Schilling and Harris (2012, p. 278-280) Power Density Spectrum F1-F2 F1+F2

   

1 2 1 2

2 ( ) sin 2 sin 2 ( ( ) )

a

x F F t t t F F   + + = −

▪ Cross-correlation of a signal with itself is called auto- correlation. ▪ Let x(k) be an N-point signal and xp(k) be its periodic

• extension. The circular auto-correlation of x(k) is

defined as ▪ The zero-lag auto-correlation is the average power of x(k) since xp(0) = x(0). That is ▪ The normalized auto-correlation is defined as

Auto-correlation

51

Schilling and Harris (2012, p. 282)

▪ Let v(k) be a zero-mean white noise signal of length N. ▪ If a random signal is ergodic, the expected value E[ f (v(k))] can be approximated by replacing the ensemble average with the time average: ▪ The circular auto-correlation of v(k) can be expressed in terms of expected value as

Auto-correlation of White Noise

52

Schilling and Harris (2012, p. 282)

▪ “The samples of a white noise signal are statistically independent of one another.” ▪ Since E[v(i)] = 0, it follows that ▪ Uncorrelated signals are statistically independent and at least one signal has a zero mean so the expected value of the product is zero. ▪ Since cvv(0) = Pv and cvv(k)  0 when k  0, we then have ▪ Applying the same analysis to linear auto-correlation yields

Auto-correlation of White Noise

53

Schilling and Harris (2012, p. 282)

Normalized linear auto-correlation of a N-point white noise signal with  = 0,  = 1, N = 1024.

Auto-correlation of White Noise

54

Schilling and Harris (2012, p. 284)

▪ The power density spectrum of an N-point signal x(k) is ▪ Applying DFT to cxx(k) yields ▪ “Thus, the DFT of the circular auto-correlation of x(k) is the power density spectrum:”

Power Density Spectrum

55

Schilling and Harris (2012, p. 284)

Compute and plot the power density spectrum of the double pulse signal with N = 512, M = 8:

Exercise

56

Schilling and Harris (2012, p. 285)

▪ The frequency response H( f ) of a stable linear system is the transfer function H(z) evaluated along the unit circle. ▪ H( f ) is also the spectrum of the impulse response. ▪ For a causal system, ▪ H( f ) can also be approximated by DFT. ▪ Let H(i) denote the DFT of the first N samples of h(k): ▪ Suppose we evaluate the frequency response at discrete frequencies fi = ifs/N.

Frequency Response using DFT

57

Schilling and Harris (2012, p. 291)

Recall that . We then have

Frequency Response using DFT

58

Schilling and Harris (2012, p. 291)

“For a stable filter, the impulse response is absolutely summable, so the RHS is finite and goes to 0 as N → .

Frequency Response using DFT

59

Schilling and Harris (2012, p. 291)

▪ So, for a sufficiently large N and a real system, the discrete-time frequency response can be approximated by the DFT of the impulse response: ▪ In this case, H(i) exhibits symmetry about the middle point i = N/2 and all essential information is contained in the subinterval 0  i  N/2, which corresponds to positive frequencies. ▪ This approximation is exact for some filters. ▪ Example: FIR filter of order m: h(k) = 0 for k > m. So, the upper bound on the error is zero when N > m.

Frequency Response using DFT

60

Schilling and Harris (2012, p. 292)

Consider a moving average filter with impulse response

Frequency Response: Example

61

Schilling and Harris (2012, p. 292)

Magnitude response: Phase response: Frequencies in plots are

▪ Decibel (dB) is a logarithmic unit often used for plotting the magnitude response of a filter or the magnitude spectrum of a signal. ▪ A decibel of signal x(k) is defined as 10 log10 |x(k)|2. ▪ The magnitude response in decibel is then defined as

Decibel Scale (dB)

62

Schilling and Harris (2012, p. 293-294)

▪ “The length of the signal can be increased from N to M by padding x(k) with M – N zeros as follows.” ▪ Suppose Sz(i) is the power density spectrum of xz(k) scaled by M/N: where Xz(i) = DFT{xz(k)}.” ▪ The frequency precision (space between adjacent frequencies) of Sz(i) is

Zero Padding

63

Schilling and Harris (2012, p. 295)

▪ “The effect of zero padding is to interpolate between the

• riginal points of the spectrum.”

▪ “If M = qN for some integer q, there will be q – 1 new points between each of the original N points. The

• riginal points remain unchanged.”

▪ “By decreasing  fz through zero padding, the location

• f an individual sinusoidal spectral component can be

more precisely determined as shown in the following example.”

Zero Padding

64

Schilling and Harris (2012, p. 295)

▪ Suppose fs = 1024 and N = 256. ▪ Consider a noise-free signal with a single sinusoidal spectral component at F0 = 330.5 Hz. ▪ Zero-padding x(k) by a factor of 8, then M = 8N = 2048. ▪ The frequency precision before and after zero padding is ▪ The corresponding sinusoidal frequency indices for the two cases are

Zero Padding: Example

65

Schilling and Harris (2012, p. 295)

The peaks in the power density spectra will occur at

Zero Padding: Example

66

Schilling and Harris (2012, p. 295)

▪ Zero padding can help identify the precise location of an isolated sinusoidal frequency component. ▪ However, zero padding cannot help resolve the difference between two closely spaced components. ▪ The frequency resolution or the smallest frequency difference that can be detected is limited by the sampling frequency and the number of nonzero samples. ▪ The Rayleigh limit represents the value of the frequency resolution (reciprocal of the length of the

• riginal signal):

Spectral Resolution

67

Schilling and Harris (2012, p. 296)

▪ Frequency precision is the spacing between discrete frequencies of the DFT. ▪ Frequency resolution is the smallest frequency difference that can be detected. ▪ To increase the frequency resolution, we need to add more samples to the signal. Zero padding does not help.

Spectral Resolution

68

Schilling and Harris (2012, p. 296)

▪ Suppose fs = 1024 Hz and N = 1024. ▪ Consider the noise-free signal with spectral components at F0 = 330 Hz and F1 = 331 Hz: ▪ Zero padding x(k) by a factor a two (M = 2N = 2048). The frequency precision before and after zero padding is ▪ The Rayleigh limit is ▪ To decrease the Rayleigh limit, we need to add more nonzero samples to the signal, e.g., doubling the signal length to N = 2048 or decreasing the sampling rate.

Spectral Resolution: Example

69

Schilling and Harris (2012, p. 297)

Spectral Resolution: Example

70

Schilling and Harris (2012, p. 297)

N = 1024

Spectral Resolution: Example

71

Schilling and Harris (2012, p. 298)

N = 1024 M = 2048

Spectral Resolution: Example

72

Schilling and Harris (2012, p. 297)

N = 2048

▪ Spectral characteristics of some signals can change with time, e.g., speech signals. ▪ To capture the time-varying spectral characteristic, we need to compute the DFT of overlapping segments of the signal. ▪ This is called short-time Fourier transform (STFT).

Short-Time Fourier Transform

73

Schilling and Harris (2012, p. 299)

▪ Suppose x(k) is an N-point signal with N = LM for integers L and M. ▪ Then x(k) can be decomposed into 2M – 1 overlapping subsignals of length L. ▪ The subsignals can be extracted from the signal x(k) as

Extracting Subsignals

74

Schilling and Harris (2012, p. 301)

M = 5

▪ The subsignals can also be represented as the product between original signal x(k) and a window function wR: ▪ Various window functions can be used to extract a subsignal from the original signal:

• Rectangular:
• Hanning:
• Hamming:
• Blackman

Window Functions

75

Schilling and Harris (2012, p. 301)

Window Functions

76

Schilling and Harris (2012, p. 301)

▪ Let x(k) be an N-point signal decomposed into 2M – 1

• verlapping subsignals xm(k) of length L.

▪ Let w(k) be a window function. ▪ The spectrogram of x(k) is a (2M – 1)L matrix G defined as

Spectrogram

77

Schilling and Harris (2012, p. 301)

Consider the signal x(k) composed of a recording of the vowels {A, E, I, O, U} with the sampling rate fs = 8 kHz and N = 32,000. Choosing L = 256, then M = 125.

Spectrogram: Example

78

Schilling and Harris (2012, p. 302)

Using the Hamming window yields the spectrogram:

Spectrogram: Example

79

Schilling and Harris (2012, p. 303)

▪ Schilling, R. J. and S. L. Harris, 2012, Fundamentals

• f Digital Signal Processing using MATLAB, Second

Edition, Cengage Learning.

Reference