Spectral super-resolution DS-GA 1013 / MATH-GA 2824 - PowerPoint PPT Presentation

Spectral super-resolution DS-GA 1013 / MATH-GA 2824 Optimization-based Data Analysis http://www.cims.nyu.edu/~cfgranda/pages/OBDA_fall17/index.html Carlos Fernandez-Granda

Spectral super-resolution Aim: Estimate frequencies of multisinusoidal signal s � g ( t ) := � c [ j ] exp ( − i 2 π f j t ) j = 1 where f 1 , f 2 , . . . , f s ∈ [ − 1 / 2 , 1 / 2 ] from n of samples g ( − ( n − 1 ) / 2 ) , g ( − ( n − 1 ) / 2 + 1 ) , . . . , g (( n − 1 ) / 2 )

Spectral super-resolution Signal Data

Dirac measure It’s a measure, not a function! � 1 if τ ∈ S � δ [ τ ] ( d u ) = 0 otherwise S For any function h we have � h ( τ ) if τ ∈ S � h ( u ) δ [ τ ] ( d u ) = 0 otherwise S

Spectral super-resolution Aim: Estimating the measure s � � µ g := c [ j ] δ [ f j ] j = 1 from n Fourier coefficients � 1 / 2 h k ( u ) µ g ( d u ) − 1 / 2

Spectral super-resolution Aim: Estimating the measure s � � µ g := c [ j ] δ [ f j ] j = 1 from n Fourier coefficients � 1 / 2 � 1 / 2 s � h k ( u ) µ g ( d u ) = c [ j ] � h k ( u ) δ [ f j ] ( d u ) − 1 / 2 − 1 / 2 j = 1

Spectral super-resolution Aim: Estimating the measure s � � µ g := c [ j ] δ [ f j ] j = 1 from n Fourier coefficients � 1 / 2 � 1 / 2 s � h k ( u ) µ g ( d u ) = � c [ j ] h k ( u ) δ [ f j ] ( d u ) − 1 / 2 − 1 / 2 j = 1 s � = � c [ j ] exp ( − i 2 π kf j ) j = 1

Spectral super-resolution Aim: Estimating the measure s � � µ g := c [ j ] δ [ f j ] j = 1 from n Fourier coefficients � 1 / 2 � 1 / 2 s � h k ( u ) µ g ( d u ) = c [ j ] � h k ( u ) δ [ f j ] ( d u ) − 1 / 2 − 1 / 2 j = 1 s � = � c [ j ] exp ( − i 2 π kf j ) j = 1 − n − 1 ≤ k ≤ n − 1 = g ( k ) , 2 2

The periodogram Prony’s method Subspace methods

Periodogram y ∈ C n is The periodogram of � n − 1 2 � P � y ( u ) := � y [ k ] h k ( u ) k = − n − 1 2

Spectrum of Dirichlet kernel D − k c k c

Dirichlet kernel d 2 k c + 1 1 1 − 2 k c +1 2 k c +1 − 1 / 2 1 / 2

Periodogram as a convolution The periodogram of g ( − ( n − 1 ) / 2 ) , g ( − ( n − 1 ) / 2 + 1 ) , . . . , g (( n − 1 ) / 2 ) equals s � P g ( u ) = � c [ j ] d [ f j ] ( u ) j = 1

Time shift The τ -shifted version of a function f ∈ L 2 [ − 1 / 2 , 1 / 2 ] is f [ τ ] ( t ) := f ( t − τ ) where the shift is circular (it wraps around) For any shift τ F [ τ ] [ k ] = exp ( − i 2 π k τ ) F [ k ]

Proof The Fourier coefficients of the shifted Dirichlet kernel equal � exp ( − i 2 π kf ) if | k | ≤ ( n − 1 ) / 2 D [ f ] [ k ] := 0 otherwise

Proof s � c [ j ] exp ( − i 2 π kf j ) g ( k ) := � j = 1

Proof s � c [ j ] exp ( − i 2 π kf j ) g ( k ) := � j = 1 s � = � c [ j ] D [ f j ] [ k ] j = 1

