SAMV (algorithm)

SAMV (iterative sparse asymptotic minimum variance^{[1] [2]} ) is a parameter-free superresolution algorithm for the linear inverse problem in spectral estimation, direction-of-arrival (DOA) estimation and tomographic reconstruction with applications in signal processing, medical imaging and remote sensing. The name was coined in 2013^[1] to emphasize its basis on the asymptotically minimum variance (AMV) criterion. It is a powerful tool for the recovery of both the amplitude and frequency characteristics of multiple highly correlated sources in challenging environments (e.g., limited number of snapshots and low signal-to-noise ratio). Applications include synthetic-aperture radar,^{[2] [3]} computed tomography scan, and magnetic resonance imaging (MRI).

The formulation of the SAMV algorithm is given as an inverse problem in the context of DOA estimation. Suppose an ${\displaystyle M}${\displaystyle M}-element uniform linear array (ULA) receive ${\displaystyle K}${\displaystyle K} narrow band signals emitted from sources located at locations ${\displaystyle \mathbf {\theta } =\{\theta _{a},\ldots ,\theta _{K}\}}${\displaystyle \mathbf {\theta } =\{\theta _{a},\ldots ,\theta _{K}\}}, respectively. The sensors in the ULA accumulates ${\displaystyle N}${\displaystyle N} snapshots over a specific time. The ${\displaystyle M\times 1}${\displaystyle M\times 1} dimensional snapshot vectors are

${\displaystyle \mathbf {y} (n)=\mathbf {A} \mathbf {x} (n)+\mathbf {e} (n),n=1,\ldots ,N}${\displaystyle \mathbf {y} (n)=\mathbf {A} \mathbf {x} (n)+\mathbf {e} (n),n=1,\ldots ,N}

where ${\displaystyle \mathbf {A} =[\mathbf {a} (\theta _{1}),\ldots ,\mathbf {a} (\theta _{K})]}${\displaystyle \mathbf {A} =[\mathbf {a} (\theta _{1}),\ldots ,\mathbf {a} (\theta _{K})]} is the steering matrix, ${\displaystyle {\bf {x}}(n)=[{\bf {x}}_{1}(n),\ldots ,{\bf {x}}_{K}(n)]^{T}}${\displaystyle {\bf {x}}(n)=[{\bf {x}}_{1}(n),\ldots ,{\bf {x}}_{K}(n)]^{T}} contains the source waveforms, and ${\displaystyle {\bf {e}}(n)}${\displaystyle {\bf {e}}(n)} is the noise term. Assume that ${\displaystyle \mathbf {E} \left({\bf {e}}(n){\bf {e}}^{H}({\bar {n}})\right)=\sigma {\bf {I}}_{M}\delta _{n,{\bar {n}}}}${\displaystyle \mathbf {E} \left({\bf {e}}(n){\bf {e}}^{H}({\bar {n}})\right)=\sigma {\bf {I}}_{M}\delta _{n,{\bar {n}}}}, where ${\displaystyle \delta _{n,{\bar {n}}}}${\displaystyle \delta _{n,{\bar {n}}}} is the Dirac delta and it equals to 1 only if ${\displaystyle n={\bar {n}}}${\displaystyle n={\bar {n}}} and 0 otherwise. Also assume that ${\displaystyle {\bf {e}}(n)}${\displaystyle {\bf {e}}(n)} and ${\displaystyle {\bf {x}}(n)}${\displaystyle {\bf {x}}(n)} are independent, and that ${\displaystyle \mathbf {E} \left({\bf {x}}(n){\bf {x}}^{H}({\bar {n}})\right)={\bf {P}}\delta _{n,{\bar {n}}}}${\displaystyle \mathbf {E} \left({\bf {x}}(n){\bf {x}}^{H}({\bar {n}})\right)={\bf {P}}\delta _{n,{\bar {n}}}}, where ${\displaystyle {\bf {P}}=\operatorname {Diag} ({p_{1},\ldots ,p_{K}})}${\displaystyle {\bf {P}}=\operatorname {Diag} ({p_{1},\ldots ,p_{K}})}. Let ${\displaystyle {\bf {p}}}${\displaystyle {\bf {p}}} be a vector containing the unknown signal powers and noise variance, ${\displaystyle {\bf {p}}=[p_{1},\ldots ,p_{K},\sigma ]^{T}}${\displaystyle {\bf {p}}=[p_{1},\ldots ,p_{K},\sigma ]^{T}}.

The covariance matrix of ${\displaystyle {\bf {y}}(n)}${\displaystyle {\bf {y}}(n)} that contains all information about ${\displaystyle {\boldsymbol {\bf {p}}}}${\displaystyle {\boldsymbol {\bf {p}}}} is

