No Spectral Leakage Time Domain Frequency Domain Signal

  • Slides: 55
Download presentation
No Spectral Leakage Time Domain Frequency Domain • Signal: sin(2π * 3 fs) •

No Spectral Leakage Time Domain Frequency Domain • Signal: sin(2π * 3 fs) • Three periods per sampling time period – When f=3, DFT correlation squares each value; result >0 – Result is N/2 times the time domain amplitude • The four periods per time period – Positives and negatives cancel each other – FFT spectral value for f=4 totals zero

Spectral Leakage Time Domain Frequency Domain • Signal: sin(2π * 3. 5 fs) •

Spectral Leakage Time Domain Frequency Domain • Signal: sin(2π * 3. 5 fs) • Three periods per sampling time period – When f=3, DFT correlation attenuates the total – Some of the amplitude spills across the spectrum

DFT Window-sync Distribution • Window sync function: sin(x)/x – Period is the same as

DFT Window-sync Distribution • Window sync function: sin(x)/x – Period is the same as sin(x) – Amplitude degrades as x increases – Sin(π(k-m))/{π(k-m)}; k is signal frequency, m is DFT basis function • Integer k points compute to zero because sin(πk) = 0 • Non integer k points ≠ 0 because sin(π(k-m)) is not a multiple of π

Spectral Leakage Examples

Spectral Leakage Examples

Leakage Symmetry 3. 4 cycles period 64 point DFT Circular view of the frequency

Leakage Symmetry 3. 4 cycles period 64 point DFT Circular view of the frequency spectrum

Leakage Symmetry (cont. ) Linear view of the frequency spectrum

Leakage Symmetry (cont. ) Linear view of the frequency spectrum

More Examples

More Examples

Scalloping Loss • • • Leakage means power loss in non-center frequencies Loss amount

Scalloping Loss • • • Leakage means power loss in non-center frequencies Loss amount greatest for frequencies at the intersections Loss can be significant (0. 637 in this example) Result is the scalloped pattern of frequency amplitudes Solutions: – Zero DFT creates finer frequency resolution, and therefor causes the scalloping low points to be closer together – Applying a windowing function alters the scalloping curve shape

Windowing • Issue: In practice, we apply DFT to finite signal pieces (frames) •

Windowing • Issue: In practice, we apply DFT to finite signal pieces (frames) • Impact: – The start of the truncated signal does not typically line up with the end (rough edge) – The frequency domain experiences distortion • Using a window – First Create the window: Create an array of frame length with values calculated with a windowing formula – Apply the window: multiply these values by each time domain frame amplitude

Speech Application Speech signal characteristics change as we speak • Goal: Extract short-term signal

Speech Application Speech signal characteristics change as we speak • Goal: Extract short-term signal features from speech signal frames • Approach – – – Block the speech into short overlapping segments Duration of each block varies from 10 to 30 ms Assumption: Speech within a single frame is unchanging Apply a window function to smooth the signal Apply DFT after applying the window Perform frequency analysis to extract speech characteristics

Window Function Application of windows to sinusoid signal 1. rectangular (no) window 2. Trangular

Window Function Application of windows to sinusoid signal 1. rectangular (no) window 2. Trangular window 3. Hanning window

Popular Window Types • Rectangular: wk = 1 where k = 0 … M

Popular Window Types • Rectangular: wk = 1 where k = 0 … M – Advantage: Easy to calculate, array elements unchanged – Disadvantage: Messes up the frequency domain • Ideal (window-sync): sin( 2 π f i) / (πi) – Disadvantage: Must be infinitely long – Truncation causes ripples and overshoots • Hamming: wk = 0. 54 – 0. 46 cos(2 kπ/M) – Advantage: Fast roll-off in frequency domain – Disadvantage: worse attenuation • Blackman: wk = 0. 42 – 0. 5 cos(2 kπ/M) + 0. 08 cos(4 kπ/M) – Advantage: better attenuation – Disadvantage: slower roll-off Multiply the window, point by point, to the audio signal

Other Windowing Formulae • Hanning: w[n] = 0. 5 -0. 5 cos(2πn/(N-1)) – Advantage:

