Abstract:
A scale adaptive filtering scheme is developed for underspread channels based on a model of the linear time varying channel operator as a process in scale. Recursions serve the purpose of adding detail to the filter estimate until a suitable measure of fidelity and complexity is achieved.

Description:
BACKGROUND OF THE INVENTION 
   1. Field of the Invention 
   The present invention relates generally to processing methods for adaptive filtering on a source signal for transmission, and more particularly an adaptive filtering method to identify the relationship between a source signal and a received signal. 
   2. Description of the Related Art 
   It will be readily appreciated by one skilled in the art that the use of adaptive filtering operations may be required in a variety of different signal processing applications. One example of such an area is that of acoustics and specifically sonar signaling. Adaptive filtering systems for time varying operators have broad applications. For example, echo cancellation, active and passive sonar, and equalization of communication signals employ adaptive filtering modeling regimes for extracting information from signals. These modeling regimes are often limited by available memory, computational resources, and acceptable processing delays. 
   Multiresolution models of time varying operators have been promoted with the idea of exploiting the sparsity or “economy” of representation attendant in them. This is a useful concept since many operators that are not shift invariant are often succinctly described in wavelet bases, the basic building block of multiresolution models. Conventional approaches have given a general framework for formulating the multiresolution filter structures and provided a fast least mean square-like estimators with an assumed maximal scale of representation. However, this assumption is often overly optimistic. In which case the maximal scale of representation or the most economical wavelet bases can be jointly estimated. 
   With regards to adaptive filtering, multiresolution models naturally allow estimation of the response at a given time instant based on given data. In addition, the model circumvents the need for locally stationary assumptions of time recursive algorithms. Accordingly, a larger amount of information is available for the estimation problem at each time instant. Unlike in-time estimation algorithms that are by their very nature causal or near causal, dependencies across time in the forward looking direction are not accounted for with classic in-time adaptive filtering algorithms. In-time estimation is a powerful paradigm due to its computational efficiencies and memory requirements and it continues to find broad appeal and enjoys superior performance in a variety of applications. Nevertheless in areas where fading and multipath delay do not conform to the wide sense stationary (WSS) assumption the multiresolution model is a viable solution. Other practical considerations related to using multiresolution modeling include signal processing and communication schemes involving strategies based on finite duration signaling. For example, mobile radio and underwater acoustic communications data employ a sequence of finite duration packets to transmit information. Similarly for underwater target localization by active sonar as well as radar applications source signals employ time localized “pings.” 
   Accordingly, a need exists for a scale adaptive filtering method that is able to match a received signal based on information from a source signal by estimating a variety of parameters. A channel operator is built up in scale, and the channel operator employs time delays and frequency spreads. For this purpose, each additional incrementation of the Doppler spread is hypothesized and then tested. 
   BRIEF SUMMARY OF THE INVENTION 
   An aspect of the present invention provides a scale adaptive filtering method to produce a match of the received signal based on the source signal, by estimating the channel operator at each scale along with nuisance parameters, for example. 
   This and other aspects are substantially achieved by providing a system and method for scale adaptive filtering to identify a relationship between a source signal and an actual received signal. The method includes calculating a scale wherein the scale includes a time delay spread and a frequency spread, estimating at least one wavelet coefficient, determining an estimated received signal at the scale using the wavelet coefficient, determining when the estimated received signal is sufficiently similar with respect to the actual received signal, and recalculating the scale by a factor when the estimated received signal is not sufficiently similar with respect to the actual received signal. 
   The above and still further aspects, features and advantages of the present invention will become apparent upon consideration of the following definitions, descriptions and descriptive figures of specific embodiments thereof, wherein like reference numerals in the various figures are utilized to designate like components. While these descriptions go into specific details of the invention, it should be understood that variations may and do exist and would be apparent to those skilled in the art based on the descriptions herein. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
       FIG. 1  is a flow chart representing a scale adaptive filtering method in accordance with an aspect of the present invention. 
   

   DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
   While the invention will be described in detail with reference to specific embodiments thereof, it will be apparent to one skilled in the art that various changes and modifications can be made therein without departing from the spirit and scope thereof. Accordingly, it is intended that the present invention covers modifications and variations of this invention provided they come within the scope of the appended claims and their equivalents. 
   Scale adaptive filtering is a general estimation method employed to estimate the relation between a received signal and a source signal in a variety of transmission environments. As is well known in the art, signaling information can be lost or distorted between a source signal and the received the signal. 
   The present invention is a novel method of applying transform function analysis to the received signal at a specific scale thus forming an estimated received signal. This estimated received signal is then compared to the actual received signal. If the estimated received signal is sufficiently complex to model the actual received signal, then the transform analysis has been successful. However, if the estimated signal is not sufficiently complex to model the actual received signal, then the scale is doubled and the transform analysis is applied again. This process is repeated, iteratively, until the estimated received signal is sufficiently complex to model the actual received signal. Sufficiency is achieved when a suitable measure of error is minimized without increasing the complexity of the channel operator. The complexity is limited such that the noise is not included in the modeled signal. In addition, achieving sufficiently similar signals is also determined by employing the Empirical Bayes estimator (a general statistical method), wherein the estimator provides a particular cost function. 
   The flow chart  10  illustrates the operation of an exemplary method of scale adaptive filtering in accordance with the present invention. Block  15  is representative of the start block. Block  20  includes generating a source signal from a source such as a transmitter on a submarine or other underwater vehicle, and receiving an actual received signal at a receiver. Block  25  is descriptive of the step including calculating a minimal scale, such that the scale includes time delay information and frequency spread information. 
   Block  30  is descriptive of the step of estimating a wavelet coefficient to be used in a transform analysis to determine an estimated received signal. In wavelet analysis it is helpful to view a signal in terms of a fully scalable modulated window. The window is shifted along the signal and for every position the spectrum is calculated. This process is repeated many times with a slightly shorter (or longer) window for every iteration. In the end the result is a collection of time-frequency representations of the signal. In the case of wavelets, we normally do not speak about time-frequency representations but about time-scale representations, scale being a means of quantifying the time variation in the channel operator. 
   Block  35  is descriptive of the step of determining an estimated received signal at the minimal scale, initially chosen in Block  20  and employing the wavelet coefficient of Block  30 . 
   Block  40  is descriptive of the step of determining if the estimated received signal is sufficiently similar to the actual received signal. In other words, a sufficiently similar signal is achieved when a suitable measure of error is minimized without increasing the complexity of the channel operator. This takes into account the amount of noise in the signal information. The noise or nuisance parameters in the actual received signal should not be modeled in the estimated received signal. If the estimated received signal is not sufficiently similar, then the method proceeds to Block  45  and the scale is increased by doubling the minimal scale and recalculating the wavelet coefficient by proceeding through Blocks  30 . Block  40  checks the estimated received signal against the actual received signal. If the two signals are sufficiently similar, according to the parameters outlined above, the method proceeds to Block  50  and Ends. 
   The flow chart  10  is an overview of the method described by the present invention. The following is an exemplary embodiment of this system and method. Specifically, when the maximal scale/frequency of representation of the most economical wavelet coefficients can not be assumed they are jointly estimated. Accordingly, the present invention exploits sparsity of representation under a priori uncertainty regarding Doppler spread and degree of sparsity of the operator. 
   Overall, the present invention employs three estimators. First, the wavelet coefficients of the channel response based on the flexible and adaptable Gaussian mixture model. Second, the estimation of the_minimal_scale of the multiresolution space in which the operator resides. Finally, an estimation of the nuisance parameters associated with the mixture model. 
   It is shown that the estimation of the maximal scale synonymous with the Doppler spread, is akin to model order selection and estimators for this parameter are provided. These channel estimators are recursive in scale each relying on the preceeding lower scale estimate as a starting point for the next higher resolution estimate. In this framework the channel is gradually built up in scale, rather than time, with each addition of detail a greater Doppler spread (decreased coherence time) is hypothesized and is tested against the data in a model selection framework. Thus multiple resolutions of the operator are provided in sequence with finer scale details descending from the lower scale estimates starting with the linear time invariant (LTI). The algorithm as presented is limited to filter operators that are underspread (i.e. Bτ max &lt;1 where τ max  is the delay spread of the filter, B the maximal Doppler bandwidth associated with any of the scatterers, this requirement excludes the underdetermined parameter estimation case). 
   The following provides preliminary notation to address the time varying channel as well as an introduction to the empirical Bayes approach. We review the multiresolution decomposition of a time varying operator and the in-scale maximum likelihood estimator is presented based on the conjugate gradient algorithm for in-scale recursions. We go on to introduce a suitable sparsity prior for the time varying channel and this leads to the posterior mean and variance along with estimators of the nuisance parameters associated with the sparsity model. Lastly, we present Doppler spread estimates and stopping rules based on the maximum a posterior (MAP) criteria and on Laplace&#39;s approximation as well as the simple and effective Akaike information criteria (AIC). Finally, the usefulness of these algorithms is demonstrated on a diverse set of simulated channels with various types and degrees of Doppler spread. Improved performance of the empirical Bayes posterior mean estimator over the maximum likelihood estimator is demonstrated. 
   The following definitions are provided for ease of reference, and will be used throughout. Without loss of generality all matrices and vectors are assumed real. Denote the M×N matrix A by A M×N  and the N×1 vector {right arrow over (a)} by {right arrow over (a)} N . Let {right arrow over (1)} N  denote a column vector of 1&#39;s, {right arrow over (0)} N  denote a column vector of 0&#39;s and I N  denote the identity matrix of size N by N. 
   Definitions: 
   
       
       
         
           1. For vectors {right arrow over (a)} N  define the diagonal matrix with diagonal elements {right arrow over (a)} N  by Diag({right arrow over (a)}). For square matrices A N,N , the vector of diagonal elements of the matrix A by diag(A). 
           2. Stack operator, [•] st  acts on a matrix A M×N =[{right arrow over (a)} 1 , {right arrow over (a)} 2 , . . . {right arrow over (a)} N ] yielding [A] st =[{right arrow over (a)} 1 ′, . . . {right arrow over (a)} N ′]′, a stacking of A&#39;s component columns into a supercolumn with inverse denoted by [•] ist  such that A=[A st ] ist . 
           3. Kronecker product A M×N {circle around (×)}B L×K  is the L M×N K matrix of weighted B blocks, having {m,n} th  block entry a m,n B. 
           4. Elementwise product, (A M×N {circle around (·)}B M×N ) m,n =a m,n b m,n . 
           5. Vector square root, 
         
       
     
  
   
     
       
         
           
             
               ( 
               
                 
                   a 
                   → 
                 
               
               ) 
             
             n 
           
           = 
           
             
               
                 a 
                 n 
               
             
             . 
           
         
       
     
       
       
         
           6. Matrix operator {hacek over (H)}, associated with the matrix array H M×N  with n th  column H n , has N+M−1 rows and N columns; the n th  column is {hacek over (H)} n =[{right arrow over (0)} n−1 , H n , {right arrow over (0)} N−n ].
 