${\displaystyle {\bf {R}}={\bf {A}}{\bf {P}}{\bf {A}}^{H}+\sigma {\bf {I}}.}${\displaystyle {\bf {R}}={\bf {A}}{\bf {P}}{\bf {A}}^{H}+\sigma {\bf {I}}.}

This covariance matrix can be traditionally estimated by the sample covariance matrix ${\displaystyle {\bf {R}}_{N}={\bf {Y}}{\bf {Y}}^{H}/N}${\displaystyle {\bf {R}}_{N}={\bf {Y}}{\bf {Y}}^{H}/N} where ${\displaystyle {\bf {Y}}=[{\bf {y}}(1),\ldots ,{\bf {y}}(N)]}${\displaystyle {\bf {Y}}=[{\bf {y}}(1),\ldots ,{\bf {y}}(N)]}. After applying the vectorization operator to the matrix ${\displaystyle {\bf {R}}}${\displaystyle {\bf {R}}}, the obtained vector ${\displaystyle {\bf {r}}({\boldsymbol {\bf {p}}})=\operatorname {vec} ({\bf {R}})}${\displaystyle {\bf {r}}({\boldsymbol {\bf {p}}})=\operatorname {vec} ({\bf {R}})} is linearly related to the unknown parameter ${\displaystyle {\boldsymbol {\bf {p}}}}${\displaystyle {\boldsymbol {\bf {p}}}} as

${\displaystyle {\bf {r}}({\boldsymbol {\bf {p}}})=\operatorname {vec} ({\bf {R}})={\bf {S}}{\boldsymbol {\bf {p}}}}${\displaystyle {\bf {r}}({\boldsymbol {\bf {p}}})=\operatorname {vec} ({\bf {R}})={\bf {S}}{\boldsymbol {\bf {p}}}},

where ${\displaystyle {\bf {S}}=[{\bf {S}}_{1},{\bar {\bf {a}}}_{K+1}]}${\displaystyle {\bf {S}}=[{\bf {S}}_{1},{\bar {\bf {a}}}_{K+1}]}, ${\displaystyle {\bf {S}}_{1}=[{\bar {\bf {a}}}_{1},\ldots ,{\bar {\bf {a}}}_{K}]}${\displaystyle {\bf {S}}_{1}=[{\bar {\bf {a}}}_{1},\ldots ,{\bar {\bf {a}}}_{K}]}, ${\displaystyle {\bar {\bf {a}}}_{k}={\bf {a}}_{k}^{*}\otimes {\bf {a}}_{k}}${\displaystyle {\bar {\bf {a}}}_{k}={\bf {a}}_{k}^{*}\otimes {\bf {a}}_{k}}, ${\displaystyle k=1,\ldots ,K}${\displaystyle k=1,\ldots ,K}, and let ${\displaystyle {\bar {\bf {a}}}_{K+1}=\operatorname {vec} ({\bf {I}})}${\displaystyle {\bar {\bf {a}}}_{K+1}=\operatorname {vec} ({\bf {I}})} where ${\displaystyle \otimes }${\displaystyle \otimes } is the Kronecker product.

To estimate the parameter ${\displaystyle {\boldsymbol {\bf {p}}}}${\displaystyle {\boldsymbol {\bf {p}}}} from the statistic ${\displaystyle {\bf {r}}_{N}}${\displaystyle {\bf {r}}_{N}}, we develop a series of iterative SAMV approaches based on the asymptotically minimum variance criterion. From ^[1] the covariance matrix ${\displaystyle \operatorname {Cov} _{\boldsymbol {p}}^{\operatorname {Alg} }}${\displaystyle \operatorname {Cov} _{\boldsymbol {p}}^{\operatorname {Alg} }} of an arbitrary consistent estimator of ${\displaystyle {\boldsymbol {p}}}${\displaystyle {\boldsymbol {p}}} based on the second-order statistic ${\displaystyle {\bf {r}}_{N}}${\displaystyle {\bf {r}}_{N}} is bounded by the real symmetric positive definite matrix

${\displaystyle \operatorname {Cov} _{\boldsymbol {p}}^{\operatorname {Alg} }\geq [{\bf {S}}_{d}^{H}{\bf {C}}_{r}^{-1}{\bf {S}}_{d}]^{-1},}${\displaystyle \operatorname {Cov} _{\boldsymbol {p}}^{\operatorname {Alg} }\geq [{\bf {S}}_{d}^{H}{\bf {C}}_{r}^{-1}{\bf {S}}_{d}]^{-1},}

