Image Processing Reference
In-Depth Information
where U is the spectral support matrix given that depends on { f 1 , f 2 ,… f K } through
U =  { x ( f 1 ), x ( f 2 ),. . . x ( f K )} (12)
Note that if there is no noise and if all frequencies are known exactly, eq. (11) can be
verified by substituting y = U A ( U ), which is equivalent to eq. (5), on the right hand side
of eq. (11).
Finally, starting from estimates of the frequencies and amplitudes from OMP as described
above, apply weighted NLS to get better values. This is done by finding the frequency or set
of frequencies f = { f 1 , f 2 ,… f k } that minimize the functional R ( f ) given by
R ( f ) = |{ A [ U ( f )] U ( f ) - y } H W { A [ U ( f )] U ( f ) - y }|, (13)
which is the same as the weighted least squares estimator given by eq. (8) with the
substitution of A [ U(f) ] defined by eq. (11) for the amplitude vector and U ( f ) defined by eq.
(12) for the mixed sinusoids (see the analogous equations in Stoica and Moses, 1997, eqs.
4.3.7 and 4.3.8). The product A [ U(f) ] U ( f ) in eq. (13) is the same as Φ Z in eq. (8).
3.2 Algorithm description
As described in Table 1, the first step in the first iteration of the Do loop is estimation of the
first frequency in the spectral support of the signal X s . This is given by the frequency of the
sinusoid whose image after multiplication by Φ has the maximum correlation with the
observation vector y (see, for example, Tropp and Gilbert, 2007 Algorithm 3, step 2). Here
we use the equivalent form, the argmax of G ( f , y ) with respect to f to obtain the first estimate
of the frequency of the first sinusoid f 1 in eq. (4). At this point previous implementations of
discrete OMP use the amplitude estimator eq. (11) to estimate the amplitude of the first
sinusoid a 1 = A [ Φ x ( f 1 )], multiply this amplitude estimate times x ( f 1 ) , given by eq. (9), and
by the measurement or mixing matrix Φ and subtract this vector from the measurement
vector y to form the first residual r 1 .
In our algorithm, we proceed differently by improving the precision of the frequency
estimates using NLS before finding the amplitude estimate. We take the frequency f 1 from
the argmax of G ( f , y ) evaluated on a discrete set of frequencies and use that as the starting
value to solve the NLS problem given by eq. (13). We have used several methods and
several different software packages to solve the NLS problem. A simple decimation routine
[i.e., tabulating R ( f ) from f 1 - f to f 1 + f ( f is the over-complete grid spacing) in 10 steps,
finding the argmin, decreasing  f by a factor of 10, tabulating and finding the argmin of R ( f )
again until the specified precision is reached] works well but is not very efficient. Powell's
method in Python ( “scipy.optimize.fmin_powell”) and one of the Newton methods, the
PrincipalAxis method, and the Conjugate Gradient method in Mathematica 's minimizer
“FindMinimum” all work and take less time than the decimation routine. A detailed
investigation of minimizers for the NLS step in our version of OMP is beyond the scope of
this chapter. The oversampling N f required for our method and that required for
conventional OMP are nearly identical as discussed below in Section 6.
Given the better value of f 1 , we compute a 1 from eq. (11) and a new value of the residual r
with the NLS estimate of the first signal removed from y as in OMP
r 1 = y - A ( U 1 ) U 1 = y - a 1 Φ x ( f 1 ).
(14)
Search WWH ::




Custom Search