Signal Representations DS-GA 1013 / MATH-GA 2824 Mathematical Tools - - PowerPoint PPT Presentation

Signal Representations

DS-GA 1013 / MATH-GA 2824 Mathematical Tools for Data Science

https://cims.nyu.edu/~cfgranda/pages/MTDS_spring20/index.html Carlos Fernandez-Granda

Motivation

Limitation of frequency representation: no time resolution Limitation of Wiener filtering: cannot adapt to noisy signal

Windowing Short-time Fourier transform Multiresolution analysis Denoising via thresholding

Speech signal

0.9 1.0 1.1 1.2 1.3

Time (s)

7500 5000 2500 2500 5000 7500

Beyond Fourier

Problem: How to capture local fluctuations First segment signal, then compute DFT Naive segmentation: multiplication by rectangular window

Rectangular window

Rectangular window π ∈ CN with width 2w:

• π [j] :=
• 1

if |j| ≤ w,

• therwise.

Is this a good choice?

Signal

100 75 50 25 25 50 75 100

Time

1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

Window

100 75 50 25 25 50 75 100

Time

0.0 0.2 0.4 0.6 0.8 1.0

Windowed signal

100 75 50 25 25 50 75 100

Time

1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

DFT of windowed signal

100 75 50 25 25 50 75 100

Frequency

10 10 20 30 40

DFT of signal

100 75 50 25 25 50 75 100

Frequency

25 50 75 100 125 150 175 200

Multiplication in time is convolution in frequency

Let y := x1 ◦ x2 for x1, x2 ∈ CN The DFT of y equals ˆ y = 1 N ˆ x1 ∗ ˆ x2, ˆ x1 and ˆ x2 are the DFTs of x1 and x2 respectively

Proof

ˆ y [k] :=

N

• j=1

x1(j)x2(j) exp

• −i2πkj

N

• =

N

• j=1

1 N

n

• l=1

ˆ x1(l) exp i2πlj N

• x2(j) exp
• −i2πkj

N

• = 1

N

N

• l=1

ˆ x1(l)

N

• j=1

x2(j) exp

• −i2π(k − l)j

N

• = 1

N

N

• l=1

ˆ x1(l)ˆ x ↓l

2 [k]

Rectangular window

• π [j] :=
• 1

if |j| ≤ w,

• therwise,

DFT of rectangular window

100 75 50 25 25 50 75 100

Frequency

20 20 40 60 80

DFT of signal

100 75 50 25 25 50 75 100

Frequency

25 50 75 100 125 150 175 200

DFT of windowed signal

100 75 50 25 25 50 75 100

Frequency

10 10 20 30 40

Hann window

100 75 50 25 25 50 75 100

Time

0.0 0.2 0.4 0.6 0.8 1.0

Hann window

The Hann window h ∈ CN of width 2w equals h [j] := 1

2

• 1 + cos
• πj

w

• if |j| ≤ w,
• therwise

DFT of Hann window

100 75 50 25 25 50 75 100

Frequency

10 20 30 40

Signal

100 75 50 25 25 50 75 100

Time

1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

Hann window

100 75 50 25 25 50 75 100

Time

0.0 0.2 0.4 0.6 0.8 1.0

Windowed signal

100 75 50 25 25 50 75 100

Time

1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

DFT of signal

100 75 50 25 25 50 75 100

Frequency

25 50 75 100 125 150 175 200

DFT of Hann window

100 75 50 25 25 50 75 100

Frequency

10 20 30 40

DFT of windowed signal

100 75 50 25 25 50 75 100

Frequency

5 10 15 20

Time-frequency resolution

Time resolution governed by width of window Can we just make the window arbitrarily narrow?

Compressing in time dilates in frequency and vice versa

x ∈ L2 [−T/2, T/2] is nonzero in a band of width 2w around zero Let y be such that y(t) = x(αt), for all t ∈ [−T/2, T/2] , for some positive real number α such that w/α < T The Fourier series coefficients of y equal ˆ y [k] = 1 α

• x, φk/α

Proof

ˆ y [k] = T/2

t=−T/2

y(t) exp

• −i2πkt

T

• dt

= w/α

t=−w/α

x(αt) exp

• −i2πkt

T

• dt

= 1 α w

τ=−w

x(τ) exp

• −i2πkτ

αT

• dτ

= 1 α T/2

τ=−T/2

x(τ) exp

• −i2πkτ

αT

• dτ

w = 90

Time Frequency

100 75 50 25 25 50 75 100

Time

0.0 0.2 0.4 0.6 0.8 1.0 100 75 50 25 25 50 75 100