Paragraphs 58 through 61 list a few useful properties and identities that follow directly from these definitions.
 
         
       
     
  
   The empirical Bayes (EB) approach is useful primarily for its computational efficiency. Consider computation of E[h|r], var[h|r] and E[φ|r] where r represents received data, h a parameter set of interest and φ a set of nuisance parameters. Expectations over the joint density p(h, φ|r) are often not easily computed. EB methods leverage a crude assumption p(φ|r)≈δ(E[φ|r]) to approximate expectations via iteration E[h|r] l ≈E[h|r, φ=E[φ|r] l−1 ] and are synonymous with approximations of the posterior marginal density via p[h|r]≈p[h|r, φ=E[φ|r]]. The method is versatile with other suitable estimators in place of {circumflex over (φ)}=E[φ|r] l−1  where their computation is faster or warranted by knowledge of the distribution. The penalty of this approximation is its underestimation of variance since var[h|r]=E φ [var[h|φ, r]]+var φ [E[h|φ, r]]. EB approaches therefore provide a lower bound, var[h|r]≧E φ [var[h|r, φ]] that is useful when var φ [E[h|φ, r]] is relatively small. For adaptive filtering problems where the computational demands of Bayesian analysis are presently out of reach the approach has merit. 
   Define two multiresolution (M R) spaces one over time and the other over delay. Let {Ω j   ξ } j∈Z  and {Ω k   ψ } k∈Z  represent these respectively; thus Ω J   ξ ={Σ j=0   J Σ l∈Z c j,l ξ j,l (t) ∥{right arrow over (c)}∥ 2 &lt;∞} where ξ 0,0 (t)=Σf l   ξ ξ 0,0 (2t−l) is the scaling function associated with this MR space, and the lowpass f ξ  and highpass g ξ  are a QMF pair. The scale index j increases with increasing bandwidth in agreement with the convention on the wavelets that ξ j,l (t)=√{square root over (2 j−1 )}ξ 1,0 (2 j−1 t−l) and ξ 1,0 (t)=Σg l   ξ ξ 0,0 (2t−l). 
   Following the filter model of Doroslovacki, H B (t, τ)∈Ω J   ξ ×Ω K   ψ  implies 
                       H   B     ⁡     (     t   ,   τ     )       =       ∑     j   ,   k       J   ,   K       ⁢       ∑     l   ,   m         L   j     ,     M   k         ⁢       w     j   ,   k   ,   l   ,   m       ⁢       ξ     j   ,   l       ⁡     (   t   )       ⁢       ψ     k   ,   m       ⁡     (   τ   )               ⁢     
     ⁢               L   j     =     BT   ⁢           ⁢     2     j   -   J   -   1                   M   k     =     F   ⁢           ⁢     τ   max               ⁢     2     k   -   K   -   1                 (   1   )               
has approximate Doppler spread B∝2 J /T and represents the time varying filter perfectly on the domain (t, τ)=(0, T)×(0, τ max ). The model (1) implicitly specifies the maximal Doppler spread and is therefore denoted with a subscript B.
 
   Frequency selectivity in the band [2 K−k−1 , 2 K−k ]F is modeled by the Fτ max 2 k−K  basis functions {ψ k,m } k,m  at delay shifts 2 K−k+1 m/F. In this way the frequency selectivity of the channel due to diverse scatterer locations and delay spreads is modeled as a superpostion of wavelet bases. 
   For moving scatterers Doppler spread is induced and is synonymous with an effective channel modulation over time corresponding to an imparted bandwidth to the impulse response process over time at any given delay. Modulations at time t≈2 −j−1 l T with durations ∝T/2 j  are modeled with the basis functions {ξ j,l (t)} j,l . These modulations correspond to each scatterers motion such that scatterers corresponding to greater accelerations will yield a channel operator with projections onto bases at fine scales (i.e. large J). 
   Let H B (n, m)=H(n/2 q B, m/F) q&gt;1 be a suitably sampled in time t and delay τ version of H(t, τ). Representing this 2-D array as a single stacked column {right arrow over (h)} B  of time invariant filters via paragraph 21, definition 2 each operating to produce F/2 q B samples of the output, express the wavelet coefficients of H(t, τ) (nonzero up to scale J) by
 
 {right arrow over (w)}   J   =U   J ( D   ξ   −q   {circle around (×)}I   2     k   ) {right arrow over (h)}   B               W   J =Ψ′ K   H   B   D   ξ   q Ξ J .  (2)
 
The matrix D ξ   q  is the 2 J+q  by 2 J  matrix decimation operator associated with the M R space of ξ. Its adjoint D ξ   q′ =D ξ   −q  is the associated interpolation operator. The columns of Ξ J  and Ψ K  are respectively the expansion coefficients of the wavelets ξ j,l (t) and ψ k,m (τ) in the bases of the associated scaling functions ξ 0,0 (t) and ψ 0,0 (τ) respectfully. The matrix Ξ J ∈R BT×BT  and Ψ K ∈R Fτ     max     ×Fτ     max    have maximal scales of J=└log 2  B T┘ and K =└ log 2  Fτ max ┘ respectively. The operator U J  is the Kronecker product of the two wavelet transforms U J =Ξ′ J {circle around (×)}Ψ′ K . The wavelet coefficients {right arrow over (w)} J ={w j,k,l,m , j=1 . . . J, k=1 . . . . K, l, m∈Z 2 } are computable via the fast wavelet transform. The operator H is synthesized from {right arrow over (w)} J  with the fast inverse wavelet transform to J+q scales (q interpolation/scaling operators added via zero appending W J )
 
{right arrow over (h)} B =( D   ξ   q   {circle around (×)}I   2     K   ) U′   J   {right arrow over (w)}   J             H   B =Ψ K   W   J Ξ′ J   D   ξ   −q   (3)
 
implying that H B  is synthesized as constant at time scales less than 1/2 q B. Clearly then q must satisfy q&lt;log 2  F/B since the channel is LTI at the sampling rate. The LHS of (3) follows from paragraph 58, property 4, with Ξ and Ψ unitary. Compute {right arrow over (h)} B  from {right arrow over (w)} J  via the RHS of (3) and denote the scaling function expansion of the channel by
 
{right arrow over (h)} B   0 =U′ J {right arrow over (w)} J .  (4)
 
Lastly define the time varying convolution matrix operator associated with the time varying array structure H B  or synonymously its stacked representation {right arrow over (h)} B  as {hacek over (H)} B .

   For channels with sparse arrivals that are dispersed and non-stationary the wavelet model offers a represention that is parsimonious. However for filters with peaky features in frequency, Fourier bases or autoregressive models will be more suitable. Fourier bases modeling Doppler spread will be optimal as well for the WSS channel as T→∞. Fourier bases will fail however to capture discontinuities or time localized phenomena. In addition where finite duration segments must be modeled boundaries pose severe penalties to Fourier methods, wavelets on the other hand parsimoniously capture these features. At lower frequencies associated with scatterers with smaller Doppler spread wavelets will offer near Fourier-like performance for the WSS condition. Nevertheless wavelets provide a means to model abrupt changes in channel conditions, localized phenomena in time and delay as well as minimizing boundry effects associated with finite duration signaling. Wavelet models of course are not universally appropriate for delay localized operators. For instance in the case of channels consisting of arrivals that have no dispersion in delay the standard Euclidean bases (in delay) are optimal (i.e. H(t, τ)=Σ n α n (t)δ τ−τ     n     (t) ). In this case the path delays and amplitudes (associated with the time changing geometry of the environment) are modeled. 
   The following paragraphs introduce an in-scale likelihood model from which a recursive in-scale MLE follows. A sparsity prior is then assumed and from this in-scale posterior expectation, variance and Doppler spread B estimates are derived. Empirical Bayes methods are employed throughout. 
   The channel response H(t, τ) of (1) with B Hz Doppler spread is observed via the source, s(t) and received signal r(t) through the linear model
 
 r ( t )=∫ H ( t , τ) s ( t −τ) dτ+n ( t )  (5)
 
where n(t) is a white Gaussian noise process of known power σ 2 F.
 
   Let N x (μ, Σ)=√{square root over ((2π) −N |Σ| −1 )}exp(−(x−μ)′Σ −1 (x−μ)/2). Assuming H is time-invariant over durations 1/2 q B and that s(t) and r(t) are suitably sampled at the Nyquist rate F let N=F(T+τ max ) the received signal dimension and M=Fτ max  the delay spread dimension. The discrete time model associated with (5) using the time varying operator of paragraph 21, definition 6, is
 
 p ( {right arrow over (r)}|{right arrow over (h)}   B   ,{right arrow over (s)} ,σ)= N   {right arrow over (r)} ( {hacek over (H)}   B   {right arrow over (s)},Iσ   2 )  (6)
 
where {hacek over (H)} B  represents the N by N−M+1 block convolution operator associated with {right arrow over (h)} B  in (3). Synonymously using (3) express (6) as a set of stacked wavelet coefficients
 
 p ( {right arrow over (r)}|{right arrow over (w)}   J   ,{right arrow over (s)} ,σ)= N   {right arrow over (r)} ( S   2     q     B ( D   ξ   q   {circle around (×)}I   2     K   ) U′   J   {right arrow over (w)}   J   , Iσ   2 )  (7)
 
where J is interchangeable with B as in (3) and for notational simplicity the conditioning on families ξ and ψ are assumed for time and delay respectively. With P=BT the Doppler dimension, S 2     q     B  is the N×2 q P M block convolution operator associated with the sampled source signal {right arrow over (s)}. The blocks of S 2     q     B  are proportional to the coherence time of 1/2 q B over which the channel is LTI. S 2     q     B  operates on a 2 q P long list of stacked M length LTI channel vectors to yield a channel output. To illustrate consider q=0 then S B  has the form
 
                   S   B     =     [           L   1         0       ⋮       ⋮       ⋮       0             X   1         0                                           ⋮             Y   1           L   2                                             ⋮           0         X   2                                             ⋮           0         Y   2         •                               ⋮           ⋮       0                                                           ⋮                               •                   0           ⋮                                                       0           ⋮                                           •         L   P                                                             0         X   P             0       ⋮       ⋮       ⋮       0         Y   P           ]             (   8   )               
The L&#39;s (Y&#39;s) are lower (upper) triangular, and are each of size M by M. The X&#39;s are full and of size F/B−M by M. Together each [L′, X′, Y′]′ block represents a time invariant filter operator over a duration of 1/2 q B seconds derived from the source signal of this block. The operator (D ξ   q {circle around (×)}I 2     K   ) maps the N by 2 K+J+q  block convolution operator S 2     q     B  of source samples, to the N by 2 K+J  block convolution operator of modulated blocks associated with coefficients in a scaling expansion of ξ. Some useful properties are exposed by considering the case of q large. Define the limiting matrix
   S   B   ξ =lim q→∞   S   2     q     B ( D   ξ   q   {circle around (×)}I   2     K   )  (9) 
