Introduction to Time-Frequency Analysis Speaker Po-Hong Wu

An Introduction to Time-Frequency Analysis Speaker: Po-Hong Wu Advisor: Jian-Jung Ding Digital Image and Signal Processing Lab GICE, National Taiwan University NTU, GICE, MD531, DISP Lab 1

Outline Introduction Short-Time Fourier Transform Gabor Transform Wigner Distribution Function Spectrogram S Tranform Cohen s Class Time-Frequency Distribution Fractional Fourier Transform Motion on Time-Frequency Distributions Hilbert-Huang Transform Conclusion Reference NTU, GICE, MD531, DISP Lab 2

Introduction Fourier transform (FT) ( ) X f x t e ( ) = t varies from ~ 2 j f t dt Time-Domain Frequency Domain 1 2 f(t) 0.5 1 Fourier transform 0 0 -0.5 -1 [A1] -1 -2 0 5 10 15 20 25 30 -5 0 5 Why do we need time-frequency transform? NTU, GICE, MD531, DISP Lab 3

Example: x(t) = cos( t) when t < 10, x(t) = cos(3 t) when 10 t < 20, x(t) = cos(2 t) when t 20 -5 -4 -3 -2 -1 0 1 2 [B2] 4 3 5 0 5 10 15 20 25 30 NTU, GICE, MD531, DISP Lab 4

Short Time Fourier Transform ( ) ( ) ( ) x = 2 j f , X t f w t e d w(t): mask function windowed Fourier transform or time-dependent Fourier transform NTU, GICE, MD531, DISP Lab 5

When w(t) is a rectangular function w(t) = 1 for |t| B , w(t) = 0 , otherwise ( ) , t B t B + ( ) = 2 j f X t f x e d STFT 14 12 10 8 Hz 6 4 [B3] 0 2 0 1 2 3 4 5 6 7 8 9 10 11 sec NTU, GICE, MD531, DISP Lab 6

Advantage: less computation time Disadvantage: worse representaion Application: deal with large data Ex: real time processing NTU, GICE, MD531, DISP Lab 7

Gabor Transform A specail case of the STFT where ( ) , x G t f e = 2 ( ) exp( ) w t t ( ) = 2 ( ) 2 t j f e x d Other definition 2 ( ) t t 1 ( ) j ( ) ( ) = , G t e e x d 2 2 x 2 Gabor 14 12 10 8 Hz 6 4 [B4] 0 1 2 3 4 5 6 sec 2 0 7 8 9 10 11 NTU, GICE, MD531, DISP Lab 8

Why do we choose the Guassian function? Among all functions of w(t), the Gaussian function has area in time-frequency distribution is minimal than other STFT. Gaussian function is an eigenfunction of Fourier transform, so the Gabor transform has the same properties in time domain and in frequency domain. NTU, GICE, MD531, DISP Lab 9

Approximation of the Gabor Transform 2 Because of when |a|>1.9143 0.00001 e a 1.9143 + t ( ) ( ) 2 ( ) 2 t j f , G t f e e x d x 1.9143 t 2/2 Because of when |a|>4.7985 0.00001 e a 2 ( ) t t 1 + ( ) 4.7985 j t ( ) ( ) = , G t e e x d 2 2 ,3 x 2 4.7985 t NTU, GICE, MD531, DISP Lab 10

Generalization of the Gabor Transform ( ) ( ) = 2 ( ) 2 t j f , G t f e e x d x For larger : higher resolution in the time domain but lower resolution in the frequency domain For smaller : higher resolution in the frequency domain but lower resolution in the time domain NTU, GICE, MD531, DISP Lab 11

Resolution Using the generalized Gabor transform with larger Using other time unit instead of second NTU, GICE, MD531, DISP Lab 12

Wigner Distribution Function ( ) ( ) , /2 x W t f x t x t ( ) = + * 2 j f /2 e d Other definition = + * 2 j ( , ) ( / 2) ( / 2) W t f X f X f e d x Wigner distribution 6 4 frequency (Hz) 2 0 -2 [B5] -4 -6 0 5 10 15 20 25 30 time (sec) NTU, GICE, MD531, DISP Lab 13

Signal auto-correlation function ( ) ( , x C t x t = + ) ( ) / 2 / 2 x t Spectrum auto-correlation function ( ) ( , f X f = ) ( ) + /2 /2 x S X f Ambiguity function (AF) ( ) , ( ) ( ) = + x t * 2 j t /2 /2 x A x t e dt IFTf FTt Cx(t, ) IFTf FTt Ax( , ) Wx(t, f ) [B6] Sx( , f ) FTt IFTf NTU, GICE, MD531, DISP Lab 14