where ${\displaystyle {\bf {S}}_{d}={\rm {d}}{\bf {r}}({\boldsymbol {p}})/{\rm {d}}{\boldsymbol {p}}}${\displaystyle {\bf {S}}_{d}={\rm {d}}{\bf {r}}({\boldsymbol {p}})/{\rm {d}}{\boldsymbol {p}}}. In addition, this lower bound is attained by the covariance matrix of the asymptotic distribution of ${\displaystyle {\hat {\bf {p}}}}${\displaystyle {\hat {\bf {p}}}} obtained by minimizing

${\displaystyle {\hat {\boldsymbol {p}}}=\arg \min _{\boldsymbol {p}}f({\boldsymbol {p}}),}${\displaystyle {\hat {\boldsymbol {p}}}=\arg \min _{\boldsymbol {p}}f({\boldsymbol {p}}),}

where ${\displaystyle f({\boldsymbol {p}})=[{\bf {r}}_{N}-{\bf {r}}({\boldsymbol {p}})]^{H}{\bf {C}}_{r}^{-1}[{\bf {r}}_{N}-{\bf {r}}({\boldsymbol {p}})].}${\displaystyle f({\boldsymbol {p}})=[{\bf {r}}_{N}-{\bf {r}}({\boldsymbol {p}})]^{H}{\bf {C}}_{r}^{-1}[{\bf {r}}_{N}-{\bf {r}}({\boldsymbol {p}})].}

Therefore, the estimate of ${\displaystyle {\boldsymbol {\bf {p}}}}${\displaystyle {\boldsymbol {\bf {p}}}} can be obtained iteratively.

The ${\displaystyle \{{\hat {p}}_{k}\}_{k=1}^{K}}${\displaystyle \{{\hat {p}}_{k}\}_{k=1}^{K}} and ${\displaystyle {\hat {\sigma }}}${\displaystyle {\hat {\sigma }}} that minimize ${\displaystyle f({\boldsymbol {p}})}${\displaystyle f({\boldsymbol {p}})} can be computed as follows. Assume ${\displaystyle {\hat {p}}_{k}^{(i)}}${\displaystyle {\hat {p}}_{k}^{(i)}} and ${\displaystyle {\hat {\sigma }}^{(i)}}${\displaystyle {\hat {\sigma }}^{(i)}} have been approximated to a certain degree in the ${\displaystyle i}${\displaystyle i}th iteration, they can be refined at the ${\displaystyle (i+1)}${\displaystyle (i+1)}th iteration by

${\displaystyle {\hat {p}}_{k}^{(i+1)}={\frac {{\bf {a}}_{k}^{H}{\bf {R}}^{-1{(i)}}{\bf {R}}_{N}{\bf {R}}^{-1{(i)}}{\bf {a}}_{k}}{({\bf {a}}_{k}^{H}{\bf {R}}^{-1{(i)}}{\bf {a}}_{k})^{2}}}+{\hat {p}}_{k}^{(i)}-{\frac {1}{{\bf {a}}_{k}^{H}{\bf {R}}^{-1{(i)}}{\bf {a}}_{k}}},\quad k=1,\ldots ,K}${\displaystyle {\hat {p}}_{k}^{(i+1)}={\frac {{\bf {a}}_{k}^{H}{\bf {R}}^{-1{(i)}}{\bf {R}}_{N}{\bf {R}}^{-1{(i)}}{\bf {a}}_{k}}{({\bf {a}}_{k}^{H}{\bf {R}}^{-1{(i)}}{\bf {a}}_{k})^{2}}}+{\hat {p}}_{k}^{(i)}-{\frac {1}{{\bf {a}}_{k}^{H}{\bf {R}}^{-1{(i)}}{\bf {a}}_{k}}},\quad k=1,\ldots ,K}

${\displaystyle {\hat {\sigma }}^{(i+1)}=\left(\operatorname {Tr} ({\bf {R}}^{-2^{(i)}}{\bf {R}}_{N})+{\hat {\sigma }}^{(i)}\operatorname {Tr} ({\bf {R}}^{-2^{(i)}})-\operatorname {Tr} ({\bf {R}}^{-1^{(i)}})\right)/{\operatorname {Tr} {({\bf {R}}^{-2^{(i)}})}},}${\displaystyle {\hat {\sigma }}^{(i+1)}=\left(\operatorname {Tr} ({\bf {R}}^{-2^{(i)}}{\bf {R}}_{N})+{\hat {\sigma }}^{(i)}\operatorname {Tr} ({\bf {R}}^{-2^{(i)}})-\operatorname {Tr} ({\bf {R}}^{-1^{(i)}})\right)/{\operatorname {Tr} {({\bf {R}}^{-2^{(i)}})}},}