Proof s � c [ j ] exp ( − i 2 π kf j ) g ( k ) := � j = 1 s � = c [ j ] D [ f j ] [ k ] � j = 1 n − 1 2 � P g ( u ) = g ( k ) h k ( u ) k = − n − 1 2

Proof s � c [ j ] exp ( − i 2 π kf j ) g ( k ) := � j = 1 s � = � c [ j ] D [ f j ] [ k ] j = 1 n − 1 2 � P g ( u ) = g ( k ) h k ( u ) k = − n − 1 2 n − 1 s 2 � � = � c [ j ] D [ f j ] [ k ] h k ( u ) j = 1 k = − n − 1 2

Proof s � c [ j ] exp ( − i 2 π kf j ) g ( k ) := � j = 1 s � = � c [ j ] D [ f j ] [ k ] j = 1 n − 1 2 � P g ( u ) = g ( k ) h k ( u ) k = − n − 1 2 n − 1 s 2 � � = c [ j ] � D [ f j ] [ k ] h k ( u ) j = 1 k = − n − 1 2 s � = � c [ j ] d [ f j ] ( u ) j = 1

Periodogram as a convolution Periodogram Data 0 5 10

Periodogram as a convolution Periodogram Data 0 5 10 15 20

Periodogram as a convolution Periodogram Data 0 10 20 30 40

Periodogram as a convolution Periodogram Data 0 10 20 30 40 50

Problem Signal (magnitude) Periodogram

Windowed periodogram y ∈ C n is The windowed periodogram of � n − 1 2 � y ( u ) := W [ k ] g ( k ) h k ( u ) . P w ,� k = − n − 1 2 � n − 1 − n − 1 � � � W , . . . , W are the Fourier coefficients of a 2 2 window function

Windowing Data Window function Windowed data

Windowed periodogram n − 1 2 � y ( u ) = W [ k ] g ( k ) h k ( u ) P w ,� k = − n − 1 2

Windowed periodogram n − 1 2 � y ( u ) = W [ k ] g ( k ) h k ( u ) P w ,� k = − n − 1 2 n − 1 s 2 � � = � c [ j ] W [ k ] exp ( − i 2 π f j k ) h k ( u ) j = 1 k = − n − 1 2

Windowed periodogram n − 1 2 � y ( u ) = W [ k ] g ( k ) h k ( u ) P w ,� k = − n − 1 2 n − 1 s 2 � � = � c [ j ] W [ k ] exp ( − i 2 π f j k ) h k ( u ) j = 1 k = − n − 1 2 n − 1 s 2 � � = � c [ j ] W [ f j ] [ k ] h k ( u ) j = 1 k = − n − 1 2

Windowed periodogram n − 1 2 � y ( u ) = W [ k ] g ( k ) h k ( u ) P w ,� k = − n − 1 2 n − 1 s 2 � � = � c [ j ] W [ k ] exp ( − i 2 π f j k ) h k ( u ) j = 1 k = − n − 1 2 n − 1 s 2 � � = � c [ j ] W [ f j ] [ k ] h k ( u ) j = 1 k = − n − 1 2 s � = � c [ j ] w [ f j ] ( u ) j = 1

Problem Signal (magnitude) Periodogram

Windowed periodogram Signal (magnitude) Windowed periodogram

Periodogram ∆ = 1 . 2 ∆ = 2 . 4 n n

Gaussian periodogram ∆ = 1 . 2 ∆ = 2 . 4 n n

The periodogram Prony’s method Subspace methods

Prony polynomial Signal (magnitude) Prony polynomial (magnitude)

Prony polynomial Given any f 1 , f 2 , . . . , f s , there exists a nonzero complex polynomial of order s s � P [ k ] z k p ( z ) := k = 0 such that its s roots are equal to exp ( i 2 π f 1 ) , exp ( i 2 π f 2 ) , . . . , exp ( i 2 π f s )

Proof s � ( 1 − exp ( − i 2 π f j ) z ) p ( z ) := j = 1

Proof s � ( 1 − exp ( − i 2 π f j ) z ) p ( z ) := j = 1 is of order s, so has at most s roots