Frequency

20 40 60 80

w = 30

Time Frequency

100 75 50 25 25 50 75 100

Time

0.0 0.2 0.4 0.6 0.8 1.0 100 75 50 25 25 50 75 100

Frequency

5 10 15 20 25 30

w = 5

Time Frequency

100 75 50 25 25 50 75 100

Time

0.0 0.2 0.4 0.6 0.8 1.0 100 75 50 25 25 50 75 100

Frequency

1 2 3 4 5

Time-frequency resolution

Fundamental trade-off Uncertainty principle: cannot resolve in time and frequency simultaneously

Windowing Short-time Fourier transform Multiresolution analysis Denoising via thresholding

Short-time Fourier transform

• 1. Segment in overlapping intervals of length ℓ
• 2. Multiply by window vector
• 3. Compute DFT of length ℓ

Short-time Fourier transform

The short-time Fourier transform of x ∈ CN is STFT[ℓ](x)[k, s] :=

• x, ξ↓ s(1−αov)ℓ

k

• ,

0 ≤ k ≤ ℓ − 1, 0 ≤ s ≤ N (1 − αov)ℓ, ξk[j] :=

• w[ℓ](j) exp
• i2πkj

ℓ

• if 1 ≤ j ≤ ℓ
• therwise

Overlap between adjacent segments equals αovℓ

k = 3 s = 2 (N := 500, ℓ := 128, αov := 0.5)

Basis vector

100 200 300 400 500 Time 1.0 0.5 0.0 0.5 1.0

STFT coefficient

1 2 3 4 5 6 7 8 9

Time

10 20 30 40 50 60

Frequency

0.0 0.2 0.4 0.6 0.8 1.0

k = 3 s = 3 (N := 500, ℓ := 128, αov := 0.5)

Basis vector

100 200 300 400 500 Time 1.0 0.5 0.0 0.5 1.0

STFT coefficient

1 2 3 4 5 6 7 8 9

Time

10 20 30 40 50 60

Frequency

0.0 0.2 0.4 0.6 0.8 1.0

k = 8 s = 2 (N := 500, ℓ := 128, αov := 0.5)

Basis vector

100 200 300 400 500 Time 1.0 0.5 0.0 0.5 1.0

STFT coefficient

1 2 3 4 5 6 7 8 9

Time

10 20 30 40 50 60

Frequency

0.0 0.2 0.4 0.6 0.8 1.0

k = 23 s = 6 (N := 500, ℓ := 128, αov := 0.5)

Basis vector

100 200 300 400 500 Time 1.0 0.5 0.0 0.5 1.0

STFT coefficient

1 2 3 4 5 6 7 8 9

Time

10 20 30 40 50 60

Frequency

0.0 0.2 0.4 0.6 0.8 1.0

Matrix representation of STFT

         

F[ℓ]

· · · · · ·

F[ℓ]

· · · · · ·

F[ℓ]

· · · · · · · · · · · · · · · · · · · · · · · · · · ·                     diag

• w[ℓ]
• · · ·

· · · diag

• w[ℓ]
• · · ·

· · · diag

• w[ℓ]
• · · ·

· · · · · · · · · · · · · · · · · ·           x,

Computing the STFT

Number of segments: nseg := N/(1 − αov)ℓ Complexity of multiplication by window: nsegℓ Complexity of applying DFT: nsegℓ log ℓ (FFT) Total complexity: O(N log ℓ) (overlap is a fixed fraction)

Inverting the STFT

Apply inverse DFT to each segment Combine segments Same complexity

Speech signal (window length = 62.5 ms)

1 2 3 4 5 6

Time (s)

10000 5000 5000 10000

Speech signal (window length = 62.5 ms)

1 2 3 4 5 6

Time (s)

1000 2000 3000 4000

Frequency (Hz)

10

3

10

2

10

1

100 101 102 103

Windowing Short-time Fourier transform Multiresolution analysis Denoising via thresholding

Image

Vertical line (column 135)

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0

Multiresolution analysis

Scale / resolution at which information is encoded is not uniform Goal: Decompose signals into components at different resolutions

Multiresolution decomposition

Let N := 2K for some K, a multiresolution decomposition of CN is a sequence of nested subspaces VK ⊂ VK−1 ⊂ . . . ⊂ V0 satisfying: ◮ V0 = CN ◮ Vk is invariant to translations of scale 2k for 0 ≤ k ≤ K. If x ∈ Vk then x ↓ 2k ∈ Vk for all l ∈ Z ◮ For any x ∈ Vj that is nonzero only between 1 and N/2, the dilated vector x↔2 belongs to Vj+1