where the estimate of ${\displaystyle {\bf {R}}}${\displaystyle {\bf {R}}} at the ${\displaystyle i}${\displaystyle i}th iteration is given by ${\displaystyle {\bf {R}}^{(i)}={\bf {A}}{\bf {P}}^{(i)}{\bf {A}}^{H}+{\hat {\sigma }}^{(i)}{\bf {I}}}${\displaystyle {\bf {R}}^{(i)}={\bf {A}}{\bf {P}}^{(i)}{\bf {A}}^{H}+{\hat {\sigma }}^{(i)}{\bf {I}}} with ${\displaystyle {\bf {P}}^{(i)}=\operatorname {Diag} ({\hat {p}}_{1}^{(i)},\ldots ,{\hat {p}}_{K}^{(i)})}${\displaystyle {\bf {P}}^{(i)}=\operatorname {Diag} ({\hat {p}}_{1}^{(i)},\ldots ,{\hat {p}}_{K}^{(i)})}.

The resolution of most compressed sensing based source localization techniques is limited by the fineness of the direction grid that covers the location parameter space.^[4] In the sparse signal recovery model, the sparsity of the truth signal ${\displaystyle \mathbf {x} (n)}${\displaystyle \mathbf {x} (n)} is dependent on the distance between the adjacent element in the overcomplete dictionary ${\displaystyle {\bf {A}}}${\displaystyle {\bf {A}}}, therefore, the difficulty of choosing the optimum overcomplete dictionary arises. The computational complexity is directly proportional to the fineness of the direction grid, a highly dense grid is not computational practical. To overcome this resolution limitation imposed by the grid, the grid-free SAMV-SML (iterative Sparse Asymptotic Minimum Variance - Stochastic Maximum Likelihood) is proposed,^[1] which refine the location estimates ${\displaystyle {\boldsymbol {\bf {\theta }}}=(\theta _{1},\ldots ,\theta _{K})^{T}}${\displaystyle {\boldsymbol {\bf {\theta }}}=(\theta _{1},\ldots ,\theta _{K})^{T}} by iteratively minimizing a stochastic maximum likelihood cost function with respect to a single scalar parameter ${\displaystyle \theta _{k}}${\displaystyle \theta _{k}}.

A typical application with the SAMV algorithm in SISO radar/sonar range-Doppler imaging problem. This imaging problem is a single-snapshot application, and algorithms compatible with single-snapshot estimation are included, i.e., matched filter (MF, similar to the periodogram or backprojection, which is often efficiently implemented as fast Fourier transform (FFT)), IAA,^[5] and a variant of the SAMV algorithm (SAMV-0). The simulation conditions are identical to:^[5] A ${\displaystyle 30}${\displaystyle 30}-element polyphase pulse compression P3 code is employed as the transmitted pulse, and a total of nine moving targets are simulated. Of all the moving targets, three are of ${\displaystyle 5}${\displaystyle 5} dB power and the rest six are of ${\displaystyle 25}${\displaystyle 25} dB power. The received signals are assumed to be contaminated with uniform white Gaussian noise of ${\displaystyle 0}${\displaystyle 0} dB power.

The matched filter detection result suffers from severe smearing and leakage effects both in the Doppler and range domain, hence it is impossible to distinguish the ${\displaystyle 5}${\displaystyle 5} dB targets. On contrary, the IAA algorithm offers enhanced imaging results with observable target range estimates and Doppler frequencies. The SAMV-0 approach provides highly sparse result and eliminates the smearing effects completely, but it misses the weak ${\displaystyle 5}${\displaystyle 5} dB targets.

An open source MATLAB implementation of SAMV algorithm could be downloaded here.

• Array processing – Area of research in signal processing
• Matched filter – Filters used in signal processing that are optimal in some sense
• Periodogram – Estimate of the spectral density of a signal
• Filtered backprojection – Integral transform in mathematics (Radon transform)
• MUltiple SIgnal Classification – Algorithm used for frequency estimation and radio direction finding (MUSIC), a popular parametric superresolution method
• Pulse-Doppler radar – Type of radar system
• Super-resolution imaging – Any technique to improve resolution of an imaging system beyond conventional limits
• Compressed sensing – Signal processing technique
• Inverse problem – Process of calculating the causal factors that produced a set of observations
• Tomographic reconstruction – Estimate object properties from a finite number of projections

• Signal estimation
• Fourier analysis
• Frequency-domain analysis
• Trigonometry
• Wave mechanics
• Medical imaging
• Inverse problems
• Multidimensional signal processing
• Signal processing
• Tomography

• Articles with short description
• Short description is different from Wikidata