Proof s � ( 1 − exp ( − i 2 π f j ) z ) p ( z ) := j = 1 is of order s, so has at most s roots s � P [ k ] z k p ( z ) = 1 + k = 1

Proof s � ( 1 − exp ( − i 2 π f j ) z ) p ( z ) := j = 1 is of order s, so has at most s roots s � P [ k ] z k p ( z ) = 1 + k = 1 is nonzero since p ( 0 ) = 1

Proof s � ( 1 − exp ( − i 2 π f j ) z ) p ( z ) := j = 1 is of order s, so has at most s roots s � P [ k ] z k p ( z ) = 1 + k = 1 is nonzero since p ( 0 ) = 1 Reveals the frequencies p ( exp ( i 2 π f j )) = 0 1 ≤ j ≤ s

Prony system Let s − n − 1 ≤ k ≤ n − 1 � � c [ j ] exp ( − i 2 π kf j ) g ( k ) := 2 2 j = 1 For any integer b s � P [ l ] g [ l − b ] = 0 l = 0 The equation only involves the data as long as − n − 1 ≤ − b ≤ s − b ≤ n − 1 2 2

Proof s � P [ l ] g [ l − b ] l = 0

Proof � 1 / 2 s s � � P [ l ] g [ l − b ] = P [ l ] h l − b ( u ) µ g ( d u ) − 1 / 2 l = 0 l = 0

Proof � 1 / 2 s s � � P [ l ] g [ l − b ] = P [ l ] h l − b ( u ) µ g ( d u ) − 1 / 2 l = 0 l = 0 � 1 / 2 s � = exp ( i 2 π bu ) P [ l ] exp ( − i 2 π lu ) µ g ( d u ) − 1 / 2 l = 0

Proof � 1 / 2 s s � � P [ l ] g [ l − b ] = P [ l ] h l − b ( u ) µ g ( d u ) − 1 / 2 l = 0 l = 0 � 1 / 2 s � = exp ( i 2 π bu ) P [ l ] exp ( − i 2 π lu ) µ g ( d u ) − 1 / 2 l = 0 � 1 / 2 = exp ( i 2 π bu ) p ( exp ( − i 2 π u )) µ g ( d u ) − 1 / 2

Proof � 1 / 2 s s � � P [ l ] g [ l − b ] = P [ l ] h l − b ( u ) µ g ( d u ) − 1 / 2 l = 0 l = 0 � 1 / 2 s � = exp ( i 2 π bu ) P [ l ] exp ( − i 2 π lu ) µ g ( d u ) − 1 / 2 l = 0 � 1 / 2 = exp ( i 2 π bu ) p ( exp ( − i 2 π u )) µ g ( d u ) − 1 / 2 s � � c [ j ] exp ( i 2 π bf j ) p ( exp ( − i 2 π f j )) = j = 1

Proof � 1 / 2 s s � � P [ l ] g [ l − b ] = P [ l ] h l − b ( u ) µ g ( d u ) − 1 / 2 l = 0 l = 0 � 1 / 2 s � = exp ( i 2 π bu ) P [ l ] exp ( − i 2 π lu ) µ g ( d u ) − 1 / 2 l = 0 � 1 / 2 = exp ( i 2 π bu ) p ( exp ( − i 2 π u )) µ g ( d u ) − 1 / 2 s � � c [ j ] exp ( i 2 π bf j ) p ( exp ( − i 2 π f j )) = j = 1 = 0

Prony’s method 1. Solve the system of equations  g ( 1 ) g ( 2 ) · · ·   g ( 0 )  g ( s ) g ( 0 ) g ( 1 ) · · · g ( s − 1 ) g ( − 1 )  �     P = −     · · · · · · · · · · · · · · ·    g ( − s + 2 ) g ( − s + 3 ) · · · g ( 1 ) g ( − s + 1 ) 2. Compute the roots z 1 , . . . , z s of the polynomial s � � P [ k ] z k p ( z ) := 1 + k = 1 3. For every root on the unit circle z j = exp ( i 2 πτ ) include τ in the set of estimated frequencies

Without noise it works! Signal (magnitude) Prony polynomial (magnitude)