Other Windowing Formulae • Hanning: w[n] = 0. 5 -0. 5 cos(2πn/(N-1)) – Advantage: Adding 50% overlapped frames produces no gain – Advantage: Excellent for pitch manipulation algorithms • Bartlett: w[n] = 2/(N-1)/2 - |(n–(N-1)/2|) • Triangle: w[n] = 2/N. (N/2 - |(n–(N-1)/2|) • Blackman: w[n] = a 0 – a 1 cos(2πn/(N-1)) – a 2 cos(4πn/(N-1)) where: a 0 = (1 -α)/2 ; a 1 = ½ ; a 2 = α /2; α =0. 16 • Blackman-Harris: w[n] = a 0–a 1 cos(2πn/(N-1))–a 2 cos(4πn/(N-1)) - a 3 cos(6πn/(N-1)) where: a 0=0. 35875; a 1=0. 48829; a 2=0. 14128; a 3=0. 01168 Note: Small changes can cause large frequency domain impact

Code: Hamming Window Create Window double[] window = new double[window. Size]; double c =

Code: Hamming Window Create Window double[] window = new double[window. Size]; double c = 2*Math. PI / (window. Size - 1); for (int h=0; h<window. Size; h++) window[h] = 0. 54 - 0. 46*Math. cos(c*h); Apply window for(int i=0; i<window. length; i++) frame[i]=frame[i]*window[i];

Rectangular Window Frequency Response Time Domain Filter Top Right: Amplitude linear scale, Bottom Right:

Rectangular Window Frequency Response Time Domain Filter Top Right: Amplitude linear scale, Bottom Right: Power logarithmic scale

Understanding Digital Signal Processing, Third Edition, Richard Lyons (0 -13 -261480 -4) © Pearson

Understanding Digital Signal Processing, Third Edition, Richard Lyons (0 -13 -261480 -4) © Pearson Education, 2011.

Windowing Frequency Response • Main Lobe: Narrow implies better frequency resolution. As the window

Windowing Frequency Response • Main Lobe: Narrow implies better frequency resolution. As the window length grows, the main lobe width narrows (less initial spectral leakage). • Side lobe: Higher for abrupt time domain window transitions from 1 to 0. • Roll-off rate: Slower for abrupt window discontinuities. Side Lobe Roll-off Rate Main Lobe

Rectangular Window • Narrow main lobe width: 4π/M, Side lobe width: 2 π /M

Rectangular Window • Narrow main lobe width: 4π/M, Side lobe width: 2 π /M • Spectral leakage directed to places that hurt analysis • Poor roll off and high initial lobe

Convoluting Rectangular Window with a particular frequency Note: Windows with more points narrow the

Convoluting Rectangular Window with a particular frequency Note: Windows with more points narrow the central lobe

Triangular (Bartlett Window)

Triangular (Bartlett Window)

Compare: Hamming to Hanning Hamming Hanning: Faster roll of; Hamming: lower first lobe

Compare: Hamming to Hanning Hamming Hanning: Faster roll of; Hamming: lower first lobe

Blackman & Hamming Frequency Response

Blackman & Hamming Frequency Response

Comparing Windows

Comparing Windows

Spectral Leakage using Hanning Window Understanding Digital Signal Processing, Third Edition, Richard Lyons (0

Spectral Leakage using Hanning Window Understanding Digital Signal Processing, Third Edition, Richard Lyons (0 -13 -261480 -4) © Pearson Education, 2011.

Spectral Leakage using Hamming Window Understanding Digital Signal Processing, Third Edition, Richard Lyons (0

Spectral Leakage using Hamming Window Understanding Digital Signal Processing, Third Edition, Richard Lyons (0 -13 -261480 -4) © Pearson Education, 2011.

Efficiency Issues • Calculation of Sin(x), Cos(x), ex – Sin(x) = x – x

Efficiency Issues • Calculation of Sin(x), Cos(x), ex – Sin(x) = x – x 3/3! + x 5/5! – x 7/7! … – Cos(x) = x – x 2/2! + x 4/4! – x 6/6! … – ex = x + x 2/2! + x 3/3! + x 3/4! … • Suggestions – Use table look up – x 3 = x * x – Consider the use of the cache when using loops

Fourier Analysis Real-time DSP Algorithms must be as efficient as possible • The Naïve

Fourier Analysis Real-time DSP Algorithms must be as efficient as possible • The Naïve DFT requires a double loop – Outer loop is T long (length of the time domain frame) – Inner loop is N long (number of FFT frequency bins) – Total calculations: O(N*T), or O(N 2) when N = T • Fast Fourier Transform (FFT) – N must be an even power of 2, if not, pad with zeros – Total calculations O(N lg. N) calculations • Compare efficiencies N N 2 lg N N lg N Improvement 128 16, 384 7 896 18. 2857 1, 024 1, 048, 576 10 10, 240 102. 4000 8, 192 67, 108, 864 13 106, 496 630. 1538

Understanding Digital Signal Processing, Third Edition, Richard Lyons (0 -13 -261480 -4) © Pearson

Understanding Digital Signal Processing, Third Edition, Richard Lyons (0 -13 -261480 -4) © Pearson Education, 2011.

Optimization of DFT 1. Replace the double loop by a recursive approach 2. Eliminate

Optimization of DFT 1. Replace the double loop by a recursive approach 2. Eliminate the recursion for further efficiency 3. Perform obscure optimizations to extract even more speed How do we know if the optimized version works? 1. Observe the frequency domain results a. b. c. d. Value at index N/2 -1 = value at index N/2+1 Value at index N/2 -k = value at index N/2+k The 0 index contains a constant multiplier that when removed causes the signal’s mean to be zero Value at N/2 is unrelated to the rest of the spectrum 2. Perform an inverse FFT to restore the original signal 3. Compare the results to the slower DFT algorithm 4. Apply algorithm to simple sinusoid signals and check results

Why is FFT fast • Evaluate polynomial at n roots of unity • Why

Why is FFT fast • Evaluate polynomial at n roots of unity • Why is this important? – Dealing with exponents is easy • eaeb=ea+b, (ea)b = eab – The ith root of unity ωin = e 2πi/n. – For even numbers, wnn/2=w 21=-1 – For even numbers, squares of the n complex roots of unity are the same as the n/2 complex roots of unity – Summing all n roots of unity equals zero – The calculation of e 4πi/n=e 2 π i/n because of the periodicity • We can avoid many calculations, and use a divide and conquer recursive solution – Solve smaller parts of the linear equations and easily combine results to solve the larger equation

The Fast Fourier Transform • • The FFT is an efficient divide and conquer

The Fast Fourier Transform • • The FFT is an efficient divide and conquer algorithm For polynomials, If n is even, we can divide into two smaller polynomials Peven(x) = a 0 + a 2 x 2 + a 4 x 4 + … + an-2 xn-2 Podd(x) = a 1 + a 3 x 2 + … + an-1 xn-1 and we can write the equation as follows P(x) = peven(x) + x * podd(x) • This suggests recursion of lg(n) steps. The above equation is the recursive relationship between levels

Theory for Optimization Base Case: x[0] Recursive Relationship (2 t)/N = t/(N/2) Fk =

Theory for Optimization Base Case: x[0] Recursive Relationship (2 t)/N = t/(N/2) Fk = ∑t=0 N-1 x[t] e-j 2πkt/N = ∑t=0 N/2 -1 x[2 t] e-j 2πk(2 t)/N + ∑ -j 2πk(2 t+1)/N t=0 N/2 -1 x[2 t+1] e = ∑t=0 N/2 -1 x[2 t] e-j 2πkt/(N/2) + e-j 2πk/N∑ -j 2πkt/(N/2) x[2 t+1] e t=0 N/2 -1 = Fkeven + e-j 2πk/N * Fkodd

The FFT Algorithm The recursive step has O(N) computations The recursion requires O(lg(N)) levels

The FFT Algorithm The recursive step has O(N) computations The recursion requires O(lg(N)) levels Total calculations are O(N lg(N))

FFT Recursive Efficiency Figure to the right – FFT of 16 points – O(N)

FFT Recursive Efficiency Figure to the right – FFT of 16 points – O(N) work at each level – lg N + 1 = 5 levels – Total calculations: O(N lg(N)

The FFT Algorithm Overview • Initial Recursive Call (FFT(time) time contains signal amplitudes, and

The FFT Algorithm Overview • Initial Recursive Call (FFT(time) time contains signal amplitudes, and is a power of 2 length (N) • Base Case: n=1, then return the time array • Divide Step separate time into even (timeeven) and odd (timeodd) sub-arrays • Recursive Step: yeven = FFT(timeeven), yodd= FFT(timeodd) • Combine Step For loop of complex number multiplications

Simple Recursive Solution public static Complex[] FFT(Complex[] time) { int N = time. length;

Simple Recursive Solution public static Complex[] FFT(Complex[] time) { int N = time. length; Complex[] fourier = new Complex[N]; if (N==1) { fourier[0] = time[0]; return fourier; } Complex[] even = new Complex[N/2], odd = new Complex[N/2]; for (int m=0; m<N/2; m++) { even[m] = time[2*m]; odd[m] = time[2*m+1]; } } Complex[] q = FFT(even), r = FFT(odd); for (int k=0; k<N/2; k++) { double exp = -2*k* Math. PI /N; Complex wk = new Complex(Math. cos(exp), Math. sin(exp)); fourier[k] = q[k]. plus(wk. times(r[k])); fourier[k+N/2] = q[k]. minus(wk. times(r[k])); } return fourier; -2 kπ/N+N/2 Note: e = -e

Inefficiencies Computations still are an order of magnitude slower than necessary • Repetitive calculations

Inefficiencies Computations still are an order of magnitude slower than necessary • Repetitive calculations of sines and cosines are extremely slow • Recursive activation record overhead is huge • Declaring and copying arrays at every step slows things down at least by half • N>>1 is ten times faster than N/2 • The Complex class methods cause machine code jumps that put pressure on the hardware cache

Eliminating the Recursion Butterfly algorithm 000 000 001 010 100 011 100 101 110

Eliminating the Recursion Butterfly algorithm 000 000 001 010 100 011 100 101 110 111 101 111 100 110 001 010 110 001 101 011 111 • The numbers in the rectangles are the array indices • Notice how the indices change during the recursion • Can you see a pattern ?

Bit flipping implementation • Determine if a number is an even power of two:

Bit flipping implementation • Determine if a number is an even power of two: n & (n-1) • Determine the upper most bit number, B int B=1, x = n-1; while (x&(x-1)) B++; Result = 0 FOR each bit b from 0 to B IF bit set, shift B – b and OR into result

Butterfly Code Flip bits from left to right int j = N>>1, k; for

Butterfly Code Flip bits from left to right int j = N>>1, k; for (int i=1; i<N-1; i++) { if (i < j) { swap (x[i], x[j]); } k = N>>1; while (k>=2 & j>=k) { j -= k; k >>= 1; } j += k; } • Most Significant Bit Swap. Bit ( x, x + lg. N) • Second most significant bit Swap. Bit(x, x + lg(N/2) • Third most significant bit Swap. Bit(x, x + lg(N/4) • kth most significant bit Swap. Bit(x, x + lg(N/2 k))

Sin and Cosine Table Look Up Compute the values ahead of time and save

Sin and Cosine Table Look Up Compute the values ahead of time and save repetitive calculations • ei 2πk/N = cos(2πk/N) + i sin(2πk/N) • We can store in an array (sin. X[]) sin(2π0/N), sin(2π1/N), sin(2π2/N) sin(2π3/N), … sin(2π(N-1)/N) • cos(2πk/N) = cos. X[(k+(N>>2))%N]

Wk. N= e 2πk/N This relationship allows table lookups for precomputed sines and cosines

Wk. N= e 2πk/N This relationship allows table lookups for precomputed sines and cosines Understanding Digital Signal Processing, Third Edition, Richard Lyons (0 -13 -261480 -4) © Pearson Education, 2011.

Optimized FFT // Perform the fft calculations. for (int stage=1; stage<=M; stage++) // M

Optimized FFT // Perform the fft calculations. for (int stage=1; stage<=M; stage++) // M = lg N { // Remember that complex numbers require pairs of doubles fft. Sub. Group. Gap = 2<<stage; // 4, 8, 16, . . . – subgroup distance gap = fft. Sub. Group. Gap>>1; // 2, 4, 8, . . . – odd/even distance k. Inc = N>>1; // Number of 2 PIki/N steps for odd/even entries. // Outer loop: each sub-fft group; inner loop: combine group elements for (int even=0; even<complex. length; even+=fft. Sub. Group. Gap) { k = 0; // Index into the trigonometric lookup table. for (int element=even; element<(even+gap); element+=2) { // ***** See Next Slide ***** k += k. Inc; } } k. Inc >>= 1; } // position for next look up.

Multiplication Portion // Look up e^2 PIki/N avoiding trig calculations here. real. W =

Multiplication Portion // Look up e^2 PIki/N avoiding trig calculations here. real. W = sines[(k+(N>>2))%N]; // cos(2 PIk/N); imag. W = -sines[k%N]; // -sin(2 PIk/N); // Complex multiplication of the odd entry of the subgroup // with (e^2 PIi/N)^k = (cos(2 PI/N) - i * sin(2*PI/N)^k j = (element + gap); temp. Real = real. W * complex[j] - imag. W * complex[j+1]; temp. Imag = real. W * complex[j+1] + imag. W * complex[j]; // Adjust the odd entry (subtract: the fft is periodic). complex[j] = complex[element] - temp. Real; complex[j+1] = complex[element+1] -temp. Imag; //Adjust the even entry. complex[element] += temp. Real; complex[element+1] += temp. Imag;

Final Notes • The above analysis – Assumes that N is a power of

Final Notes • The above analysis – Assumes that N is a power of 2 – If not, pad the array with zeroes – For speech recognition, we always use power of 2 frames • What if not a power of 2 – The same analysis is possible with other radices – When hit a prime, fall back to the slower O(n 2) algorithm

Final Notes • Standard Fast Fourier Transform – requires N to be a power

Final Notes • Standard Fast Fourier Transform – requires N to be a power of 2 for recursion to work – Can pad the array with zeroes to extend frequency domain • Can it work if N is not a power of 2? – Yes, but special slower processing is needed • How do we know if it works? – Point N/2 -1 = Point N/2+1, Point N/2 -2 = Point N/2+2, Point N/2 -k = Point N/2 + k, etc. – Note: Points 0 and N/2 don't match, so don’t check these – The FFT Inverse should restore the time domain signal – Compare to the slower correlation DFT calculation – Try some simple impulses and check the results