Dilation

Let x ∈ CN be such that x[j] = 0 for all j ≥ N/M, where M is a positive integer The dilation of x by a factor of M is x↔M[j] = x j M

How to build a multiresolution decomposition

◮ Set the coarsest subspace to be spanned by a low-frequency vector ϕ, called a scaling vector or father wavelet VK := span (ϕ) .

How to build a multiresolution decomposition

◮ Decompose the finer subspaces into the direct sum Vk := Wk ⊕ Vk+1, 0 ≤ k ≤ K − 1, where Wk captures the finest resolution available at level k ◮ Set Wk to be spanned by shifts of a vector µ dilated to have the appropriate resolution: Vk := Wk ⊕ Vk+1, 0 ≤ k ≤ K − 1, WK−1 :=

N−1 2k+1

• m=0

span

• µ ↓ m2k+1

↔2k

• .

The vector µ is called a mother wavelet

Challenge

How to choose mother and father wavelets? If chosen appropriately, basis vectors can be orthonormal

Haar wavelet basis

The Haar father wavelet ϕ is a constant vector, such that ϕ[j] := 1 √ N , 1 ≤ j ≤ N The mother wavelet µ satisfies µ[j] :=        − 1

√ 2,

j = 1,

1 √ 2,

j = 2, 0, j > 2 Other options: Meyer, Daubechies, coiflets, symmlets, etc.

Haar wavelets

Scale Basis functions 20 21 22

Multiresolution decomposition

VK W2 W1 W0 PVk x is an approximation of x at scale 2k

Vertical line (column 135)

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0

Scale 29

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

0.04 0.02 0.00 0.02 0.04 2 4 6 8 10 12 14 16

Scale 28

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

0.04 0.02 0.00 0.02 0.04 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Scale 27

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8

Scale 26

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

Scale 25

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

2 4 6 0.0 0.1 0.2 0.3 0.4 0.5 0.6

Scale 24

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

5 10 15 0.2 0.0 0.2 0.4 0.6 0.8

Scale 23

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

10 20 30 0.1 0.0 0.1 0.2 0.3 0.4

Scale 22

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

20 40 60 0.05 0.00 0.05 0.10 0.15 0.20 0.25

Scale 21

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

25 50 75 100 125 0.10 0.05 0.00 0.05 0.10 0.15 0.20

Scale 20

Approximation Coefficients

100 200 300 400 500 0.2 0.4 0.6 0.8 1.0 Data Approximation

50 100 150 200 250 0.2 0.1 0.0 0.1 0.2

Haar wavelets in the frequency domain

200 150 100 50 50 100 150 200

Frequency

0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16

Width: 200 Width: 100 Width: 50

Time-frequency support of basis vectors

STFT Wavelets

2D Wavelets

Extension to 2D by using outer products of 1D basis vectors To build a 2D basis vector at scale (m1, m2) and shift (s1, s2) we set v2D

[s1,s2,m1,m2] := v1D [s1,m1]

• v1D

[s2,m2]

T , where v1D can refer to 1D father or mother wavelets Nonseparable designs: steerable pyramid, curvelets, bandlets...

2D Haar wavelet basis vectors

Image

2D Haar wavelet decomposition

Approximation Coefficients

300 310 320 330 340 350

2D Haar wavelet decomposition

Approximation Coefficients

20 40 60 80 100 20 40 60 80 100 20 40 60 80 100

2D Haar wavelet decomposition

Approximation Coefficients

5 5 10 15 20 25 30 5 5 10 15 20 25 30 5 5 10 15 20 25 30

2D Haar wavelet decomposition

Approximation Coefficients

5 5 10 15 5 5 10 15 5 5 10 15

2D Haar wavelet decomposition

Approximation Coefficients

6 4 2 2 4 6 6 4 2 2 4 6 6 4 2 2 4 6

2D Haar wavelet decomposition

Approximation Coefficients

6 4 2 2 4 6 4 2 2 4 6 4 2 2 4

2D Haar wavelet decomposition

Approximation Coefficients

2 1 1 2 3 2 1 1 2 3 2 1 1 2 3

2D Haar wavelet decomposition

Approximation Coefficients

2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5

2D Haar wavelet decomposition

Approximation Coefficients

1.0 0.5 0.0 0.5 1.0 1.0 0.5 0.0 0.5 1.0 1.0 0.5 0.0 0.5 1.0

2D Haar wavelet decomposition

Approximation Coefficients

0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.6 0.4 0.2 0.0 0.2 0.4 0.6