which operates on the stacked coefficients of a scaling function expansion of the channel (i.e. {right arrow over (h)} B   0 =U′ J {right arrow over (w)} J ) to produce the expected value of the output of the time varying filter according to (7). The following scaling equation for such block convolution operators is informative S B/2   ξ =S B   ξ (D ξ   1 {circle around (×)}I 2     K   ). Further observe that for ξ Haar the scaling equation for the operator S B  takes the form
   S   B   =S   2B ( D   Haar   1   {circle around (×)}I   2     K   )  (10) 
and while it requires no multiplies it exhibits abrupt discontinuities between blocks. In practice (9) is not computable, the upper bound on q is determined by the bandwidth of the source q&lt;log 2 F/B such that modeling of the channel as anything but LTI over durations smaller than the sampling period is infeasible and unnecessary. For channels with reasonable regularity q&gt;3 is sufficient in practice. For instance with q=5, a B Hz Doppler spread channel is modeled as LTI over durations &lt;1/32 B sec.
 
   Returning back to (7) the implied likelihood function on 
               w   →     J     =       U   J     ⁢       h   →     B   0             
is
 
                     p   ⁡     (         r   →     |     w   →       ,     s   →     ,     J   =       log   2     ⁢   BT         )       ∝       ??       w   ~     J       ⁡     (           w   →     ~     J     ,     Δ   J       )         ⁢     
     ⁢                   w   →     ~     J     =       U   J     ⁢         h   →     ~     B   0         ,                   h   →     ~     B   0     =         (     C   B   ξ     )       -   1       ⁢     S   B     ξ   ′       ⁢     r   →         ,             Δ   J     =       σ   2     ⁢     U   J     ⁢     C   B   ξ     ⁢     U   J   ′                       (   11   )               
where the 2 K+J ×2 K+J  matrix C B   ξ =S B   ξ S B   ξ  is the source covariance matrix associated with the time-varying filter model. To illustrate consider again the simple case q=0 and ξ Haar. It follows that C B =S′ B S B  and takes the form
 
                   C   B     =     [           T   1           V   1         0                   …       0             V   1   ′           T   2                                 ⋰       ⋮           0                   •                               0           0                               •                   0           ⋮       ⋰                                 T     P   -   1             V     P   -   1               0       …                   0         V     P   -   1     ,           T   P           ]             (   12   )               
with T p  M×M Toeplitz, V p  lower triangular Hankel with zeros on the main diagonal having rank M−1 and P=BT. C B  is represented either by the 2P vectors of length M corresponding to the first columns of T p  and V p , or by the 2M+1 diagonals of C B  each being of length P M. Thus C B  is non-Toeplitz, banded Hermitian with 2M+1 nonzero diagonals. For ξ≠Haar C B   ξ  will differ from C B  in that it will be banded with width &gt;2M+1. The vector
 
               b   →     B     =       S   b     ξ   t       ⁢     r   →             
in (11) is a block overlapped cross-correlation of source and received data vectors associated with the windowing scaling function ξ 0,0 (t). The weakness of the in-scale likelihood approach is simply that the parameter set grows exponentially with scale, (i.e. PM∝2 J ) a severe penalty for a linear model. Nevertheless for the very limited set of spectrally flat signals the approach has some merit computationally. Signals for which ∥C B   ξ −C 0 I PM ∥ is small for some scalar C 0  (e.g. Gaussian white signals, maximal length sequences) ∥{right arrow over (b)} B −{right arrow over (h)} B ∥ will be small as well and close approximations to {right arrow over (h)} B  are attained with the fast block correlation.
 
   
     
       
         
           
             
               b 
               → 
             
             B 
           
           = 
           
             
               S 
               B 
               
                 ξ 
                 ′ 
               
             
             ⁢ 
             
               
                 r 
                 → 
               
               . 
             
           
         
       
     
   
   The in-scale computation is done using the conjugate gradient method. A direct computation of 
               h   →     ~     B   0         
in (11) requires O((PM) 3 ) a computationally impractical feat. A feasible approach assumes that
 
               h   →     ~       B   /   2     0         
is known to a tolerable precision. Compute
 
               h   →     ~     B   0         
starting from
 
             (       D   ξ   1     ⊗     I     2   K         )     ⁢         h   →     ~       B   /   2     0           
via the conjugate gradient (CG) method. Initialization requires only O(M 2 ) since for J=0 the system is Toeplitz. The CG algorithm has found wide application in solving banded and sparse systems. Indeed the CG algorithm has been shown to be synonymous with the multistage Wiener filter. Table 1 lists the operations necessary along with the associated computational demands. Note that only a single matrix multiplication by the sparse banded matrix C B   ξ  is required per iteration. A multiply is accomplished efficiently from the 2M+1 stored diagonals of length P M. In this framework the inversion of large matrices is avoided by solving much lower (1/2) dimension problems to a given precision then using these as starting points for the next higher hypothesized dimension. The resulting maximum likelihood estimate when the CG algorithm is iterated to machine precision is termed
 
             h   →     B   mle         
and is computed via interpolation from
 
               h   →     ~     B   0         
as
 
   
     
       
         
           
             
               
                 h 
                 → 
               
               ~ 
             
             B 
             mle 
           
           = 
           
             
               
                 ( 
                 
                   
                     D 
                     ξ 
                     q 
                   
                   ⊗ 
                   
                     I 
                     
                       2 
                       K 
                     
                   
                 
                 ) 
               
               ⁢ 
               
                 
                   
                     
                       h 
                       → 
                     
                     ~ 
                   
                   B 
                   0 
                 
                 ⟺ 
                 
                   H 
                   B 
                   mle 
                 
               
             
             = 
             
               
                 
                   [ 
                   
                     
                       h 
                       → 
                     
                     B 
                     0 
                   
                   ] 
                 
                 ist 
               
               ⁢ 
               
                 
                   D 
                   ξ 
                   
                     - 
                     q 
                   
                 
                 . 
               
             
           
         
       
     
   
   The CG algorithm also allows for flexibilty in the number of iterations. Computation of {right arrow over (h)} B   mle  to machine precision may not be necessary and a solution slightly more coherent (closer to 
               h   →       B   /   2     mle     )         
may be tolerable. For less computation iterate to within a ball of radius proportional to the covariance of the MLE and call such a solution {right arrow over (h)} B   a , an approximate MLE. To determine this iteration number an approximate ad-hoc approach is taken. Define κ the ratio of minimum to maximum eigenvalues of C B   ξ  then at iteration i a , where
 
               h   →     B     {     i   α     }       =       h   →     B   α           
the CG solution is within a ball of radius
 
                      h   →     B   α     -       h   →     B   mle            =                ∑   PM         i   α     +   1       ⁢       ρ   i     ⁢       d   →     i              ≤     2   ⁢            h   →     B   mle          ⁢     κ     ⁢     ω     i   α                         where   ⁢           ⁢   ω     =       (       κ     -   1     )     /       (       κ     +   1     )     .             
The computation of κ is not trivial and must be approximated for the class of source signals that are considered. Further in the case of in-scale computation proposed here the extreme eigenvalues of C B   ξ  at each B would have to be approximated, at present this is impractical and an approximation is used for ω. The implementation of the CG algorithm presented makes an additional ad-hoc assumption to eliminate
 
                     ⁢       h   →     B   mle     ⁢                   
from the stopping criteria based on the fact that the CG algorithm gives fairly smooth exponential convergence due to dominant eigenvectors not dominating the convergence. For this reason in practice the tightness of the
 
           ≤     2   ⁢            h   →     B   mle          ⁢     κ     ⁢     ω     i   α               
bound does not vary greatly from iteration to iteration. Express the radius as
 
                      ∑   PM         i   α     +   1       ⁢       ρ   i     ⁢       d   →     i              =       G   ⁢              ρ   ⁢                 i   α     +   1       ⁢       d   →         i   α     +   1                +              ∑   PM         i   α     +   2       ⁢       ρ   i     ⁢       d   →     i                      
for some constant 0&lt;G&lt;1 (since the innovation components are not orthogonal). Assuming the convergence is fairly well behaved approximate
 
                         ⁢     ρ     i   a         +     1   ⁢       d   →       i   a         +     1   ⁢                 ≈                 ⁢       h   →     B   a       -         h   →     B   mle     ⁢             ⁢       (     1   -   ω     )     .               