Modified Wigner Distribution Wigner Ville Distribution For compressing inner interference = + * / 2) ( 2 j f ( , ) ( / 2) W t f x t x t e d x Analytic signal ( ) ( ) x t x t = + ( ) jx t NTU, GICE, MD531, DISP Lab 15

Pseudo Wigner Distribution For surpressing outer interference ( , ) Y t f Pseudo L-Wigner distribution 6 = ( / 2) ( w / 2) ( + * 2 j f ( , ) / 2) ( / 2) W t f w x t x t e d x = + = + * ( , / 2) ( , / 2) ) e Y t f Y t f ( ) ( w d where 2 j f x t d 4 2 frequency (Hz) 0 [B7] -2 -4 -6 0 5 10 15 20 25 30 NTU, GICE, MD531, DISP Lab 16 time (s)

Gabor-Wigner Distribution ( , ) ( ( , ), f f C t w p G t w W t w = ( , )), f 10 10 (a) (b) 5 5 0 0 -5 -5 [B8] ( , ) p x y xy = -10 -10 -10 -5 0 5 10 -10 -5 0 5 10 2 = ( , ) p x y min( , ) y x NTU, GICE, MD531, DISP Lab 17

Spectrogram = 2 ( ) ( 2 j f ( , ) ) SP t f x h t e d Another form = ( , ) d d ( , ) ( , ) SP t f W t f W h x Spectrogram 14 12 10 8 Hz 6 4 2 [B9] 0 0 1 2 3 4 5 6 7 8 9 10 11 sec NTU, GICE, MD531, DISP Lab 18

S-Transform Original S-Transform = f 2 2 ( ) f t ( , ) 2 i ft ( ) x t exp[ ] ST f e dt 2 2 f 2 2 ( ) f t exp[ ] Where w(t)= 2 2 S Transform 4 3 cos( cos(3 cos(2 t) ) ) t when when 1 10 20 10 t when t 2 = 1 frequency ( ) x t 20 30 t 0 t -1 -2 -3 [B10] -4 0 5 10 15 20 25 30 time (sec) NTU, GICE, MD531, DISP Lab 19

Generalized S-Transform ( , , ) GS f p = 2 i t ( ) ( x t w , , ) t f p e dt Another definition = + ( , , ) f p e 2 i ( , , ) f p ( ) GS X f W d Ristriction = ( , , ) t f p d 1 w NTU, GICE, MD531, DISP Lab 20

Novel S-Transform with the Special Varying Window 2 ( ) 2 F f ( )( F f ) t = ( , ) t f exp( ) w S 2 Restriction = ( , ) t f d 1 s w = 2 ( ) 1/ F f When , it becomes the Gabor transform. When , it becomes the original S-trnasform. ( ) F f f = 2 NTU, GICE, MD531, DISP Lab 21

Cohens Class Time-Frequency Distribution = ( , ) ( , )exp( 2 ( x d d ( , ) t f )) C A j t f x Ambiguity function ( , ) = / 2) ( + * 2 j t ( / 2) x A x t x x e dt C t ( , ) x FTt IFTf FTt IFTf x A ( , ) ( , ) WD t f + ( / 2) ( / 2) x t x t x [B11] FTt ( , ) IFTf x S f NTU, GICE, MD531, DISP Lab 22

For the ambiguity function The auto terms are always near to the origin. The cross terms are always from the origin. [B12] NTU, GICE, MD531, DISP Lab 23

( , ) Kernel function Choi-Williams Distribution ( ) , exp Choi-Williams distribution 10 ( ) 2 = 8 6 4 2 eta 0 -2 [B13] -6 -4 -8 -10 -15 -10 -5 0 5 10 15 24 NTU, GICE, MD531, DISP Lab tau (sec)

Cone-Shape Distribution ( , ) sin ( c = 2 >0 2 )e Cone Shape distribution 10 8 6 4 2 eta 0 -2 -4 -6 -8 -10 -15 -10 -5 0 5 10 15 tau (sec) NTU, GICE, MD531, DISP Lab 25

Fractional Fourier Transform ( ) FT x t X f = ( ) ( ) FT FT x t x t = ( ) ( ) FT FT FT x t X f = = ( ) ( ) FT FT FT FT x t x t ( ) ( ) ( ) IFT f t ( ) = How to rotate the time-frequency distribution by the angle other than /2, , and 3 /2? 26 NTU, GICE, MD531, DISP Lab

Zero rotation: 0 R = I Consistency with Fourier transform: = FT Additivity of rotation: R R R = R /2 + rotation: 2 R I = NTU, GICE, MD531, DISP Lab 27

( ) x t dt e e e 2 2 csc cot ( ) u j j t t 2 u = cot j 1 cot X u j 2 2 2 = = = f(t): rectangle 1 1 1 0 0 0 -1 -1 -1 -5 0 5 -5 0 5 -5 0 5 2 2 2 = = = 1 1 1 0 0 0 F(w): sinc function [A3] -1 -1 -1 -5 0 5 -5 0 5 -5 0 5 NTU, GICE, MD531, DISP Lab 28

Application Decomposition in the time-frequency distribution 1.5 2 x(t) = signal + noise Fourier transform of x(t) 1 1 0.5 0 0 (non-separable) -0.5 -1 -10 -5 0 5 10 -10 -5 0 5 10 1.5 2 fractional Fourier transform of x(t) = recovered signal 1 1 0.5 0 0 (separable) -1 -0.5 -10 -5 0 5 10 -10 -5 0 5 10 NTU, GICE, MD531, DISP Lab 29

f-axis Signa l noise noise noise Signal Signa l FRFT FRFT t-axis cutoff line cutoff line NTU, GICE, MD531, DISP Lab 30

Modulation and Multiplexing 5 5 unfilled T-F slot 0 0 -5 -5 -20 0 20 -20 0 20 (c) WDF of G(u) (d) GWT of G(u) NTU, GICE, MD531, DISP Lab 31

Time domain Frequency domain fractional domain Modulation Shifting Modulation + Shifting Shifting Modulation Modulation + Shifting Differentiation j2 f Differentiation and j2 f j2 f Differentiation Differentiation and j2 f NTU, GICE, MD531, DISP Lab 32

Motion on Time-Frequency Distributions Horizontal Shifting 2 j f t 0 ( ) ( , ) ,STFT, Gabor x t t S t t f e 0 0 x ( , ) ,Wigner W t t f 0 x Vertical Shifting 2 j f t 0 ( ) x t ( , t f ) ,STFT,Gabor e S f 0 x ( , t f ) ,Wigner W f 0 x NTU, GICE, MD531, DISP Lab 33

Dilation 1 x a t t ( ) ( , ) ,STFT,Gabor S af x | | a a t ( , ) ,WDF W af x a Shearing ( ) x t e = 2 j at ( ) y t ( , ) t f ( , t f ) ,STFT,Gabor S S at x y = ( , ) t f ( , t f ) ,WDF W W 2 t a at x y j = ( ) x t ( ) y t e ( , ) t f ( , ) ,STFT,Gabor af f S S t x y = ( , ) t f ( , ) ,WDF af f W W t x y NTU, GICE, MD531, DISP Lab 34

Rotation If F{x(t)}=X(f), then F{X(t)}=x(-f). We can derive: | ( , )| | ( , )| ,STFT S t f S f t x X j G t f G f t e x X W t f W f t x X 2 ft = ( , ) ( , ) ,Gabor = ( , ) ( , ) ,WDF NTU, GICE, MD531, DISP Lab 35

Hilbert-Huang Transform Introduction Most of distribution are designed for stationary and linear signals, but, In the real world, most of signals are non-stationary and non-linear. HHT consists two parts: empirical mode decomposition (EMD) Hilbert spectral analysis (HSA) NTU, GICE, MD531, DISP Lab 36

Empirical decomposition function Any complicated data can be decomposed into a finite and small number of intrinsic mode functions (IMF) by sifting processing. Intrinsic mode function (1)In the whole data set, the number of extrema and the number of zero-crossing must either equal or differ at most by one. (2)At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. NTU, GICE, MD531, DISP Lab 37

Sifting Process (1) First, find all the local maxima extrema of x(t). IMF 1; iteration 0 2 1 0 -1 -2 10 20 30 40 50 60 70 80 90 100 110 120 (2) Interpolate (cubic spline fitting) between all the maxima extrema ending up with some upper envelope . e t max( ) IMF 1; iteration 0 2 1 0 -1 NTU, GICE, MD531, DISP Lab 38 -2 10 20 30 40 50 60 70 80 90 100 110 120

(3) Find all the local minima extrema. IMF 1; iteration 0 2 1 0 -1 -2 10 20 30 40 50 60 70 80 90 100 110 120 (4) Interpolate (cubic spline fitting) between all the minima extrema ending up with some lower envelope . e t min( ) IMF 1; iteration 0 2 1 0 -1 -2 10 20 30 40 50 60 70 80 90 100 110 120 NTU, GICE, MD531, DISP Lab 39

(5) Compute the mean envelope between upper envelope and lower envelope. + ( ) t ( ) t min e max e = ( ) m t 2 IMF 1; iteration 0 2 1 0 -1 -2 10 20 30 40 50 60 70 80 90 100 110 120 = ( ) ( ) x t ( ) h t m t (6) Compute the residue residue 1.5 1 0.5 0 -0.5 -1 -1.5 10 20 30 40 50 60 70 80 90 100 110 120 NTU, GICE, MD531, DISP Lab 40

(7) Repeat the above procedure (step (1) ~ step (6)) on the residue until the residue is a monotonic function or constant. The original signal equals the sum of the various IMFs plus the residual trend. n x t c t r t = = + ( ) ( ) ( ) k n 1 k NTU, GICE, MD531, DISP Lab 41

EX: 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 1 2 3 IMF1 4 5 6 x 104 0.1 0.05 0 -0.05 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 IMF2 0.2 0.1 0 -0.1 -0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Time NTU, GICE, MD531, DISP Lab 42

IMF3 0.4 0.2 0 -0.2 -0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 IMF4 0.4 0.2 0 -0.2 -0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Time IMF5 0.1 0.05 0 -0.05 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 IMF6 0.04 0.02 0 -0.02 -0.04 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 NTU, GICE, MD531, DISP Lab 43 Time

IMF7 0.01 0.005 0 -0.005 -0.01 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 IMF8 0.01 0.005 0 -0.005 -0.01 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Time IMF9 x 10-3 10 5 0 -5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x 10-3 IMF10 10 5 0 -5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Time NTU, GICE, MD531, DISP Lab 44

-3 IMF11 x 10 4 2 0 -2 -4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -3 IMF12 x 10 0 -0.5 -1 -1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Time NTU, GICE, MD531, DISP Lab 45

Hilbert Spectral Anaysis ( ) 1 1 x t = = = ( ) y t { ( )} H x t ( ) x t PV d t a t e = + = ( ) t i ( ) z t ( ) x t ( ) ( ) iy t ( ) ( ) x t y t = = 1 ( ) t arg( ( )) z t tan ( ) 1 1 ( ) t dt d = = ( ) ( ) t f t 2 2 NTU, GICE, MD531, DISP Lab 46

Advantage Disadvantage STFT 1. Low computation 1. Complex value and 2. The range of the integration is limited 2. Low resolution Gabor transform 3. No cross term 4. Linear operation Wigner distribution 1. Real 1. High computation function 2. High resolution 2. Cross term 3. If the time/frequency limited, time/frequency of the 3. Non-linear operation WDF is limited with the same range Cohen s 1. Avoid the cross term 1. High computation class 2. Higher clarity 2. Lack of well mathematical properties distribution Gabor-Wigner distribution 1. Combine the advantage of the WDF and the Gabor 1. High computation function transform 2. Higher clarity 3. No cross-term NTU, GICE, MD531, DISP Lab 47

Conclusion We introduce many distributions here and put most attention on computation time and representations. We can find that the representation with higher clarity cost more computation time for all methods. Resolution Computation time The Hilbert-Huang transform is the most power method to deal with non-linear and non-stationary signals but lacks of physical background. NTU, GICE, MD531, DISP Lab 48

Reference [1][A]J. J. Ding, Time-Frequency Analysis and Wavelet Transform, National Taiwan University, 2009. [Online].Available: http://djj.ee.ntu.edu.tw/TFW.htm. [2][B]W. F. Wang, Time-Frequency Analyses and Their Fast Implementation Algorithm, Master Thesis, National Taiwan University, June, 2009. [3]Luis B. Almeida, Member, IEEE, The Fractional Fourier Transform and Time-Frequency Representations, IEEE Transaction On Signal Processing, vol. 42, no. 11, November 1994. [4]M. R. Spiegel, Mathematical Handbook of Formulas and Tables, McGraw- Hill, 1990. [5]N. E. Huang, Z. Shen and S. R. Long, et al., The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time Series Analysis , Proc. Royal Society, vol. 454, pp.903-995, London, 1998. NTU, GICE, MD531, DISP Lab 49