Windowing Short-time Fourier transform Multiresolution analysis Denoising via thresholding

Denoising

Aim: Estimate signal x from data of the form y = x + z

Motivation

STFT coefficients of audio and wavelet coefficients of images are sparse Coefficients of noise are dense Idea: Get rid of small entries in Ay = Ax + A z

Why are noise coefficients dense?

If ˜ z is Gaussian with mean µ and covariance matrix Σ, then for any A, A˜ z is Gaussian with mean Aµ and covariance matrix AΣA∗ If A is orthogonal, iid zero mean noise is mapped to iid zero mean noise

Example

Data Signal

Thresholding

Hard-thresholding operator Hη (v) [j] :=

• v [j]

if |v [j]| > η

• therwise

Denoising via hard thresholding

Estimate Signal

Denoising via hard thresholding

Given data y and a sparsifying linear transform A

• 1. Compute coefficients Ay
• 2. Apply the hard-thresholding operator Hη : Cn → Cn to Ay
• 3. Invert the transform

xest := L Hη (Ay) , where L is a left inverse of A

Speech signal

4.30 4.35 4.40 4.45 4.50 4.55 4.60 4.65 4.70 Time (s) 7500 5000 2500 2500 5000 7500

STFT coefficients 2 4 6 Time (s) 1000 2000 3000 4000 Frequency (Hz) 10

5

10

4

10

3

10

2

10

1

Noisy signal

4.30 4.35 4.40 4.45 4.50 4.55 4.60 4.65 4.70 Time (s) 7500 5000 2500 2500 5000 7500 10000

STFT coefficients 2 4 6 Time (s) 1000 2000 3000 4000 Frequency (Hz) 10

4

10

3

10

2

10

1

Thresholded STFT coefficients 2 4 6 Time (s) 1000 2000 3000 4000 Frequency (Hz) 10

5

10

4

10

3

10

2

10

1

Denoised signal

4.30 4.35 4.40 4.45 4.50 4.55 4.60 4.65 4.70 Time (s) 7500 5000 2500 2500 5000 7500

Denoised signal

4.3600 4.3625 4.3650 4.3675 4.3700 4.3725 4.3750 4.3775 Time (s) 3000 2000 1000 1000 2000 3000 Signal STFT thresholding Noisy data

Denoised signal (Wiener filtering)

4.3600 4.3625 4.3650 4.3675 4.3700 4.3725 4.3750 4.3775 Time (s) 3000 2000 1000 1000 2000 3000 Signal Wiener denoising Noisy data

Image

Wavelet coefficients

Noisy signal

Wavelet coefficients

Thresholded wavelet coefficients

Denoised signal

Comparison

Clean Noisy Wiener filtering Wavelet thresholding

Coefficients are structured 2 4 6 Time (s) 1000 2000 3000 4000 Frequency (Hz) 10

5

10

4

10

3

10

2

10

1

Coefficients are structured

Block thresholding

Assumption: Coefficients are group sparse, nonzero coefficients cluster together Threshold according to block of surrounding coefficients Ij Bη (v) [j] :=

• v [j]

if j ∈ Ij such that

• vIj
• 2 > η, ,
• therwise,

Denoising via block thresholding

Given data y and a sparsifying linear transform A

• 1. Compute coefficients Ay
• 2. Apply the block-thresholding operator Hη : Cn → Cn to Ay
• 3. Inverting the transform

xest := L Bη (Ay) , where L is a left inverse of A

Noisy STFT coefficients 2 4 6 Time (s) 1000 2000 3000 4000 Frequency (Hz) 10

4

10

3

10

2

10

1

Thresholded STFT coefficients 2 4 6 Time (s) 1000 2000 3000 4000 Frequency (Hz) 10

5

10

4

10

3

10

2

10

1

Block-thresholded STFT coefficients (block of length 5) 2 4 6 Time (s) 1000 2000 3000 4000 Frequency (Hz) 10

5

10

4

10

3

10

2

10

1

Thresholding

4.3600 4.3625 4.3650 4.3675 4.3700 4.3725 4.3750 4.3775 Time (s) 3000 2000 1000 1000 2000 3000 Signal STFT thresholding Noisy data

Block thresholding

4.3600 4.3625 4.3650 4.3675 4.3700 4.3725 4.3750 4.3775 Time (s) 3000 2000 1000 1000 2000 3000 Signal STFT block thresholding Noisy data

Noisy wavelet coefficients

Thresholded wavelet coefficients

Block thresholded wavelet coefficients

Denoised signal

Comparison

Clean Noisy Wiener filtering Wavelet thresholding Wavelet block thresholding