A new sparsity-enabled signal separation method based on signal resonance PDF Free Download

1 / 35
0 views35 pages

A new sparsity-enabled signal separation method based on signal resonance PDF Free Download

A new sparsity-enabled signal separation method based on signal resonance PDF free Download. Think more deeply and widely.

A NEW SPARSITY-ENABLED SIGNAL SEPARATION METHOD
BASED ON SIGNAL RESONANCE
Ivan Selesnick
Polytechnic Institute of New York University
Brooklyn, New York
1
Introduction
Problem: Decomposition of a signal into the sum of two components:
1. Oscillatory (rhythmic, tonal) component
2. Transient (non-oscillatory) component
Outline
1. Signal resonance and Q-factors
2. Morphological component analysis (MCA)
3. Rational-dilation wavelet transform (RADWT)
4. Split augmented Lagrangian shrinkage algorithm (SALSA)
5. Examples
References (MDCT, etc)
S. N. Levine and J. O. Smith III. A sines+transients+noise audio representation for data compression and time/pitch scale
modications. (1998)
L. Daudet and B. Torr´esani. Hybrid representations for audiophonic signal encoding. (2002)
S. Molla and B. Torr´esani. An hybrid audio coding scheme using hidden Markov models of waveforms. (2005)
M. E. Davies and L. Daudet. Sparse audio representations using the MCLT. (2006)
2
Oscillatory (rhythmic) and Transient Components in EEG
Many measured signals have both an oscillatory and a non-oscillatory component.
0 1 2 3 4 5 6 7 8 9 10
−40
−20
0
20
40
EEG SIGNAL
TIME (SECONDS)
Rhythms of the EEG:
Delta 0 - 3 Hz
Theta 4 - 7 Hz
Alpha 8 - 12 Hz
Beta 12 - 30 Hz
Gamma 26 - 100 Hz
Transients in EEG due to:
1) unwanted measurement artifacts
2) non-rhythmic brain activity (spikes, spindles, and vertex waves)
3
Signal resonance and Q-factor
0 100 200 300 400 500
−1
0
1
PULSE 1
0 100 200 300 400 500
−1
0
1
PULSE 2
0 100 200 300 400 500
−1
0
1
PULSE 3
0 100 200 300 400 500
−1
0
1
PULSE 4
TIME (SAMPLES)
0 0.02 0.04 0.06 0.08 0.1 0.12
0
5
10
SPECTRUM
Q−factor = 1.15
0 0.02 0.04 0.06 0.08 0.1 0.12
0
20
40
Q−factor = 4.6
0 0.02 0.04 0.06 0.08 0.1 0.12
0
10
20
Q−factor = 1.15
0 0.02 0.04 0.06 0.08 0.1 0.12
0
50
100
Q−factor = 4.6
FREQUENCY (CYCLES/SAMPLE)
(a) Signals. (b) Spectra.
Figure 1: The resonance of an isolated pulse can be quantified by its Q-factor, defined as the ratio of its center frequency to its bandwidth. Pulses 1 and 3,
essentially a single cycle in duration, are low-resonance pulses. Pulses 2 and 4, whose oscillations are more sustained, are high-resonance pulses.
4
Resonance-based signal decomposition
0 100 200 300 400 500
−2
−1
0
1
2
TEST SIGNAL
0 100 200 300 400 500
−2
−1
0
1
2HIGH−RESONANCE COMPONENT
0 100 200 300 400 500
−2
−1
0
1
2LOW−RESONANCE COMPONENT
0 100 200 300 400 500
−2
−1
0
1
2RESIDUAL
TIME (SAMPLES)
0 100 200 300 400 500
−2
−1
0
1
2
TEST SIGNAL
0 100 200 300 400 500
−2
−1
0
1
2LOW−PASS FILTERED
0 100 200 300 400 500
−2
−1
0
1
2BAND−PASS FILTERED
0 100 200 300 400 500
−2
−1
0
1
2HIGH−PASS FILTERED
TIME (SAMPLES)
(a) Resonance-based decomposition. (b) Frequency-based filtering.
Figure 2: Resonance- and frequency-based filtering. (a) Decomposition of a test signal into high- and low-resonance components. The high-resonance signal
component is sparsely represented using a high Q-factor RADWT. Similarly, the low-resonance signal component is sparsely represented using a low Q-factor
RADWT. (b) Decomposition of a test signal into low, mid, and high frequency components using LTI discrete-time filters.
5
Resonance-based signal decomposition must be nonlinear
Figure 3: Resonance-based signal decomposition must be nonlinear: The signal in the bottom left panel is the sum of the signals above it; however, the low-resonance
component of a sum is not the sum of the low-resonance components. The same is true for the high-resonance component. Neither the low- nor high-resonance
components satisfy the superposition property.
6
Rational-dilation wavelet transform (RADWT)
Prior work on rational-dilation wavelet transforms addresses the critically-sampled case.
1. K. Nayebi, T. P. Barnwell III and M. J. T. Smith (1991)
2. P. Auscher (1992)
3. J. Kovacevic and M. Vetterli (1993)
4. T. Blu (1993, 1996, 1998)
5. A. Baussard, F. Nicolier and F. Truchetet (2004)
6. G. F. Choueiter and J. R. Glass (2007)
New RADWT (2009) gives a solution for the overcomplete case.
7
Rational-dilation wavelet transform (RADWT)
p H(ω)q
G(ω)s
q H(ω)p
+
sG(ω)
Figure 4: Analysis and synthesis filter banks for the implementation of the rational-dilation wavelet transform (RADWT). The dilation factor is q/p and the
redundancy is (s(1 p/q))1assuming iteration of the filter bank on its low-pass (upper) branch ad infinitum.
Y(ω) =
q1
X
k=0
Lk(ω)Xω+p k2π
q+
s1
X
k=0
Mk(ω)Xω+k2π
s,
where
Lk(ω) = 1
p q
p1
X
n=0
Hω
p+k2π
q+n2π
pHω
p+n2π
p,
Mk(ω) = 1
sGω+k2π
sG(ω).
There are no perfect reconstruction filters with rational-transfer functions unless the filter bank is or-
thonormal (p+ 1 = q=s).
=For overcomplete case, use filters with non-rational transfer functions for PR.
8
Perfect reconstruction can be attained by
H(ω) =
pq |ω| 11
sπ
p
0|ω| π
q, π
G(ω) =
0|ω| 11
sπ
s|ω| p
qπ, π
where transition-bands are chosen so as to satisfy 1
p q Hω
p
2+1
s|G(ω)|2= 1.
0π
¯
H(ω)
qq
p
s
(s1)π
s
q
G(ω)
ω
With q= 3,p= 2,s= 2:
0!"2 2!"3 !
0
"
H(ω)G(ω)2
!3/2
9
Low Q-factor RADWT
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
1.5
FREQUENCY
p = 2, q = 3, s = 1
FREQUENCY RESPONSES
0 50 100 150 200 250 300
−0.1
−0.05
0
0.05
0.1
0.15 WAVELET
TIME
Figure 5: Low Q-factor rational-dilation wavelet transforms (RADWT) with p= 2, q= 3, s= 1. The wavelet is approximately the Mexican hat function. The
RADWT is 3-times overcomplete.
10
High Q-factor RADWT
0 0.1 0.2 0.3 0.4 0.5
0
1
2
3
FREQUENCY
p = 5, q = 6, s = 2
FREQUENCY RESPONSES
0 50 100 150 200 250 300
−0.1
0
0.1
WAVELET
TIME
Figure 6: High Q-factor rational-dilation wavelet transforms (RADWT) with p= 5, q= 6, s= 2. The dilation factor is 1.2, much closer to 1 than the dyadic
wavelet transform. The RADWT is 3-times overcomplete.
11
Low Q-factor vs High Q-factor RADWT
0 100 200 300 400 500
9
8
7
6
5
4
3
2
1
1.89%
6.14%
11.72%
16.33%
22.82%
26.80%
13.92%
0.38%
0.00%
SUBBANDS (LOW−Q RADWT)
TIME (SAMPLES)
SUBBAND
P = 2, Q = 3, S = 1, Levels = 10 Dilation = 1.50, Redundancy = 3.00
0 100 200 300 400 500
17
16
15
14
13
12
11
10
9
8
7
2.43%
7.97%
3.52%
0.54%
7.30%
16.03%
4.72%
1.74%
20.12%
29.69%
5.75%
SUBBANDS (HIGH−Q RADWT)
TIME (SAMPLES)
SUBBAND
P = 5, Q = 6, S = 2, Levels = 20 Dilation = 1.20, Redundancy = 3.00
12
Low Q-factor vs High Q-factor RADWT after sparsification
0 100 200 300 400 500
9
8
7
6
5
4
3
2
1
0.00%
2.07%
8.78%
34.39%
0.79%
53.29%
0.68%
0.00%
0.00%
SUBBANDS (LOW−Q RADWT)
TIME (SAMPLES)
SUBBAND
P = 2, Q = 3, S = 1, Levels = 10 Dilation = 1.50, Redundancy = 3.00
0 100 200 300 400 500
17
16
15
14
13
12
11
10
9
8
7
0.01%
15.64%
0.11%
0.00%
0.02%
38.66%
0.02%
0.02%
6.26%
39.25%
0.00%
SUBBANDS (HIGH−Q RADWT)
TIME (SAMPLES)
SUBBAND
P = 5, Q = 6, S = 2, Levels = 20 Dilation = 1.20, Redundancy = 3.00
13
Low Q-factor vs High Q-factor RADWT
0 100 200 300 400 500
9
8
7
6
5
4
3
2
1
2.58%
6.31%
10.70%
16.11%
21.44%
22.50%
13.51%
4.87%
1.42%
SUBBANDS (LOW−Q RADWT)
TIME (SAMPLES)
SUBBAND
P = 2, Q = 3, S = 1, Levels = 10 Dilation = 1.50, Redundancy = 3.00
0 100 200 300 400 500
19
18
17
16
15
14
13
12
11
10
9
8
7
6
1.82%
2.75%
3.02%
3.82%
5.65%
6.82%
6.95%
8.72%
12.37%
13.98%
12.04%
8.44%
5.41%
3.07%
SUBBANDS (HIGH−Q RADWT)
TIME (SAMPLES)
SUBBAND
P = 5, Q = 6, S = 2, Levels = 20 Dilation = 1.20, Redundancy = 3.00
14
Low Q-factor vs High Q-factor RADWT after sparsification
0 100 200 300 400 500
9
8
7
6
5
4
3
2
1
0.00%
3.12%
10.99%
27.98%
10.60%
40.59%
3.53%
2.64%
0.23%
SUBBANDS (LOW−Q RADWT)
TIME (SAMPLES)
SUBBAND
P = 2, Q = 3, S = 1, Levels = 10 Dilation = 1.50, Redundancy = 3.00
0 100 200 300 400 500
19
18
17
16
15
14
13
12
11
10
9
8
7
6
2.43%
3.78%
2.35%
2.16%
4.12%
5.87%
6.25%
4.01%
18.71%
12.72%
9.69%
19.51%
1.36%
3.11%
SUBBANDS (HIGH−Q RADWT)
TIME (SAMPLES)
SUBBAND
P = 5, Q = 6, S = 2, Levels = 20 Dilation = 1.20, Redundancy = 3.00
15
Constant-Q vs Constant-BW
Constant-Q Constant-BW
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
0
0.5
1
FREQUENCY (CYCLES/SAMPLE)
FREQUENCY RESPONSES
0 50 100 150 200 250 300
−1
0
1
ANALYSIS FUNCTION
TIME (SAMPLES)
0 50 100 150 200 250 300
−1
0
1
ANALYSIS FUNCTION
TIME (SAMPLES)
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
0
0.5
1
FREQUENCY (CYCLES/SAMPLE)
FREQUENCY RESPONSES
0 50 100 150 200 250 300
−1
0
1
TIME (SAMPLES)
ANALYSIS FUNCTION
0 50 100 150 200 250 300
−1
0
1
TIME (SAMPLES)
ANALYSIS FUNCTION
fixed ‘resonance’ frequency-dependent ‘resonance’
frequency-dependent temporal duration fixed temporal duration
16
Rational-dilation wavelet transform (RADWT)
Summary:
1. Fully-discrete, modestly overcomplete
2. Exact perfect reconstruction (‘self-inverting’)
3. Adjustable Q-factor:
Can attain higher Q-factors than (or same low Q-factor of) the dyadic WT.
=Can achieve higher-frequency resolution needed for oscillatory signals.
4. Samples the time-frequency plane more densely in both time and frequency.
=Exactly invertible, fully-discrete approximation of the continuous WT.
5. FFT-based implementation
Reference:
Bayram, Selesnick. Frequency-domain design of overcomplete rational-dilation wavelet transforms. IEEE Trans. on Signal
Processing, 57, August 2009.
17
Morphological Component Analysis (MCA)
Given an observed signal
x=x1+x2,with x,x1,x2RN,
the goal of MCA is to estimate/determine x1and x2individually. Assuming that x1and x2can be
sparsely represented in bases (or frames) S1and S2respectively, they can be estimated by minimizing the
objective function,
J(w1,w2) = kxS1w1S2w2k2
2+λ1kw1k1+λ2kw2k1
with respect to w1and w2. Then MCA provides the estimates
ˆ
x1=S1w1
and
ˆ
x2=S2w2.
Reference:
Starck, Elad, Donoho. Image Decomposition via the Combination of Sparse Representations and a Variational Approach,
IEEE Trans. on Image Processing, Oct 2005.
18
Why not an `2-norm penalty?
If the `2-norm is used for the penalty term,
J(w1,w2) = kxS1w1S2w2k2
2+λ1kw1k2
2+λ2kw2k2
2,
then, using S1St
1=S2St
2=I, the minimizing w1and w2can be found in closed form:
w1=λ1
λ1+λ2+λ1λ2
St
1x
w2=λw
λ1+λ2+λ1λ2
St
2x
and the estimated components ˆ
x1=S1w1and ˆ
x2=S2w2are given by
ˆ
x1=λ1
λ1+λ2+λ1λ2
x
ˆ
x2=λ2
λ1+λ2+λ1λ2
x
Both ˆ
x1and ˆ
x2are just scaled versions of x.
=No separation at all!
19
MCA as a linear inverse problem
The objective function
J(w1,w2) = kxS1w1S2w2k2
2+λ1kw1k1+λ2kw2k1
can be written as
J(w) = kxHwk2
2+kλtwk1
where
H=hS1S2i,w=
w1
w2
.
An `1-regularized linear inverse problem . . .
Non-differentiable
Convex
=Use Iterative Soft Thresholding Algorithm (ISTA) or another algorithm to minimize J(w).
Other algorithms include FISTA, SALSA, TwIST, etc.
20
Split augmented Lagrangian shrinkage algorithm (SALSA)
SALSA is an algorithm for minimizing
J(w) = kxHwk2
2+λkwk1
SALSA is based on the minimization of
min
uf1(u) + f2(u) (1)
by the alternating split augmented Lagrangian algorithm:
u(k+1) = arg min
uf1(u) + µkuv(k)d(k)k2
2(2)
v(k+1) = arg min
vf2(v) + µku(k+1) vd(k)k2
2(3)
d(k+1) =d(k)u(k+1) +v(k+1) (4)
Reference:
Figueiredo, Bioucas-Dias, Afonso. Fast Frame-Based Image Deconvolution Using Variable Splitting and Constrained
Optimization. IEEE Workshop on Statistical Signal Processing, 2009.
21
SALSA
Applying SALSA to the MCA problem yields the iterative algorithm:
b(k)
1=St
1x+µ(w(k)
1+d(k)
1) (5)
b(k)
2=St
2x+µ(w(k)
2+d(k)
2) (6)
c(k)=S1b(k)
1+S2b(k)
2(7)
u(k+1)
1=1
µb(k)
11
µ(µ+ 2) St
1c(k)(8)
u(k+1)
2=1
µb(k)
21
µ(µ+ 2) St
2c(k)(9)
w(k+1)
1= softu(k+1)
1d(k)
1,λ1
2µ(10)
w(k+1)
2= softu(k+1)
2d(k)
2,λ2
2µ(11)
where soft(x, T )is the soft-threshold rule with threshold T,
soft(x, T ) = xmax(0,1T/|x|).
Note: no matrix inverses; only forward and inverse transforms.
22
0 10 20 30 40 50 60 70 80 90 100
5.5
6
6.5
7
7.5
8
8.5
9
9.5
ITERATION
OBJECTIVE FUNCTION ISTA
SALSA
Figure 7: Reduction of objective function during the first 100 iterations. SALSA converges faster than ISTA.
23
Constant bandwidth + Constant Q-factor
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
CONSTANT BANDWITH
FREQUENCY
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
CONSTANT Q−FACTOR
FREQUENCY
A constant bandwidth and a constant Q-factor decomposition can have high coherence due to some
analysis functions, from each decomposition, having similar frequency support. This can degrade the
results of MCA in principle.
24
Constant Q-factor: High Q-factor + Low Q-factor
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
CONSTANT Q−FACTOR (HIGH Q−FACTOR)
FREQUENCY
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
CONSTANT Q−FACTOR (LOW Q−FACTOR)
FREQUENCY
Two constant Q-factor decompositions with markedly different Q-factors will have low coherence because
no analysis functions from the two decompositions will have similar frequency support. This is beneficial
for the operation of MCA.
25
Small coherence between low and high Q-factor RADWTs
0f1
pQ1/f1Q1/f1
Ψ1(f)
f2
pQ2/f2Q2/f2
Ψ2(f)
f
Figure 8: For reliable resonance-based decomposition, the inner product between the low-Q and high-Q wavelets should be small for all dilations and translations.
The computation of the maximum inner product is simplified by assuming the wavelets are ideal band-pass functions and expressing the inner product in the
frequency domain.
The inner products can be defined in the frequency domain,
ρ(f1, f2) := ZΨ1(f2(f)df,
as a function of their center frequencies (equivalently, dilation).
The maximum value of the inner product, ρ(f1, f2), occurs when f2=f1(2 + 1/Q1)/(2 + 1/Q2)and is
given by
ρmax =sQ1+ 1/2
Q2+ 1/2, Q2> Q1.(12)
26
Constant Bandwidth: Narrow-band + Wide-band
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
CONSTANT BANDWITH (NARROW−BAND)
FREQUENCY
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
CONSTANT BANDWITH (WIDE−BAND)
FREQUENCY
Two constant bandwidth decompositions with markedly different bandwidths will also have low coherence
and are therefore also suitable transform for MCA-based signal decomposition. This gives a bandwidth-
based decomposition, rather than a resonance-based decomposition.
27
Example: Resonance-selective nonlinear band-pass filtering
0 100 200 300 400
−2
−1
0
1
2
(a) TEST SIGNAL
TIME (SAMPLES)
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
FREQUENCY (CYCLES/SAMPLE)
(b) TWO BAND−PASS FILTERS
BAND−PASS FILTER 1
BAND−PASS FILTER 2
0 100 200 300 400
−2
−1
0
1
2
(c) OUTPUT OF BAND−PASS FILTER 1
TIME (SAMPLES)
0 100 200 300 400
−2
−1
0
1
2(d) OUTPUT OF BAND−PASS FILTER 2
TIME (SAMPLES)
Figure 9: LTI band-pass filtering. The test signal (a) consists of a sinusoidal pulse of frequency 0.1 cycles/sample and a transient. Band-pass filters 1 and 2 in
(b) are tuned to the frequencies 0.07 and 0.10 cycles/second respectively. The output signals, obtained by filtering the test signal with each of the two band-pass
filters, are shown in (c) and (d). The output of band-pass filter 1, illustrated in (c), contains oscillations due to the transient in the test signal. Moreover, the
transient oscillations in (c) have a frequency of 0.07 Hz even though the test signal (a) contains no sustained oscillatory behavior at this frequency.
28
0 100 200 300 400
−2
−1
0
1
2
(a) HIGH−RESONANCE COMPONENT
TIME (SAMPLES)
0 100 200 300 400
−2
−1
0
1
2(b) LOW−RESONANCE COMPONENT
TIME (SAMPLES)
0 100 200 300 400
−2
−1
0
1
2
(c) OUTPUT OF BAND−PASS FILTER 1
TIME (SAMPLES)
0 100 200 300 400
−2
−1
0
1
2(d) OUTPUT OF BAND−PASS FILTER 2
TIME (SAMPLES)
Figure 10: Resonance-based decomposition and band-pass filtering. When resonance-based analysis method is applied to the test signal in Fig. 9a, it yields the
high- and low-resonance components illustrated in (a) and (b). The output signals, obtained by filtering the high-resonance component (a) with each of the two
band-pass filters shown in Fig. 9b, are illustrated in (c) and (d). The transient oscillations in (c) are substantially reduced compared to Fig. 9c.
29
Example: Resonance-based decomposition of speech
0 0.05 0.1 0.15
−0.4
−0.2
0
0.2
0.4
(a) SPEECH SIGNAL [Fs = 16,000 SAMPLES/SECOND]
0 0.05 0.1 0.15
−0.4
−0.2
0
0.2
0.4 (b) HIGH−RESONANCE COMPONENT
0 0.05 0.1 0.15
−0.4
−0.2
0
0.2
0.4 (c) LOW−RESONANCE COMPONENT
TIME (SECONDS)
Figure 11: Decomposition of a speech signal (“I’m”) into high- and low-resonance components. The high-resonance component (b) contains the sustained oscillations
present in the speech signal, while the low-resonance component (c) contains non-oscillatory transients. (The residual is not shown.)
30
0 500 1000 1500 2000 2500 3000
0
0.005
0.01
(a) ORIGINAL SPEECH
0 500 1000 1500 2000 2500 3000
0
0.005
0.01 (b) HIGH−RESONANCE COMPONENT
0 500 1000 1500 2000 2500 3000
0
0.005
0.01 (c) LOW−RESONANCE COMPONENT
FREQUENCY (Hz)
Figure 12: Frequency spectra of the speech signal in Fig. 11 and of the extracted high- and low-resonance components. The spectra are computed using the 50
msec segment from 0.05 to 0.10 seconds. The energy of each resonance component is widely distributed in frequency and their frequency-spectra overlap.
31
0 0.05 0.1 0.15
−0.2
0
0.2
(a) RECONSTRUCTION FROM SUBBANDS 9, 10
0 0.05 0.1 0.15
−0.2
0
0.2
(b) RECONSTRUCTION FROM SUBBANDS 18, 19
0 0.05 0.1 0.15
−0.2
0
0.2
(c) RECONSTRUCTION FROM SUBBANDS 21−24
0 0.05 0.1 0.15
−0.2
0
0.2
(D) SUM OF ABOVE 3 SIGNALS
TIME (SECONDS)
Figure 13: Frequency decomposition of high-resonance component in Fig. 11. Reconstructing the high-resonance component from a few subbands of the high
Q-factor RADWT at a time, yields an efficient AM/FM decomposition.
32
0 500 1000 1500 2000 2500
13
12
11
10
9
8
7
6
5
4
3
2
1
RECONSTRUCTION FROM ONE SINGLE SUBBAND (LOW Q!FACTOR RADWT)
n
SUBBAND
33
0 500 1000 1500 2000 2500
26
25
24
23
22
21
20
19
18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
RECONSTRUCTION FROM ONE SINGLE SUBBAND (HIGH Q!FACTOR RADWT)
n
SUBBAND
34
Conclusion: Resonance-based signal decomposition
Low Q-factor RADWT used for sparse representation of the transient component.
High Q-factor RADWT used for sparse representation of the oscillatory (rhythmic) component.
Morphological component analysis (MCA) used to separate the two signal components.
Oscillatory component not necessarily high-pass may contain both low and high frequencies.
Transient component not necessarily a low-pass signal may contains sharp bumps and jumps.
35