A common scenario is that of a white non-fading source signal with variance ζ 2 . The variance of each coefficient w j,k,l,m  is simply δ J   2 =Δ j,k =σ 2 /ζ 2 ×B/F. With PM=BTFτ max  coefficients this implies {right arrow over (w)} J  has variance P×Bτ max /SNR. In this scenario, with G=10 iterate while ∥ρ i+1 {right arrow over (d)} i+1 ∥&gt;10(1−ω)PBτ max /SNR. A conservative approximation of κ=1000 corresponding to ω=0.93 can be made.
 
   
     
       
             
           
             
             
             
           
             
           
             
             
             
             
           
             
           
             
             
             
           
         
             
               TABLE 1 
             
             
                 
             
             
               
                 
                   
                     
                       Computation 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       of 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           
                             w 
                             ⇀ 
                           
                           ~ 
                         
                         J 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       via 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       conjugate 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         gradient 
                         . 
                       
                     
                   
                 
               
             
             
                 
             
           
           
             
               P = BT, M = Fτ max , N c  = F/B 
             
             
               m 0  ∝ length(g ξ ) 
             
             
               Initialize 
             
           
        
         
             
                 
               
                 
                   
                     
                       
                         
                           x 
                           ⇀ 
                         
                         0 
                       
                       = 
                       
                         
                           ( 
                           
                             
                               D 
                               ξ 
                               1 
                             
                             ⊗ 
                             
                               I 
                               K 
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           
                             
                               h 
                               ⇀ 
                             
                             ~ 
                           
                           
                             B 
                             / 
                             2 
                           
                           0 
                         
                       
                     
                   
                 
               
               O(PM × m o   2 ) 
             
             
                 
                 
             
             
                 
               
                 
                   
                     
                       
                         C 
                         B 
                       
                       , 
                       
                         
                           
                             b 
                             ⇀ 
                           
                           B 
                         
                         = 
                         
                           
                             S 
                             B 
                             ′ 
                           
                           ⁢ 
                           
                             r 
                             ⇀ 
                           
                         
                       
                     
                   
                 
               
               O(PN c logN c ) 
             
             
                 
                 
             
             
                 
               
                 
                   
                     
                       
                         
                           
                             d 
                             ⇀ 
                           
                           0 
                         
                         = 
                         
                           
                             e 
                             ⇀ 
                           
                           0 
                         
                       
                       , 
                       
                         
                           
                             e 
                             ⇀ 
                           
                           0 
                         
                         = 
                         
                           
                             b 
                             ⇀ 
                           
                           - 
                           
                             
                               C 
                               B 
                             
                             ⁢ 
                             
                               
                                 x 
                                 ⇀ 
                               
                               0 
                             
                           
                         
                       
                     
                   
                 
               
             
             
                 
                 
             
             
                 
               
                 
                   
                     
                       
                         
                           f 
                           ⇀ 
                         
                         0 
                       
                       = 
                       
                         
                           C 
                           B 
                         
                         ⁢ 
                         
                           
                             d 
                             ⇀ 
                           
                           0 
                         
                       
                     
                   
                 
               
               O(PM 2 ) 
             
           
        
         
             
               
                 
                   
                     
                       
                         while 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           
                              
                             
                               
                                 ρ 
                                 i 
                               
                               ⁢ 
                               
                                 
                                   d 
                                   ⇀ 
                                 
                                 i 
                               
                             
                              
                           
                           2 
                         
                       
                       &gt; 
                       
                         
                           ( 
                           
                             1 
                             - 
                             ω 
                           
                           ) 
                         
                         ⁢ 
                         PB 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           τ 
                           max 
                         
                         ⁢ 
                         
                           δ 
                           J 
                           2 
                         
                       
                     
                   
                 
               
             
           
        
         
             
                 
               
                 
                   
                     
                       
                         η 
                         = 
                         
                           
                             
                               e 
                               ⇀ 
                             
                             i 
                             ′ 
                           
                           ⁢ 
                           
                             
                               
                                 f 
                                 ⇀ 
                               
                               i 
                             
                             / 
                             
                               
                                 d 
                                 ⇀ 
                               
                               i 
                               ′ 
                             
                           
                           ⁢ 
                           
                             
                               f 
                               ⇀ 
                             
                             i 
                           
                         
                       
                       , 
                     
                   
                 
               
               
                 
                   
                     
                       
                         
                           d 
                           ⇀ 
                         
                         
                           i 
                           + 
                           1 
                         
                       
                       = 
                       
                         
                           
                             e 
                             ⇀ 
                           
                           i 
                         
                         - 
                         
                           η 
                           ⁢ 
                           
                             
                               d 
                               ⇀ 
                             
                             i 
                           
                         
                       
                     
                   
                 
               
               O(PM) 
             
             
                 
                 
             
             
                 
               
                 
                   
                     
                       
                         
                           f 
                           ⇀ 
                         
                         
                           i 
                           + 
                           1 
                         
                       
                       = 
                       
                         
                           C 
                           B 
                         
                         ⁢ 
                         
                           
                             d 
                             ⇀ 
                           
                           
                             i 
                             + 
                             1 
                           
                         
                       
                     
                   
                 
               
                 
               O(PM 2 ) 
             
             
                 
                 
             
             
                 
               
                 
                   
                     
                       
                         
                           ρ 
                           
                             i 
                             + 
                             1 
                           
                         
                         = 
                         
                           
                             
                               e 
                               ⇀ 
                             
                             i 
                             ′ 
                           
                           ⁢ 
                           
                             
                               
                                 d 
                                 ⇀ 
                               
                               
                                 i 
                                 + 
                                 1 
                               
                             
                             / 
                             
                               
                                 d 
                                 ⇀ 
                               
                               
                                 i 
                                 + 
                                 1 
                               
                               ′ 
                             
                           
                           ⁢ 
                           
                             
                               f 
                               ⇀ 
                             
                             
                               i 
                               + 
                               1 
                             
                           
                         
                       
                       , 
                     
                   
                 
               
               
                 
                   
                     
                       
                         
                           v 
                           ⇀ 
                         
                         
                           i 
                           + 
                           1 
                         
                       
                       = 
                       
                         
                           
                             v 
                             ⇀ 
                           
                           i 
                         
                         + 
                         
                           
                             ρ 
                             
                               i 
                               + 
                               1 
                             
                           
                           ⁢ 
                           
                             
                               f 
                               ⇀ 
                             
                             
                               i 
                               + 
                               1 
                             
                           
                         
                       
                     
                   
                 
               
               O(PM) 
             
             
                 
                 
             
             
                 
               
                 
                   
                     
                       
                         
                           
                             x 
                             ⇀ 
                           
                           
                             i 
                             + 
                             1 
                           
                         
                         = 
                         
                           
                             
                               x 
                               ⇀ 
                             
                             i 
                           
                           + 
                           
                             
                               ρ 
                               
                                 i 
                                 + 
                                 1 
                               
                             
                             ⁢ 
                             
                               
                                 d 
                                 ⇀ 
                               
                               
                                 i 
                                 + 
                                 1 
                               
                             
                           
                         
                       
                       , 
                     
                   
                 
               
               
                 
                   
                     
                       
                         
                           e 
                           ⇀ 
                         
                         
                           i 
                           + 
                           1 
                         
                       
                       = 
                       
                         
                           b 
                           ⇀ 
                         
                         - 
                         
                           
                             v 
                             ⇀ 
                           
                           
                             i 
                             + 
                             1 
                           
                         
                       
                     
                   
                 
               
               O(PM) 
             
           
        
         
             
               i = i + 1 
             
             
               end 
             
           
        
         
             
                 
               
                 
                   
                     
                       
                         
                           
                             h 
                             ⇀ 
                           
                           ~ 
                         
                         B 
                       
                       = 
                       
                         
                           x 
                           ⇀ 
                         
                         
                           i 
                           + 
                           1 
                         
                       
                     
                   
                 
               
             
             
                 
                 
             
             
                 
               
                 
                   
                     
                       
                         
                           
                             w 
                             ⇀ 
                           
                           ~ 
                         
                         J 
                       
                       = 
                       
                         
                           U 
                           J 
                         
                         ⁢ 
                         
                           
                             x 
                             ⇀ 
                           
                           
                             i 
                             + 
                             1 
                           
                         
                       
                     
                   
                 
               
               O(PM × m o   2 ) 
             
             
                 
                 
             
           
        
       
     
   
   The Gaussian (normal) mixture can be chosen as a sparsity prior for the doppler spread channel. The previous maximum likelihood in-scale estimate of the channel assumes that all delay lags of the channel at each time contribute to the output a priori equally. Synonymously there is in the likelihood model an implicit prior of large and equal variance for each and every time-delay component of the channel operator. This assumption is simplistic for sparse channels and a modification to the likelihood approach is warranted. The prior chosen should give the model flexibility to choose components, associated with frequency selectivity and Doppler spread (time modulation) that are significant and reject those that are not. This problem is synonymous with variable selection and mixture prior assignments have been useful for such problems in both the Bayesian and empircal Bayes methodologies. For this reason the binary discrete mixture normal model
 
 p ( w |π,λ,␣)= M (π,λ,␣)             w|π,λ,␣˜πN   w (0,λ 2 )+(1−π) N   w (0,∈ 2 )  (13)
 
is chosen as a prior. Specify the prior conditional density of {right arrow over (w)} J  as

                       w   →     J     ⁢     ❘     ⁢     π     j   ,   k         ,     λ     j   ,   k       ,       ε     j   ,   k       ~       ∏     j   ,   k       J   ,   K       ⁢       ∏     l   ,   m         L   j     ,     M   k         ⁢         M     w     j   ,   k   ,   l   ,   m         ⁡     (       π     j   ,   k       ,     λ     j   ,   k       ,     ε     j   ,   k         )       .                   (   14   )               
The mixture normal model can be useful because it captures an attribute of sparsity. For each of the significant arrival paths associated with an operator we have a normal model. For those delays for which no arrival energy is present we model these with small (∈) variance. The hyperparameters necessary to specify the model are summarized in a single parameter {right arrow over (φ)} J ={π j,k , λ j,k , ∈ j,k } j,k≦J,K . With the prior and likelihood specified the resulting posterior distribution to be maximized jointly over {right arrow over (w)} J  and J (synonymous with {right arrow over (h)} B  and B) and hyperparameters {right arrow over (φ)} J  is
 p({right arrow over (w)}, {right arrow over (φ)},J|{right arrow over (r)},{right arrow over (s)}, σ 2 )∝N {right arrow over (w)}     J   ({right arrow over ({tilde over (w)} J ({right arrow over (r)},{right arrow over (s)}), Δ J ({right arrow over (s)},σ 2 ))×p({right arrow over (w)}|{right arrow over (φ)},σ 2 ,J)×p({right arrow over (φ)}|J)×p(J)  (15) 
The first term, the likelihood assumes independence of wavelet coefficients via the near diagonalizing property of U J  on the source covariance C B   ξ  and the second term the prior, is predicated on this assumption.
 
   For ease of notation let *=l, m represent an arbitrary wavelet coefficient at time-delay location associated with l, m according to the model of Equation (1) so that w j,k,* =w j,k,l,m . In the spirit of empirical Bayes consider the posterior distribution of wavelet coefficients given all other parameters. This conditional density is shown in Equation 33 to be 
                     w     j   ,   k   ,   *       |       w   ~       j   ,   k   ,   *         ,     σ   2     ,             ϕ   →     J     ~     Π   ⁡     (       w   ~       j   ,   k   ,   *       )         ⁢       ??     w     j   ,   k   ,   *         ⁡     (         γ   ⁡     (     λ     j   ,   k       )       ⁢       w   ~       j   ,   k   ,   *         ,       γ   ⁡     (     λ     j   ,   k       )       ⁢     δ     j   ,   k     2         )         +       (     1   -     Π   ⁡     (       w   ~       j   ,   k   ,   *       )         )     ⁢       ??     w     j   ,   k   ,   *         ⁡     (         γ   ⁡     (     ε     j   ,   k       )       ⁢       w   ~       j   ,   k   ,   *         ,       γ   ⁡     (     ε     j   ,   k       )       ⁢     δ     j   ,   k     2         )                   (   16   )               
where the posterior mixture coefficient Π({right arrow over (w)} j,k,* ) and the gains γ(·) are defined in the derivation of Equation 33. This posterior (16) has the same mixture normal form as the prior (13). There are however two important distinctions, first the posterior mixture coefficents Π(·) are functions of each individual empirical coefficient {right arrow over (w)} j,k,l,m . Secondly the means of the mixture densities are not equal. For this reason the posterior density is not symetric and therefore its first moment ŵ=E[w|{tilde over (w)},{right arrow over (φ)} J ] is not synonymous with the argument maximizing (16).
 
   Since the posterior is a mixture the mean is simple to assess as 
                             w   ^       j   ,   k   ,   *       =     E   ⁡     [         w     j   ,   k   ,   *       |       w   ~       j   ,   k   ,   *         ,     σ   2     ,       ϕ   →     J       ]                   =       [       Π   ⁢           ⁢     (       w   ~       j   ,   k   ,   *       )     ⁢     γ   ⁡     (     λ     j   ,   k       )         +       (     1   -     Π   ⁡     (       w   ~       j   ,   k   ,   *       )         )     ⁢     γ   ⁡     (     ε     j   ,   k       )           ]     ⁢       w   ~       j   ,   k   ,   *                 .           (   17   )               
The resulting operator
 
               w   →     J     =     E   ⁡     [         w   →     ❘       w     →   ~       J       ,     σ   2     ,       ϕ   →     J       ]             
attenuates smaller coefficients and leaves larger ones relatively unchanged. It is similar in form to a Wiener filter with the noticeable difference of a mixture of Wiener gains, each gain being proportional to the implicit empirical posterior probability of the associated model in the mixture.
 
   The posterior mean of the channel {right arrow over (h)} B  sampled at the rate 2 q B Hz follows directly from (17) and the linearity of the wavelet transform as
 
 {right arrow over (ĥ)}   B   =E[{right arrow over (h)}|{right arrow over (r)}, {right arrow over (s)}, σ   2 ,{right arrow over (φ)} J ]=( D   q   {circle around (×)}I   2     K   ) U   J {right arrow over (ŵ)} J               Ĥ   B =Ψ K   [{right arrow over (ŵ)}   J ] ist Ξ′ J   D   ξ   −q .  (18)
 
The posterior mean is useful since for the ensemble of channels associated with the mixture prior on wavelet coefficients no other estimator has smaller mean square error. It is however a biased estimator. The MAP estimate for the posterior is computationally intensive requiring iterations over every wavelet coefficient and for this reason an approximate MAP is worth considering.

                   w     j   ,   k   ,   *     amap     =     {             γ   ⁡     (     λ     j   ,   k       )       ⁢       w   ~       j   ,   k   ,   *                 if   ⁢           ⁢     Π   ⁡     (       w   ~       j   ,   k   ,   l   ,   m       )         &gt;     1   /   2                   γ   ⁡     (     ε     j   ,   k       )       ⁢       w   ~       j   ,   k   ,   *                 if   ⁢           ⁢     Π   ⁡     (       w   ~       j   ,   k   ,   l   ,   m       )         &lt;     1   /   2                       (   19   )               
is a simple approximation and has been tested against the actual MAP, computed via Nelder-Mead iteration, and gives comparable performance in terms of MSE on simulated channels.
 
   The conditional variance of the channel {right arrow over (h)} B  is computable directly from the posterior density (16) and gives measure of the uncertainty associated with the channel after observation of the data. The derivation of Equation 34 presents proof of the following 
                         v     j   ,   k   ,   *     2     =     var   ⁡     [         w     j   ,   k   ,   *       |       w   ~       j   ,   k   ,   *         ,       ϕ   →     J       ]                   =         (         Π   ⁡     (       w   ~       j   ,   k   ,   *       )       ⁢     γ   ⁡     (     λ     j   ,   k       )         +       (     1   -     Π   ⁡     (       w   ~       j   ,   k   ,   *       )         )     ⁢     γ   ⁡     (     ε     j   ,   k       )           )     ⁢     δ     j   ,   k     2       +                 =       Π   ⁡     (       w   ~       j   ,   k   ,   *       )       ⁢     (     1   -     Π   ⁡     (       w   ~       j   ,   k   ,   *       )         )     ⁢       (       γ   ⁡     (     λ     j   ,   k       )       -     γ   ⁡     (     ε     j   ,   k       )         )     2     ⁢         w   ~       j   ,   k   ,   *     2     .                     (   20   )               
The above shows that the conditional variance for each coefficient is the sum of two components. The first is a weighted average of the variances associated with the two models of the mixture. The last component represents the uncertainty of the coefficient being associated definitively with either of the two models in the mixture and is maximum for Π(w)=½. For this reason large variance is associated with coefficients with magnitude approximately
 
               w   ~       j   ,   k   ,   *       =     ±               (       λ     j   ,   k     2     +     δ     j   ,   k     2       )     ⁢     (       ε     j   ,   k     2     +     δ     j   ,   k     2       )           λ   2     -     ε   2         ⁢     log   ⁡     (         (       λ     j   ,   k     2     +     δ     j   ,   k     2       )       (       ε     j   ,   k     2     +     δ     j   ,   k     2       )       ⁢         (     1   -     π     j   ,   k         )     2       π     j   ,   k     2         )           .             
Define {right arrow over (v)} J  as the vector of posterior variances (20) associated with {right arrow over (w)} J ; then under the assumption that the wavelet coefficients are uncorrelated, cov[{right arrow over (w)} J ]=Diag[{right arrow over (v)} J ]. The posterior covariance of {right arrow over (h)} B  is then, using (3) and paragraph 58, property 2,
 
                         Γ   B     =     cov   ⁡     [       h   →     B     ]                   =       (       D   ξ   q     ⊗     I     2   κ         )     ⁢     U   J   ′     ⁢     Diag   ⁡     (       v   →     J     )       ⁢       U   J     ⁡     (       D   ξ     -   q       ⊗     I     2   κ         )                     =       (       D   ξ   q     ⁢       Ξ   J     ⊗     Ψ   K         )     ⁢     Diag   ⁡     (       v   →     J     )       ⁢       (       Ξ   J   ′     ⁢       D   ξ     -   q       ⊗     Ψ   K   ′         )     .                     (   21   )               
Determine the variance of H B  at any give t, τ by computing V B =[diag(Γ B )] ist  using paragraph 58 to express diag(Γ B )=[(D ξ   q Ξ J {circle around (·)}D ξ   q Ξ J ){circle around (×)}(Ψ K {circle around (·)}Ψ K )]{right arrow over (v)} J . This implies that the posterior variance field V B  defined over t, τ is
   V   B =(Ψ K {circle around (·)}Ψ K )[{right arrow over (v)} J ] ist ( D   ξ   q Ξ J   {circle around (·)}D   ξ   q Ξ J ).  (22) 
Since the operators (Ψ K {circle around (·)}Ψ K ) and (D ξ   q Ξ J {circle around (·)}D ξ   q Ξ J ) do not obey a scaling recursion they must be computed directly an impractical feat for channels of high dimension BT×Fτ max . For this reason a fast approximate (1−α)% confidence interval for H B  is worth pursuing. Proceed as follows; define
 
               w   →     J   α     =         w     →   ^       J     +       z   α     ⁢       v   J       ⁢           ⁢   and   ⁢           ⁢       w   →     J     1   -   α       ⁢       w     →   ^       J       -       z   α     ⁢         v   →     J                 
and transform via (3). Here of course P[z&gt;z α ]=α/2, with z˜N(0, 1). These wavelet domain bounds produce an approximate EB high probability region for the channel as
 
   
     
       
         
           
             
               
                 
                   
                     P 
                     ⁡ 
                     
                       [ 
                       
                         
                           
                             
                               
                                 
                                   min 
                                   ⁢ 
                                   
                                     ( 
                                     
                                       
                                         
                                           
                                             H 
                                             ^ 
                                           
                                           B 
                                           α 
                                         
                                         ⁡ 
                                         
                                           ( 
                                           
                                             t 
                                             , 
                                             τ 
                                           
                                           ) 
                                         
                                       
                                       , 
                                       
                                         
                                           
                                             H 
                                             ^ 
                                           
                                           B 
                                           
                                             1 
                                             - 
                                             α 
                                           
                                         
                                         ⁡ 
                                         
                                           ( 
                                           
                                             t 
                                             , 
                                             τ 
                                           
                                           ) 
                                         
                                       
                                     
                                     ) 
                                   
                                 
                                 &lt; 
                               
                               ⁢ 
                               
                                   
                               
                             
                           
                         
                         
                           
                             
                               
                                 H 
                                 ⁢ 
                                 
                                   ( 
                                   
                                     t 
                                     , 
                                     τ 
                                   
                                   ) 
                                 
                               
                               &lt; 
                               
                                 max 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       
                                         
                                           H 
                                           ^ 
                                         
                                         B 
                                         α 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           t 
                                           , 
                                           τ 
                                         
                                         ) 
                                       
                                     
                                     , 
                                     
                                       
                                         
                                           H 
                                           ^ 
                                         
                                         B 
                                         
                                           1 
                                           - 
                                           α 
                                         
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           t 
                                           , 
                                           τ 
                                         
                                         ) 
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                       ] 
                     
                   
                   ≈ 
                   
                     1 
                     - 
                     α 
                   
                 
                 ⁢ 
                 
                   
 
                 
                 ⁢ 
                 
                   
                     
                       
                         
                           
                             
                               
                                 
                                   H 
                                   ^ 
                                 
                                 B 
                                 α 
                               
                               = 
                               
                                 
                                   
                                     
                                       Ψ 
                                       K 
                                     
                                     ⁡ 
                                     
                                       [ 
                                       
                                         
                                           w 
                                           → 
                                         
                                         J 
                                         α 
                                       
                                       ] 
                                     
                                   
                                   ist 
                                 
                                 ⁢ 
                                 
                                   Ξ 
                                   J 
                                   ′ 
                                 
                                 ⁢ 
                                 
                                   D 
                                   ξ 
                                   
                                     - 
                                     q 
                                   
                                 
                               
                             
                           
                           
                             
                               
                                 
                                   
                                     H 
                                     ^ 
                                   
                                   J 
                                   
                                     1 
                                     - 
                                     α 
                                   
                                 
                                 = 
                                 
                                   Ψ 
                                   K 
                                 
                               
                               ⁢ 
                               
                                   
                               
                             
                           
                         
                       
                       ⁡ 
                       
                         [ 
                         
                           
                             w 
                             → 
                           
                           J 
                           
                             1 
                             - 
                             α 
                           
                         
                         ] 
                       
                     
                     ist 
                   
                   ⁢ 
                   
                     Ξ 
                     J 
                     ′ 
                   
                   ⁢ 
                   
                     
                       D 
                       ξ 
                       
                         - 
                         q 
                       
                     
                     . 
                   
                 
               
             
             
               
                 ( 
                 23 
                 ) 
               
             
           
         
       
     
   
   In keeping with the empirical Bayes methodology, estimate the hyperparameters {right arrow over (φ)} J={π   j,k ,λ j,k ,∈ j,k } as follows: first for π j,k  use Donoho&#39;s threshold argument, namely that z i ˜N(0, δ 2 ) implies that
 
lim n→∞   P [max i≦n   |z   i |&gt;δ√{square root over (log(2 πn ))}]=0. Let {right arrow over (1)} S ( x )={1 for  x∈S,  0 otherwise}.
 
With the total number of coefficients at the j,k th  scale denoted N j,k =2 K−k+J−j BTFτ max  then an asymptotically unbiased estimate of π j,k  when λ j,k &gt;&gt;δ j,k  for the model of (13) is
 
               π   ^       j   ,   k       =       ∑       1   →       S     j   ,   k             N     j   ,   k                   where               S     j   ,   k       =     {       w     j   ,   k   ,   l   ,   m       :            w     j   ,   k   ,   l   ,   m            &gt;       δ     j   ,   k       ⁢       log   ⁡     (     2   ⁢   π   ⁢           ⁢     N     j   ,   k         )               }           
For λ j,k  use {circumflex over (λ)} j,k =max * |{right arrow over (w)} j,k,* |/3. Lastly for the ∈&#39;s use the assumptions that the channel w j,k,*  and noise processes are independent with π&lt;&lt;1/2 the median absolute deviation (MAD) estimator of zero-mean iid (independent identically distributed) variates suggests E 2 [mad * {right arrow over (w)} j,k,* ]≈δ j,k   2 +∈ j,k   2  leading to the estimator {circumflex over (∈)} j,k   2 =max[δ j,k   2 /20, [mad * {right arrow over (w)} j,k,* ] 2 −δ j,k   2 ].
 
   In this and the following two paragraphs, inference regarding Doppler spread B or maximum scale J=log 2 BT is shown to be synonymous with model order selection. An EB-MAP estimate as well as rules based on Laplace&#39;s approximation and Akaike (AIC) are provided. It is convenient in keeping with the in-scale paradigm to place a prior on B synonymous with an uninformative prior on J, the maximum scale parameter. With this assumption the MAP estimate of the Doppler spread is determined by maximizing the likelihood p({right arrow over (r)}|J); the objective then is 
               J   ^     =     arg   ⁢           ⁢       min   J     ⁢     {       -   log     ⁢           ⁢     p   ⁡     (       r   →     ❘   J     )         }           ,     
     ⁢       p   ⁡     (         r   →     ❘   B     =       2   J     ⁢     /     ⁢   T       )       =     ∫       p   ⁡     (       r   →     ❘       w   →     J       )       ⁢     p   ⁡     (         w   →     J     ❘       ϕ   →     J       )       ⁢     p   ⁡     (         ϕ   →     J     ❘   J     )       ⁢     ⅆ     w   →       ⁢       ⅆ     ϕ   →       .                 
The first term in the integrand is decomposed as
 
             p   ⁡     (       r   →     ❘       w   →     J       )       =         (       2   ⁢   π       )         2   J     ⁢     Fτ   mzx         ⁢            Δ   J            1   2       ×       N     r   →       ⁡     (           H   ⋓     B   mle     ⁢     s   →       ,     I   ⁢           ⁢     σ   2         )       ×       N       w   →     J       ⁡     (         w     →   ~       J     ,     Δ   j       )               
leading to
 
             p   ⁡     (         r   →     |   B     =       2   J     /   T       )       =         ??     r   →       ⁡     (           H   ⋓     B   mle     ⁢     s   →       ,     I   ⁢           ⁢     σ   2         )       ×       (     2   ⁢   π     )       BTF   ⁢           ⁢       τ   max     /   2         ⁢            Δ   J            1   /   2       ×     ∫       ∏     j   ,   k   ,   l   ,   m       ⁢           ⁢       (         ??     w     j   ,   k   ,   l   ,   m         ⁡     (         w   ~       j   ,   k   ,   l   ,   m       ,     δ     j   ,   k     2       )       ×     [         π     j   ,   k       ⁢       ??     w     j   ,   k   ,   l   ,   m         ⁡     (     0   ,     λ     j   ,   k     2       )         +       (     1   +     π     j   ,   k         )     ⁢       ??     w     j   ,   k   ,   l   ,   m         ⁡     (     0   ,     ε     j   ,   k     2       )           ]     ⁢     p   ⁡     (       ϕ   →     |   J     )         )     ⁢     ⅆ     w   →       ⁢     ⅆ     ϕ   →                     
Integrate over {right arrow over (w)} analytically since N w ({right arrow over (w)}, δ 2 )N w (0, λ 2 )=N w (γ(λ){right arrow over (w)}, γ(λ)δ 2 )N {right arrow over (w)} (0, δ 2 +λ 2 ) and over the 3J K parameters of {right arrow over (φ)} J  using the asymptotic result of Schwarz (i.e. ∝ (FT) −3JK/2 ). Under the near diagonalizing property of U J  approximate |Δ J |≈Π j,k  δ j,k2 , and with the EB estimate of {right arrow over (φ)} J  and {right arrow over ({tilde over (e)} r,B =r−{hacek over (H)} B   mle {right arrow over (s)}, the MAP criterion is
   {circumflex over (B)}   MAP   =arg min J   {Q   MAP ( J )}  (25) 
where
 
                     Q   MAP     ⁡     (   J   )       =               ⅇ   →     ~       r   ,   B     ′     ⁢         ⅇ   →     ~       r   ,   B           2   ⁢     σ   2         -         (       2   J     ⁢   F   ⁢           ⁢     τ   max       )     2     ⁢     log   ⁡     (     2   ⁢     πδ     j   ,   k     2       )         +         3   ⁢   JK     2     ⁢   log   ⁢           ⁢   FT     -       ∑     j   ,   k   ,   l   ,   m       J   ,   K       ⁢       log   ⁡     (       ℳ       w   ~       j   ,   k   ,   l   ,   m         ⁡     (         π   ^       j   ,   k       ,       λ   ^       j   ,   k     2     ,     +           ⁢     δ     j   ,   k     2       ,         ε   ^       j   ,   k       +     δ     j   ,   k     2         )       )       .                 (   26   )               
It is somewhat striking and counterintuitive that this empirical Bayes MAP estimate of Doppler spread (25) does not require the computation of the MAP or posterior mean (PM), the residual errors being associated with that of the MLE.
 
   To consider rules based on the residual errors associated with the posterior mean or approximate MAP channel estimate proceed as follows; let lp({right arrow over (w)})=log p({right arrow over (r)}|{right arrow over (w)})p({right arrow over (w)}|{right arrow over (φ)}, J) (noting that this is proportional to the log posterior of {right arrow over (w)} expand about {right arrow over (w)} map ; lp({right arrow over (w)})≈lp({right arrow over (w)} map )+({right arrow over (w)}−{right arrow over (w)} map )′lp″({right arrow over (w)} map )({right arrow over (w)}−{right arrow over (w)} map )/2. Approximate lp″({right arrow over (w)} map )≈−V −1  and substitute the posterior mean {right arrow over (ŵ)} J  and covariance V B  yielding 
             p   ⁡     (       r   →     ❘   J     )       ≈       p   ⁡     (       r   →     ❘       h     →   ^       B       )       ⁢            2   ⁢   π              BTFτ   max     /   2       ⁢            V   B            1   2       ⁢     ∫       p   ⁡     (           w     →   ^       J     ❘       ϕ   →     J       ,   J     )       ⁢     p   ⁡     (       ϕ   →     ❘   J     )       ⁢       ⅆ     ϕ   →       .                 
With |V B | 1/2 =Π{right arrow over (v)} J   1/2  and the asymptotic result of Scharwz applied to {right arrow over (φ)} J  the resulting
 
Laplace approximation rule is
   {circumflex over (B)}   LAP   =arg min J   {Q   LAP ( J )}  (27) 
where
 
                     Q   LAP     ⁡     (   J   )       =         -   log     ⁢           ⁢     p   ⁡     (       r   →     |         w   →     ^     J       )         -     log   ⁢           ⁢     p   ⁡     (         w   →     ^     |   J     )         +         3   ⁢   JK     2     ⁢   log   ⁢           ⁢   FT     -       1   2     ⁢       ∑     j   ,   k   ,   l   ,   m       ⁢     log   ⁢           ⁢   2   ⁢   π   ⁢           ⁢       v     j   ,   k   ,   l   ,   m     2     .                     (   28   )               
The first term
 
             -   log     ⁢           ⁢     p   ⁡     (       r   →     ❘       w     →   ^       J       )             
is a measure of the quality of the posterior mean channel estimate to predict the data {right arrow over (r)} and is synonymous with the coding complexity of the data given the estimate of the channel. From this perspective it is proportional to the length of the best code of the residuals of the data predicted by the adaptive filter estimate and source signal. Similarly the sum of the second and third terms is proportional to the information necessary to specify this particular channel estimate {right arrow over (ĥ)} B  from the prior assignment to the resolution associated with the posterior variance.
 
   A popular and heuristically simple alternative is Akaike&#39;s Information Criteria (AIC) (for a reasonable explanation of this ad-hoc measure see Kay). The associated penalty cost in this case is simply one half the number of parameters of the model
 
 {circumflex over (B)}   AIC=ar gmin   J   {Q   AIC ( J )}  (29)
 
where
 
   
     
       
         
           
             
               
                 
                   
                     Q 
                     AIC 
                   
                   ⁡ 
                   
                     ( 
                     J 
                     ) 
                   
                 
                 = 
                 
                   
                     
                       - 
                       log 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       p 
                       ⁡ 
                       
                         ( 
                         
                           
                             r 
                             → 
                           
                           | 
                           
                             
                               
                                 w 
                                 → 
                               
                               ^ 
                             
                             J 
                           
                         
                         ) 
                       
                     
                   
                   + 
                   
                     
                       3 
                       ⁢ 
                       JK 
                     
                     2 
                   
                   + 
                   
                     
                       2 
                       
                         J 
                         - 
                         1 
                       
                     
                     ⁢ 
                     F 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         τ 
                         max 
                       
                       . 
                     
                   
                 
               
             
             
               
                 ( 
                 30 
                 ) 
               
             
           
         
       
     
   
   These in-scale algorithms were tested on a diverse set of simulated multipath channels to determine their relative performance in jointly estimating the channel operator and its Doppler spread. 
   The considered test channels are of the form 
                   H   ⁡     (     t   ,   τ     )       =       ∑     m   =   1       M   p       ⁢         α   m     ⁡     (   t   )       ⁢       χ   m     ⁡     (     τ   -       d   m     ⁡     (   t   )         )                   (   31   )               
having M p  independent component paths. Each path has an associated arrival spreading function x m (•) that is time invariant. Doppler spread is induced by both the path gain processes α m (t) and the path delay processes d m (t). Arrival path delay times d m (t) and amplitudes α m (t) are modeled as correlated Gaussian variates with respective covariances c d     m   (τ)=c d     m   (0)(1−τ/T d     m   ) for |t|&lt;T d     m   , 0 otherwise and c α     m   (t)=c α     m   (0)(1−t/T α     m   for |t|&lt;T α     m   , 0 otherwise. A realization of such a channel is then generated given the 3M d  parameters: the weight functions χ m (·), and coherence times T d     m    and T α     m   . An associated Doppler bandwidth for each of these path processes is termed B α     m   =2/T α     m    and B d     m   =2/T d     m   . For each test case the Doppler spread is termed the maximum of these constituent Doppler bandwidths. Varying the correlation times associated with the processes d m &#39;s and α m &#39;s as well as the spreading shapes χ m &#39;s allows for the simulation of processes that are diverse and realistic.
 
   Four channel operators were tested, the number of paths and the maximal Doppler spread for each of the cases is listed in Table 2. 
   
     
       
             
           
             
             
             
             
           
         
             
               TABLE 2 
             
           
           
             
                 
             
             
               Test cases for time-varying channel parameters 
             
           
        
         
             
                 
               case 
               M d   
               B (Doppler spread) 
             
             
                 
                 
             
             
                 
               LTI 
               4 
               0 Hz 
             
             
                 
               TV α&#39;s 
               4 
               1 Hz 
             
             
                 
               TV delay&#39;s 
               4 
               4 Hz 
             
             
                 
               TV delay&#39;s &amp; α&#39;s 
               4 
               4 Hz 
             
             
                 
                 
             
           
        
       
     
   
   Some channels share a common feature of sparsity in that while each scatterer path exhibits Gaussian uncorrelated scattering each delay lag does not. Thus the channels are not wide sense stationary-uncorrelated scattering (WSS-US). 
   Gaussian source {right arrow over (s)} and additive noise signals of the same spectra where used in the simulations. The source signal correlation c {right arrow over (s)} (τ)=E[s(t)s(t+τ)]=ζ 2 (1−2τF) for |τ|&lt;2F, 0 otherwise. c {right arrow over (n)} (τ)=σ 2 c {right arrow over (s)} (τ)/ζ 2 . The channels were simulated for durations of approximately T=8 seconds and a bandwidth F=2 k Hz. 
   In these simulations C B  and S B  were used as approximations to C B   ξ  and S B   ξ  in the CG algorithm for the computational savings associated with its tight banded structure. In practice the number of C G iterations to compute {right arrow over ({tilde over (h)} B   a  starting from (D ξ   1 {circle around (×)}I 2     K   ){tilde over ({right arrow over (h)} B/2   a  hB/2 to a chosen precision of 0.1×P×δ J   2  is approximately 5. Further iterations are used to compute the MLE. The algorithm was implemented with Daubechies orthogonal wavelets order 3 (ψ=D(3)) in delay and ξ=D(5) in time for U J  under the folded wavelet assumption. The interpolation operator D ξ   −q  of (3) in this simulation was implemented as a windowed sinc for the simplicity associated with its linear phase and q was set at 4 so that for a given hypothesized J the estimate was assumed LTI over durations of T/2 J+4 . 
   The normalized mean square error for the channel estimates 
                    e   h          2   2     =     [           (           h   →     ^     B     -       h   →     B       )     ′     ⁢     (           h   →     ^     B     -       h   →     B       )             h   →     B   ′     ⁢       h   →     B         ]           
is computed for each of the estimators; the approximate MLE
 
             h     →   ~       B   a         
(aMLE), the MLE, {right arrow over (h)} B   mle , the posterior mean
 
             h     →   ^       B         
(PM) and the approximate MAP.
 
   Granular noise is present in the MLE and represents the cost incured by the use of this unbiased estimator. With the sparsity prior the PM is able to provide a more regular estimate that exhibits better MSE performance. However this performance comes at a cost of bias in the PM estimate. This bias is pronounced in the peaks of the channel response estimate. Indeed it is easily verified by considering near L ∞  loss functions. Tests were conducted with ∥e h ∥ p  p=8 to approximate L ∞  loss and demonstrate that the MLE often outperforms the PM even for these sparse test cases. 
   To determine the ability of the in-scale approach to accurately determine Doppler spread the algorithm can be tested on an LTI operator. One important attribute of the in-scale approach is that LTI filters can be recognized even at quite low SNR. The in-scale algorithm can identify LTI channels with only the computation of 2 iterations in the scale domain. The first iteration is the LTI estimate. The second is the scale 1 iteration from which the LTI is compared favorably and the in-scale iteration is ceased. This is in stark contrast to in-time recursions where computational resources are distributed over time without respect to the actual coherence time of the channel process. Thus for finite duration signaling the in-scale approach has a distinct advantage. The Doppler complexity measure is a minimum for the J=0 th  scale estimate for all three criteria (AIC, LAP, MAP) and does not fall below this at greater resolution estimates. Here each of the Doppler costs are normalized by E[log p({right arrow over (r)}|{right arrow over (h)} B )]=−(1+log σ 2 )FT/2. It is further noted that for sparse channels in delay the mixture normal model is an effective sparsing model affording a 2 dB improvement in performance over the MLE estimate at the 0 th  scale estimate. The PM estimate exhibits greater gains at higher scales demonstrating the robustness of the posterior mean and MAP estimates against variance of the Doppler spread rule. 
   Monte Carlo simulations reveal slightly greater variability of the MAP stopping rule over that of the AIC or LAP. We can demonstrate that a model of the time varying filter as a process in scale yields estimators that provide joint estimation of channel response and coherence time. Coherence time is modeled in this framework as model order selection. 
   An MLE, approximate MLE, approximate MAP and posterior mean were derived along with AIC, Laplace approximation and MAP Doppler estimation procedures based on empirical Bayes assumptions. These algorithms were tested and compared with one and other on a number of simulated time varying channels. As expected the posterior mean demonstrates improvement in performance over the MLE and approximate MAP for time varying channels relative to mean square error. 
   All of the estimators presented for in-scale processing rely on the conjugate gradient algorithm to compute approximate MLE estimates from the previous lower scale estimate. From the MLE and conditional variance along with a sparsity prior the posterior mean is a simple adaptive thresholding of wavelet coefficents. In this framework computation is distributed over scale in the estimator and this contrasts with in-time Kalman methods that distribute computation in time for an assumed coherence time. Joint estimation of channel parameters and coherence time is fundamental to in-scale estimation as it provides a stopping rule at which additional computation is unlikely to provide an improved estimate. Lastly covariance estimates of the time varying model are provided and approximate high probability regions are shown to be easily accessible. These present useful computationally efficient empirical Bayes lower bounds on certainty bands of the channel operator. 
   The implications for multisensor and multichannel estimation are apparent. Computational resources can be automatically and seemlessly focused on system nodes that exhibit greater Doppler bandwidth. This implies two advantages, first channels that are more coherent will not waste computational resources. Secondly performance degradation is limited since overfitting of data associated with overly complex models is less likely. 
   Important issues relating to computation and performance relative to in-time recursions must be addressed. Future work should focus on comparison of in-scale recursions with in-time recursions (e.g. Kalman, RLS) to give bounds on signal duration, channel sparsity and coherence time measures where the in-scale regime is to be favored. The in-scale approach could be broadened to include other scale stepping increments rather than the classic octave band partition presented here. Mixed scale-time methods can also be envisioned for block processing over time. 
   The following matrix and Kronecker properties are useful. 
   Properties: 
   
       
       
         
           1. [A H B] st =(B′{circle around (×)}A)[H] st    
           2. (A C{circle around (×)}B D)=(A{circle around (×)}B)(C{circle around (×)}D) 
           3. (A{circle around (·)}C){circle around (×)}(B{circle around (·)}D)=(A{circle around (×)}B){circle around (·)}((C{circle around (×)}D) where A M×N C M×N  and B L×K  D L×K . 
           4. (A{circle around (×)}B) −1 =A −1 {circle around (×)}B −1 , and (A{circle around (×)}B)′=A′{circle around (×)}B′ 
           5. diag(A B)=[A{circle around (·)}B′]{right arrow over (1)} M  where A M×N  and B N×M    
           6. diag(A{circle around (·)}B)=diag(A){circle around (·)}diag(B) where A N×N  and B N×N  and B L×LM/N    
           7. A{right arrow over (b)}=A Diag({right arrow over (b)}){right arrow over (1)} N    
           8. A Diag({right arrow over (b)})=A{circle around (·)}({right arrow over (b)}′{circle around (×)}{right arrow over (1)} M ) where A M×N  and {right arrow over (b)} N×1    
         
       
     
  
   The following result is useful when computing the posterior variance of a time varying operator
 
for: A K×L , C L×K , B M×N , D N×M , {right arrow over (v)} LN×1 
 
if:  V =( A{circle around (×)}B )Diag( {right arrow over (v)} )( C{circle around (×)}D )
 
then: diag( V )=[( A{circle around (·)}C ′){circumflex over (×)}( B{circle around (·)}D ′)] {right arrow over (v)}   (32)
 
Proof:
 
   
     
       
         
           
             
               V 
             
             
               = 
             
             
               
                 
                   ( 
                   
                     
                       ( 
                       
                         A 
                         ⊗ 
                         B 
                       
                       ) 
                     
                     ⊙ 
                     
                       ( 
                       
                         
                           v 
                           ′ 
                         
                         ⊗ 
                         
                           
                             1 
                             → 
                           
                           KM 
                         
                       
                       ) 
                     
                   
                   ) 
                 
                 ⁢ 
                 
                   ( 
                   
                     C 
                     ⊗ 
                     D 
                   
                   ) 
                 
               
             
             
               
                 property 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 8 
               
             
           
           
             
               
                 diag 
                 ⁡ 
                 
                   ( 
                   V 
                   ) 
                 
               
             
             
               = 
             
             
               
                 
                   [ 
                   
                     
                       ( 
                       
                         A 
                         ⊗ 
                         B 
                       
                       ) 
                     
                     ⊙ 
                     
                       ( 
                       
                         
                           v 
                           ′ 
                         
                         ⊗ 
                         
                           
                             1 
                             → 
                           
                           KM 
                         
                       
                       ) 
                     
                     ⊙ 
                     
                       ( 
                       
                         
                           C 
                           ′ 
                         
                         ⊗ 
                         
                           D 
                           ′ 
                         
                       
                       ) 
                     
                   
                   ] 
                 
                 ⁢ 
                 
                   
                     1 
                     → 
                   
                   
                     L 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     N 
                   
                 
               
             
             
               
                 property 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 5 
               
             
           
           
             
               
                   
               
             
             
               = 
             
             
               
                 
                   [ 
                   
                     
                       ( 
                       
                         
                           ( 
                           
                             A 
                             ⊗ 
                             B 
                           
                           ) 
                         
                         ⊙ 
                         
                           ( 
                           
                             
                               C 
                               ′ 
                             
                             ⊗ 
                             
                               D 
                               ′ 
                             
                           
                           ) 
                         
                       
                       ) 
                     
                     ⁢ 
                     
                       Diag 
                       ⁡ 
                       
                         ( 
                         
                           v 
                           → 
                         
                         ) 
                       
                     
                   
                   ] 
                 
                 ⁢ 
                 
                   
                     1 
                     → 
                   
                   
                     L 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     N 
                   
                 
               
             
             
               
                 
                   property 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   8 
                 
                 ⁢ 
                 
                     
                 
               
             
           
           
             
               
                   
               
             
             
               = 
             
             
               
                 
                   [ 
                   
                     
                       ( 
                       
                         A 
                         ⊗ 
                         B 
                       
                       ) 
                     
                     ⊙ 
                     
                       ( 
                       
                         
                           C 
                           ′ 
                         
                         ⊗ 
                         
                           D 
                           ′ 
                         
                       
                       ) 
                     
                   
                   ] 
                 
                 ⁢ 
                 
                   v 
                   → 
                 
               
             
             
               
                 property 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 7 
               
             
           
           
             
               
                   
               
             
             
               = 
             
             
               
                 
                   [ 
                   
                     
                       ( 
                       
                         A 
                         ⊙ 
                         
                           C 
                           ′ 
                         
                       
                       ) 
                     
                     ⊗ 
                     
                       ( 
                       
                         B 
                         ⊙ 
                         
                           D 
                           ′ 
                         
                       
                       ) 
                     
                   
                   ] 
                 
                 ⁢ 
                 
                   v 
                   → 
                 
               
             
             
               
                 property 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 3 
               
             
           
         
       
     
   
   The following form of the posterior distribution is useful. The prior and likelihood on wavelet coefficients are w j,k,* |π j,k , λ j,k ˜π j,k N w     j,k*   (0, λ j,k   2 )+(1−π j,k )N(0, ∈ j,k   2 ) and {tilde over (w)} j,k,* |w j,k* ˜N {tilde over (w)}     j,k,*   (w j,k,* , δ j,k   2 ) respectively. The marginal density of empirical coefficients is then 
             p   ⁡     (         w   ~       j   ,   k   ,   *       |       ϕ   →     J       )       =       ∫       p   ⁡     (         w   ~       j   ,   k   ,   *       |     w     j   ,   k   ,   *         )       ⁢     p   ⁡     (         w     j   ,   k   ,   *       |     π     j   ,   k         ,     λ     j   ,   k       ,     ε     j   ,   k         )       ⁢     ⅆ     w     j   ,   k   ,   *             =         π     j   ,   k       ⁢       ??       w   ~       j   ,   k   ,   *         ⁡     (     0   ,       λ     j   ,   k     2     +     δ     j   ,   k     2         )         +       (     1   -     π     j   ,   k         )     ⁢         ??         w   ~       j   ,   k       ,   *       ⁡     (     0   ,       ε     j   ,   k     2     +     δ     j   ,   k     2         )       .     
     ⁢   Since                             ??     w     j   ,   k   ,   *         ⁡     (         w   ~       j   ,   k   ,   *       ,     δ     j   ,   k     2       )       ⁢       ??       w     j   ,   k       ,   *       ⁡     (     0   ,     λ     j   ,   k     2       )         =         ??     w     j   ,   k   ,   *         ⁡     (         γ   ⁡     (     λ     j   ,   k       )       ⁢       w   ~       j   ,   k   ,   *         ,       γ   ⁡     (     λ     j   ,   k       )       ⁢     δ     j   ,   k     2         )       ⁢       ??       w   ~       j   ,   k   ,   *         ⁡     (     0   ,       λ     j   ,   k     2     +     δ     j   ,   k     2         )               
where γ(λ)=λ 2 /(λ 2 +δ 2 ) is akin to a Wiener gain, a direct application of Bayes theorem p(w|{tilde over (w)})=p({tilde over (w)}|w)p(w)/p({tilde over (w)}) confirms the posterior density to be
 
   
     
       
         
           
             
               
                 
                   
                     
                       w 
                       
                         j 
                         , 
                         k 
                         , 
                         * 
                       
                     
                     | 
                     
                       
                         w 
                         ~ 
                       
                       
                         j 
                         , 
                         k 
                         , 
                         * 
                       
                     
                   
                   , 
                   
                     
                       
                         
                           
                             ϕ 
                             → 
                           
                           J 
                         
                         ~ 
                         
                           Π 
                           ⁡ 
                           
                             ( 
                             
                               
                                 w 
                                 ~ 
                               
                               
                                 j 
                                 , 
                                 k 
                                 , 
                                 * 
                               
                             
                             ) 
                           
                         
                       
                       ⁢ 
                       
                         
                           ?? 
                           
                             w 
                             
                               j 
                               , 
                               k 
                               , 
                               * 
                             
                           
                         
                         ⁡ 
                         
                           ( 
                           
                             
                               
                                 γ 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     λ 
                                     
                                       j 
                                       , 
                                       k 
                                     
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               
                                 
                                   w 
                                   ~ 
                                 
                                 
                                   j 
                                   , 
                                   k 
                                   , 
                                   * 
                                 
                               
                             
                             , 
                             
                               
                                 γ 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     λ 
                                     
                                       j 
                                       , 
                                       k 
                                     
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               
                                 δ 
                                 
                                   j 
                                   , 
                                   k 
                                 
                                 2 
                               
                             
                           
                           ) 
                         
                       
                     
                     + 
                     
                       
                         ( 
                         
                           1 
                           - 
                           
                             Π 
                             ⁡ 
                             
                               ( 
                               
                                 
                                   w 
                                   ~ 
                                 
                                 
                                   j 
                                   , 
                                   k 
                                   , 
                                   * 
                                 
                               
                               ) 
                             
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         
                           ?? 
                           
                             w 
                             
                               j 
                               , 
                               k 
                               , 
                               * 
                             
                           
                         
                         ⁡ 
                         
                           ( 
                           
                             
                               
                                 γ 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     ε 
                                     
                                       j 
                                       , 
                                       k 
                                     
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               
                                 
                                   w 
                                   ~ 
                                 
                                 
                                   j 
                                   , 
                                   k 
                                   , 
                                   * 
                                 
                               
                             
                             , 
                             
                               
                                 γ 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     ε 
                                     
                                       j 
                                       , 
                                       k 
                                     
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               
                                 δ 
                                 
                                   j 
                                   , 
                                   k 
                                 
                                 2 
                               
                             
                           
                           ) 
                         
                       
                     
                   
                 
                 ⁢ 
                 
                   
 
                 
                 ⁢ 
                 where 
                 ⁢ 
                 
                   
 
                 
                 ⁢ 
                 
                   
                     Π 
                     ⁡ 
                     
                       ( 
                       
                         
                           w 
                           ~ 
                         
                         
                           j 
                           , 
                           k 
                           , 
                           * 
                         
                       
                       ) 
                     
                   
                   = 
                   
                     
                       [ 
                       
                         1 
                         + 
                         
                           
                             1 
                             - 
                             
                               
                                 π 
                                 
                                   j 
                                   , 
                                   k 
                                 
                               
                               ⁢ 
                               
                                 
                                   ?? 
                                   
                                     
                                       w 
                                       ~ 
                                     
                                     
                                       j 
                                       , 
                                       k 
                                       , 
                                       * 
                                     
                                   
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     0 
                                     , 
                                     
                                       
                                         ε 
                                         
                                           j 
                                           , 
                                           l 
                                         
                                         2 
                                       
                                       + 
                                       
                                         δ 
                                         
                                           j 
                                           , 
                                           k 
                                         
                                         2 
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                           
                             
                               π 
                               
                                 j 
                                 , 
                                 k 
                               
                             
                             ⁢ 
                             
                               
                                 ?? 
                                 
                                   
                                     w 
                                     ~ 
                                   
                                   
                                     j 
                                     , 
                                     k 
                                     , 
                                     * 
                                   
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   0 
                                   , 
                                   
                                     
                                       λ 
                                       
                                         j 
                                         , 
                                         k 
                                       
                                       2 
                                     
                                     + 
                                     
                                       δ 
                                       
                                         j 
                                         , 
                                         k 
                                       
                                       2 
                                     
                                   
                                 
                                 ) 
                               
                             
                           
                         
                       
                       ] 
                     
                     
                       - 
                       1 
                     
                   
                 
               
             
             
               
                 ( 
                 33 
                 ) 
               
             
           
         
       
     
   
   The following form of the variance is useful. For the variance, start with var[w j,k,* |{tilde over (w)} j,k,* , {tilde over (φ)} J ]=E[w j,k,*   2 |{tilde over (w)} j,k,* , {right arrow over (φ)} J ]−E 2 [w j,k,* |{right arrow over (w)} j,k,* , {right arrow over (φ)} J ]. For simplicity denote the modulus square of a coefficient by x 2 =x′x. From the posterior density the second moment E[w j,k,*,   2 |{tilde over (w)} j,k,*, {right arrow over (φ)}   J ] is easily evaluated as 
             E   ⁡     [         w     j   ,   k   ,   *     2     |       w   ~       j   ,   k   ,   *         ,       ϕ   →     J       ]       =         Π   ⁡     (       w   ~       j   ,   k   ,   *       )       ⁢     (         γ   ⁡     (   λ   )       ⁢     δ     j   ,   k     2       +         γ   2     ⁡     (   λ   )       ⁢       w   ~       j   ,   k   ,   *     2         )       +       (     1   -     Π   ⁡     (       w   ~       j   ,   k   ,   *       )         )     ⁢     (         γ   ⁡     (   ε   )       ⁢     δ     j   ,   k     2       +         γ   2     ⁡     (   ε   )       ⁢     δ     j   ,   k     2       +         γ   2     ⁡     (   ε   )       ⁢       w   ~       j   ,   k   ,   *     2         )               
and of course E 2 [w j,k,* |{right arrow over (w)} j,k,* , {right arrow over (φ)} J ] is taken from (17). The result after simplification is
 
                   var   ⁡     [         w     j   ,   k   ,   *       |       w   ~       j   ,   k   ,   *         ,       ϕ   →     J       ]       =         [         Π   ⁡     (       w   ~       j   ,   k   ,   *       )       ⁢     γ   ⁡     (   λ   )         +       (     1   -     Π   ⁡     (       w   ~       j   ,   k   ,   *       )         )     ⁢     γ   ⁡     (   ε   )           ]     ⁢     δ     j   ,   k     2       +       Π   ⁡     (       w   ~       j   ,   k   ,   *       )       ⁢         (     1   -     Π   ⁡     (       w   ~       j   ,   k   ,   *       )         )     ⁡     [       γ   ⁡     (   λ   )       -     γ   ⁡     (   ε   )         ]       2     ⁢         w   ~       j   ,   k   ,   *     2     .                 (   34   )               
Since lim x→0  Π(x)=0 and lim x→∞  Π(x)=1, it is easily confirmed that for λ&gt;&gt;δ and ∈/δ&lt;&lt;1
             lim var[w|{tilde over (w)}]=δ 2            lim var[w|{tilde over (w)}]=∈ 2   (35) 
while the maximum variance attainable is approximately var[w|{tilde over (w)}={tilde over (w)} p ]≈δ 2 (2+log(λ 2 /δ 2 π 2 ))/4 where

   
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           w 
                           ~ 
                         
                         p 
                         2 
                       
                       = 
                       
                         arg 
                         ⁢ 
                         
                           { 
                           
                             
                               Π 
                               ⁡ 
                               
                                 ( 
                                 
                                   w 
                                   ~ 
                                 
                                 ) 
                               
                             
                             = 
                             
                               1 
                               / 
                               2 
                             
                           
                           } 
                         
                       
                     
                   
                 
                 
                   
                     
                       = 
                       
                         
                           
                             
                               ( 
                               
                                 
                                   λ 
                                   2 
                                 
                                 + 
                                 
                                   δ 
                                   2 
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               ( 
                               
                                 
                                   ε 
                                   2 
                                 
                                 + 
                                 
                                   δ 
                                   2 
                                 
                               
                               ) 
                             
                           
                           
                             
                               λ 
                               2 
                             
                             - 
                             
                               ε 
                               2 
                             
                           
                         
                         ⁢ 
                         
                           log 
                           ⁡ 
                           
                             ( 
                             
                               
                                 
                                   ( 
                                   
                                     
                                       λ 
                                       2 
                                     
                                     + 
                                     
                                       δ 
                                       2 
                                     
                                   
                                   ) 
                                 
                                 
                                   ( 
                                   
                                     
                                       ε 
                                       2 
                                     
                                     + 
                                     
                                       δ 
                                       2 
                                     
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               
                                 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       π 
                                     
                                     ) 
                                   
                                   2 
                                 
                                 
                                   π 
                                   2 
                                 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                 
               
             
             
               
                 ( 
                 36 
                 )