Minimum Mean-Squared Error Reconstruction for Generalized ...
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
3237
On the other hand, we have shown that B (^b) is idempotent, that is, B (^b) = B 2 (^b), which means that it only has eigenvalues 0 and 1 and > 0. Combining the results, the zero-phase property hence [B (^b)]
Minimum Mean-Squared Error Reconstruction for Generalized Undersampling of Cyclostationary Processes
REFERENCES
Ryan S. Prendergast, Student Member, IEEE, and Truong Q. Nguyen, Fellow, IEEE
m;m
of the denominator of (13) is proved.
[1] S. M. Kay, Modern Spectral Estimation : Theory and Application. Englewood Cliffs, NJ: Prentice-Hall, 1988. [2] S. Kay and R. Nekovei, “An efficient two-dimensional frequency estimator,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 38, no. 10, pp. 1807–1809, Oct. 1990. [3] Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” IEEE Trans. Signal Process., vol. 40, no. 9, pp. 2267–2280, Sep. 1992. [4] M. P. Clark and L. L. Scharf, “Two-dimensional modal analysis based on maximum likelihood,” IEEE Trans. Signal Process., vol. 42, no. 6, pp. 1443–1452, Jun. 1994. [5] C. R. Rao, L. Zhao, and B. Zhou, “Maximum likelihood estimation of 2-D superimposed exponential signals,” IEEE Trans. Signal Process., vol. 42, no. 7, pp. 1795–1802, Jul. 1994. [6] A. Sourice, G. Plantier, and J.-L. Saumet, “Two-dimensional frequency estimation using autocorrelation phase fitting,” in Proc. Int. Conf. Acoustics, Speech, Signal Processing (ICASSP), Hong Kong, Apr. 2003, vol. 3, pp. 445–448. [7] P. Stoica and R. Moses, Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice-Hall, 2005. [8] S. A. Tretter, “Estimating the frequency of a noisy sinusoid by linear regression,” IEEE Trans. Inf. Theory, vol. 31, pp. 832–835, Nov. 1985. [9] G. W. Lank, I. S. Reed, and G. E. Pollon, “A semicoherent detection and Doppler estimation statistic,” IEEE Trans. Aerosp. Electron. Syst., vol. 9, pp. 151–165, Mar. 1973. [10] S. Kay, “A fast and accurate single frequency estimator,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 12, pp. 1987–1990, Dec. 1989. [11] H. Chien, “2-D estimation from AR models,” Ph.D. dissertation, Elect. Eng. Dept., Univ. of Rhode Island, Kingston, RI, 1981. [12] R. Kumaresan, L. L. Scharf, and A. K. Shaw, “An algorithm for polezero modeling and spectral analysis,” IEEE Trans. Acoust., Speech, Signal Process., vol. 34, no. 6, pp. 637–640, Jun. 1986. [13] Y. Bresler and A. Macovski, “Exact maximum likelihood parameter estimation of superimposed exponential signals in noise,” IEEE Trans. Acoust., Speech, Signal Process., vol. 34, no. 10, pp. 1081–1089, Oct. 1986. [14] G. C. Goodwin and R. L. Payne, Dynamic System Identification: Experiment Design and Data Analysis. New York: Academic, 1977, pp. 46–50. [15] H. C. So and K. W. Chan, “A generalized weighted linear predictor frequency estimation approach for a complex sinusoid,” IEEE Trans. Signal Processing, vol. 54, no. 4, pp. 1304–1315, Apr. 2006. [16] B. Völcker and P. Händel, “Frequency estimation from proper sets of correlations,” IEEE Trans. Signal Process., vol. 50, no. 4, pp. 791–802, Apr. 2002. [17] M. L. Fowler and J. A. Johnson, “Extending the threshold and frequency range for phase-based frequency estimation,” IEEE Trans. Signal Process., vol. 47, no. 10, pp. 2857–2863, Oct. 1999. [18] H. Lutkepohl, Handbook of Matrices. New York: Wiley, 1996.
Abstract—The generalized sampling problem is considered using a filter bank model with a stochastic process framework. A signal reconstruction solution minimizing the time-averaged mean-squared error is found using a filter bank optimization technique for cyclostationary signals. The use of a stochastic model approach to find a minimized error solution rather than a perfect reconstruction solution allows consideration of sampling problems that deterministic approaches cannot solve. For instance, a particular sampling density is not required, allowing reconstruction optimization in cases of undersampling. An extension is also presented for optimal reconstruction in the presence of additive noise and interference. The result is a generalized sampling model which can be used for an extensive variety of scenarios, combining aspects of undersampling, multiple-input multiple-output (MIMO) sampling, and periodic nonuniform sampling. Index Terms—Cyclostationary processes, generalized sampling, minimum mean-square error (MMSE) filtering, undersampling.
I. INTRODUCTION This correspondence examines a modified approach to the generalized sampling problem using a filter bank model with a discrete-time cyclostationary input. Traditional implementations of Papoulis’ generalized sampling expansion [1] have considered sampling of several linearly filtered versions of a bandlimited continuous-time signal at a collective Nyquist rate, from which perfect reconstruction (PR) can be achieved. Since the sampled signal is perfectly bandlimited, an equivalent discrete-time model is found using the multirate filter bank, and solutions can be determined based on established PR filter bank theory [2]. While previous generalized sampling approaches have mainly dealt with deterministic signals of known finite spectral support, this correspondence will consider stochastic signals with known spectral densities. The solution does not require a minimum sampling density to determine a signal reconstruction. Instead, the optimal reconstruction is found for a given sampling process. This provides a distinct advantage over the traditional approach, allowing for optimal reconstruction solutions in the event of nonuniform undersampling and in the presence of interference. As such, the result is directly applicable to undersampling problems such as image super-resolution, and may also be of interest for sampling rate reduction of analog-to-digital conversion and equalization systems. Much previous literature in the area of filter banks have shown that the spectral characteristics of stochastic signals can be taken advantage of to determine optimal filtering structures for a variety of applications [3]–[9]. The generalized sampling problem requires a different approach than many other filter bank optimization problems since it must incorporate a fixed model based on the sampling process. In other cases, typically both the analysis and synthesis filter banks (see Fig. 1) are jointly optimized. The filter bank model for a sampling problem Manuscript received March 4, 2005; revised August 28, 2005. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. H. Boelcskei. This research has been supported by a grant from Skyworks, Inc. and UC Dimi. The authors are with the Department of Electrical and Computer Engineering, University of California, San Diego, La Jolla, CA 92093 USA (e-mail: [email protected]; [email protected]). Digital Object Identifier 10.1109/TSP.2006.877649
1053-587X/$20.00 © 2006 IEEE
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
3238
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
x(nTM ) with a sufficiently short model period TM
Fig. 1. General multirate filter bank structure.
uses the analysis bank to represent a physical sampling process and any enhancement or processing is confined to the synthesis bank. As such, the general problem to be solved can be stated as follows. Given a particular generalized sampling process (represented through analysis filters H0 (ej! ); . . . ; HC 01 (ej! )), determine a linear processing technique (through synthesis filters F0 (ej! ); . . . ; FC 01 (ej! )) that is optimal for some criterion. The presented approach determines an optimal solution through minimization of the mean of je[n]j2 = jx[n] 0 y [n]j2 , the squared error. Since a stationary signal model and linear synthesis filters are used, the result is a Wiener filtering solution for the generalized sampling problem. The cyclostationary filter bank analysis techniques developed by Ohno and Sakai [5] are used to determine a frequency domain solution for the synthesis filters. For the purposes of this correspondence, the sampled signal will be considered using a discrete-time model. The limitations associated with the use of a discrete-time model for real-world signals can be dismissed since the resolution of this model can be arbitrarily increased. A considerable amount of work has been done using filter bank models of generalized sampling since Papoulis’ result. One problem of interest is a purely periodic sampling problem, which is considered through Fig. 1 when the analysis filters are simple delays. This has been considered for signals of finite spectral support in a number of different cases [10]–[18]. Primary interest is in periodic nonuniform sampling, since results for uniform sampling are well established. The use of periodic sampling patterns of sufficient average rate (equal to the Nyquist rate for lowpass signals) allows PR solutions to be found, although practical implementations can only approximate PR. Some modifications to the standard problem include undersampling [19], multi-input problems [20], [21], and reconstruction using alternative (to finite spectral support) bases [22]–[24]. This correspondence’s method is capable of dealing with a wide variety of sampling problems using a minimized mean-square error (MSE) approach. In some instances, PR solutions can be found, but optimal reconstruction solutions are still found when PR cannot be achieved. The organization of this correspondence is as follows. Section II provides background on the use of a discrete-time framework for generalized sampling problems and on filter bank analysis of cyclostationary signals. Section III presents the general optimal solution for the sampling problem with some related discussion and an extension allowing consideration of multiple inputs and additive noise. Concluding remarks are provided in Section IV, along with a brief mention of problem extensions and future work. II. PRELIMINARIES
enables representation without loss of information. In order to maintain the periodicity of the sampling process in a filter bank representation, it will be required that the sampling process period TP be divisible by TM , or DTM = TP , where D is a positive integer corresponding to the degree of decimation in the filter bank model. For a given TP , the decimation factor D is herein assumed to be sufficiently large so that x[n] = x(nTP =D) provides a lossless representation of x(t). Through this assumption, strictly discrete-time signal models can be considered. It should be noted that use of the discrete-time representation will generally not limit the modelling of a physical periodic sampling process. Referring to Fig. 1, each branch of the analysis filter bank models a physical uniform sampling with a unique delay and distortion function. The collection of these C samplers forms the periodic sampling process. In absence of distortion, each sampler is represented through the delay term Hi (ej! ) = e0j!d : Allowing the use of fractional delays enables consideration of any periodic sampling process. It is known that fractional delays can only be perfectly represented through filters of infinite length in the time domain [25]. However, as is done in this correspondence, ideal fractional delay filters are easily considered in the frequency domain over the limited range j! j < . B. Filter Bank Analysis of Cyclostationary Signals This subsection introduces definitions based on the analysis developed by Ohno and Sakai in [5] for the purpose of minimizing a filter bank’s time-averaged MSE (TAMSE) when a subband is lost. The approach considers a discrete-time cyclostationary input even though primary interest will be for wide-sense stationary (WSS) inputs. The reason for this is that even when the filter bank is provided a stationary input, it will generally produce an output that is cyclostationary with period D (corresponding to the order of decimation). Because of this, the error between the input and output signals is cyclostationary. For generalization purposes, the input is considered so as well (since stationarity is a subset of cyclostationarity). Specifically, x[n] is wide-sense cyclostationary with period D , or WSCS(D), meaning its first- and second-order expectations are periodically stationary with period D . This is expressed as
E [x[n]] = E [x[n + kD]]
and Rxx [n; m] = Rxx [n + kD; m + kD ]
(1)
for every integer k , where Rxx [n; m] = E [x[n]x 3 [m]]. Herein, a zeromean process is considered. The cyclic correlation function is defined as
[ u] = Rxx
D01
R [k + u; k]e0j 2 k D k=0 xx 1
(2)
for = n=D and integer n (the notation used in [5] is based on that of [26]). Since the exponential term is periodic with , (2) only needs to be determined for 0 (D 0 1)=D (or, n = 0; 1; . . . ; D 0 1) to be defined for all valid . The discrete-time Fourier transforms of the cyclic correlation functions produce the cyclic spectral densities (CSDs) through
A. Discrete-Time Framework for Sampling Model This section describes the discrete-time framework used to model a periodically sampled continuous-time signal. Provided this continuous-time signal x(t) is bandlimited, a discrete-time model x[n] =
(ej! ) = Sxx
1 u=01
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
[u]e0j!u : Rxx
(3)
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
Through (2) and (3), the time-averaged variance of a WSCS(D) process is defined as 1
x2 = Nlim !1
N 01
N n=0
= Rxx [0] = 0
E [x[n]x3 [n]] =
1
2
0
1
D01
R [n; n] D n=0 xx 2 0 j! Sxx (e )d!:
3239
the minimization of (9) with respect to the synthesis AC matrix can be stated as
F (4)
(
p;q)
=
(p0q )=D j! p Sxx (e W )
over
j!j D
(5)
p; q = 0; . . . ; D 0 1, using the standard definition of e0j (2=D) . This representation is only valid for j!j =D ; however, this is sufficient to completely represent all D CSDs of (3). [u] = 0 for all valid 6= 0 if and only if x[n] is WSS, it is Since Rxx easily verified that a diagonal Sxx (ej! ) is a necessary and sufficient for
W
111
H0 (ej! )
HAC (ej! ) =
.. .
..
H0 (ej! W (D01) )
HC 01 (ej! ) .. .
.
(6)
1 1 1 HC01 (ej! W (D01) ) for the analysis bank and an equivalent representation FAC (ej! ) for the synthesis bank. Like the CSD matrix, these AC representations are also only considered over j! j =D but completely and uniquely describe the analysis and synthesis filters in the frequency domain. Defining the matrix product
P(ej! ) = D1 FAC (ej! )HTAC (ej! )
)
tr(
See (ej! ))
See (ej! ) = Sxx (ej! ) 0 Syx (ej! ) 0 Sxy (ej! ) + Syy (ej! ) j! j! j! H j! H j! = Sxx (e ) 0 P(e )Sxx (e ) 0 Sxx (e )P (e ) j! j! H j! (8) + P(e )Sxx (e )P (e ) which represents the error as a function of the filter bank and the input signal statistics.
=
2
2 1
2
0
0 j! See (e )d! =
=D
0=D
tr(
1
2
See (ej! ))d!:
D01
=D
k=0
0=D
(10)
FAC;opt (ej! ) = DQH (ej! )R01 (ej! )
(11)
T (ej! )S (ej! ) Q(ej! ) = HAC xx
(12)
T (ej! )SH (ej! )H3 (ej! ): R(ej! ) = HAC xx AC
(13)
and
An expression for the minimized TAMSE is found by applying the solution (11) to (8), resulting in the optimal error CSD matrix
See;opt (ej! ) = Sxx (ej! ) 0 QH (ej! )R0H (ej! )Q(ej! ) H j! QH (ej! )FAC ;opt (e ) : j! = Sxx (e ) 0 D
(14)
As through (9), the optimal TAMSE e;opt can be determined from this expression. The remainder of this section provides further analysis and extensions on the basic solution. 2
As shown above, the time-averaged MSE is minimized through (11). However, since this error is WSCS(D), the expected squared error will generally vary periodically with the time index. This prompts the question: Is this expected squared error minimized at all times? By considering the errors of all D -fold decimated output sequences it will be shown that yes, the solution (11) does not periodically sacrifice instantaneous performance to optimize average performance. The left portion of Fig. 2 shows the block diagram producing one of these decimated sequences, defined through
ek [n] = e[nD + k] = x[nD + k] 0 y[nD + k]
2 e;k = lim N !1
1
(9)
0 j! (e ) is Thus, the TAMSE is a function of the trace of (8). Since See strictly nonnegative, the minimization of e2 is equivalent to the minimization of tr(See (ej! )). Since (8) is only defined for j! j =D ,
N 01
R [nD + k; nD + k]: N n=0 ee
(16)
The system on the left of Fig. 2 can be simplified to the equivalent system on the right, where Ei;k (ej! ) is the k th polyphase component of Fi (ej! ) based on the type-1 decomposition [2]
Fi (ej! ) =
0 j! k See (e W )d!
(15)
for k = 0; 1; . . . ; D 0 1. Each of these subsequences is WSS with an MSE determined by
III. SYNTHESIS FILTER OPTIMIZATION This section determines the frequency domain representation for the synthesis bank for minimizing the system’s MSE where the input signal statistics and sampling process (modelled through the analysis bank) are assumed known. The approach is especially of interest for the case of D > C , which generally results in a loss of information, thereby negating the possibility of a PR solution. The TAMSE is determined through an equivalent to (4) as 1
j!j D :
A. Optimality of Decimated Output Sequences (7)
CSD matrices relating the input x[n] and output y [n] are defined as Syy (ej! ) = P(ej! )Sxx (ej! )PH (ej! ), Syx (ej! ) = P(ej! )Sxx (ej! ), and Sxy (ej! ) = SHyx (ej! ) = SHxx (ej! )PH (ej! ). The superscript H denotes the conjugate transpose. The CSD matrix corresponding to e[n] = x[n] 0 y [n], the reconstruction error of the filter bank, is found to be
e2 =
for all
where
=
condition for a WSS input. The analysis and synthesis filters are represented using alias component (AC) matrices [2], defined as
e
(
The solution, detailed in the Appendix, is found to be
To determine a filter bank’s CSD input/output relationship, a matrix representation for the CSDs of (3) is defined. This D 2 D CSD matrix is expressed as Sxx (ej! ), and the (p; q )th element is defined as
Sxx (ej! )
min
D01 k=0
e0j!k Ei;k (ej!D ):
(17)
The mean-squared variance (15) can therefore be minimized through selection of an optimal set of Ei;k (ej! ). As will be proven below, combining the optimal Ei;k (ej! ) through (17) determines the same optimal synthesis bank of (11), showing FAC;opt (ej! ) minimizes the MSE for all time indexes.
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
3240
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
Fig. 3. Dual-analysis bank model used to consider the presence of an intermerging analfering source. The system can be generalized to consist of in the presence of distinct interfering sources. ysis banks to model
x[n]
Fig. 2. Equivalent reconstruction models used to determine the MSE of a decimated output sequence.
The k th set of polyphase components is represented through the vector E k (ej! ) = [E0;k (ej! ); . . . ; EC 01;k (ej! )]. Similarly to (11), the mean-square variance of ek [n] is minimized through selection of the optimal filters
E k;opt (ej!D ) = hTk (ej! )QH (ej! )R01 (ej! )
(18)
where hTk (ej! ) = ej!k [1; W k ; . . . ; W k(D01) ]. As before, (18) is only valid for j! j =D . The alias components of the polyphase representation (17) are found to be
Fi (ej! W p ) = =
D01 k=0 D01 k=0
e0j!k W 0pk Ei;k (ej!D W pD ) e0j!k W 0pk Ei;k (ej!D ):
These alias components make up the entries of plying (18) to (19), is written as
F
AC (e
j! ) =
F
AC (
(19)
ej! ) which, ap-
A(ej! )B(ej! )QH (ej! )R0 (ej! ) 1
where
1
A(ej! ) =
1
.. .
1 and
B(ej! ) = A
e0j! 0 e j! W 01 .. .
e0j! W 0(D01)
111 111
e0j!(D01) 0 j! (D 01) e W 0(D01)
.
e0j!(D01) W 0(D01)
..
111
hT(D01) (ej! )
B
:
2 e;k =
D =D Se e (ej! )d! = Se e (ej!D )d! (22) 2 0 2 0=D 1
where 1
Se e (ej!D ) = hTk (ej! ) D 1 Sxx(ej! ) 0 QH (ej! )R0H (ej! )Q(ej! ) h3k (ej! ): (23) Solving (22) for 0 k D 0 1 allows the expected squared error at any given time index to be found. An important implication of this result is that, since all D-fold decimated sequences are individually optimal, the output does not need to be implemented at the model rate TM . In fact, since selection of D is only bounded by the minimum value necessary for model accuracy, an output signal can be produced at nearly any rate (it must be an integer ratio of 1=TP ). This enables lower cost reconstruction implementations without reducing the modeling accuracy of x[n]. B. Additive Noise and Interference The dual analysis bank model of Fig. 3 is used to consider w[n], an additive noise or interference signal with known statistics. This input has the same requirements and assumptions as x[n]. The separate set of analysis filters G0 (ej! ); . . . ; GC 01 (ej! ) allows w[n] to be considered through either the same sampling process as x[n] (in which case j! ) = j! )) or a unique sampling process (e.g., to model AC (e AC (e multipath distortions of an interfering signal). One useful configuration is an advance chain, where Gl (ej! ) = ej!l . This considers the subband cross correlations using the WSCS(D) autocorrelation of w[n] through
H
E [w[nD + p]w3 [mD + q]] = E vp [n]vq3 [m]
.. .
hT0 (ej! ) .. .
tions. Additionally, through (16) the decimated output MSE is given by
G
(20)
K +1
K
(21)
Realizing (ej! ) = D 01 (ej! ), the result (20) simplifies to the original optimal solution (11), thereby showing equivalence of the solu-
(24)
for all n, m, and p; q = 0; 1; . . . ; C 0 1. For C p; q D 0 1, (24) can be set to zero. This configuration models noise that is introduced subsequent to the sampling process. To further generalize, K distinct interfering sources (w1 [n]; . . . ; wK [n]) can be considered using K + 1 merging banks of analysis filters. This framework is useful for considering generalized multiple-input multiple-output (MIMO) problems, since a separate synthesis bank can be determined for each output. For this general problem, a separate alias component matrix is constructed for each analysis bank. The matrix corresponding to the k th interfering source wk [n] is defined as through (6) and labeled k;AC (ej! ).
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
G
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
Using this generalization, the optimal reconstruction filter bank and error CSD matrix are still determined, respectively, through (11) and (14). However, the component matrices Q(ej! ) and R(ej! ) must be modified to account for the multiple source structure. The new definitions for multiple sources are
Q(ej! ) = HTAC (ej! )Sxx (ej! ) +
K k=1
GTk;AC (ej! )Sxw
(e
T (ej! )SH (ej! )H3 (ej! ) R(ej! ) = HAC xx AC K
+
+
+
l=1 K
T (ej! )SH HAC xw
(e
j! )
k=1 l=1
IV. CONCLUDING REMARKS AND PROBLEM EXTENSIONS
j! )G3 (ej! ) l;AC
GTk;AC (ej! )SwH w (ej! )Gl;3 AC (ej! ):
An alternate numerical solution was presented for this problem in the earlier work of Ohno and Sakai [5]. Their approach examined the error caused by missing subband of a biorthogonal filter bank in the correlation domain. Numerical optimization techniques were then used to find synthesis filters of predetermined length. Although their solution could not guarantee the globally optimal solution, it did have the advantage of allowing specific design constraints (e.g., linear phase synthesis filters).
(25)
3 (ej! ) GTk;AC (ej! )SwH x (ej! )HAC
k=1 K K
3241
(26)
It is worth noting that for a given set of analysis banks modeling, a MIMO scenario R(ej! ) does not change with respect to selection of the desired signal. This is useful for reducing the design complexity since each distinct synthesis bank is determined by applying a unique Q(ej! ) to a common R01 (ej! ).
A filter bank framework for generalized sampling without strict bandlimiting conditions was demonstrated. Optimization techniques were used to determine a set of filters to minimize the TAMSE of reconstruction for a given signal and sampling model. A generalized form of the problem was also presented for cases of additive noise and interfering sources. An example sampling scenario using this presented approach can be found in [27]. There also exist several extensions using the foundation presented in this correspondence. A similar method was used in [28] for signal reconstruction in the case of unknown random sampler timing offsets. A multidimensional approach is used in [29] for image super-resolution. Current research investigating the optimization of periodic sampling patterns is ongoing. APPENDIX CALCULATION OF (11)
C. Numerical Considerations While an analytical solution has been considered so far, a feasible solution requires a numerical approach. For this, the solution (11) is determined for varying frequencies over j! j =D , providing sampled versions of the synthesis filters’ frequency responses through piecewise combinations of the alias components. The sampled frequency responses should be of sufficiently fine resolution from which accurate filter realizations can be determined. A high-resolution model will 2 generally also allow an accurate calculation of e; opt through numerical integration of (9). For the solution (11) to exist at a particular ! , the matrix R(ej! ) of (13) must be invertible. As mentioned in Section II-B, Sxx (ej! ) is diagonal for a WSS x[n], with real nonnegative terms representing j! j! portions of the power spectrum. This means SH xx (e ) = Sxx (e ), allowing the (p; q)th entry of R(ej! ) to be found through
R(ej! )
D01 (
p;q)
=
k=0
Sxx (e
j! W k )H
p (e
j! W k )H 3 (ej! W k ) (27) q
where (p; q) = 0; . . . ; C 0 1. To be invertible, R(ej! ) must be of full rank C , requiring diagonal Sxx (ej! ) to have at least C nonzero entries. This is not true at frequencies for which there are fewer nonzero aliased components than samplers. Although R(ej! ) is not invertible in these cases, the reduced number of aliased components actually indicates oversampled frequencies which can be made free of reconstruction error. To circumvent this problem, the zero-valued portions of Sxx (ej! ) can be replaced with a small positive value (much smaller than the non-zero valued spectral components), thereby providing an invertible approximation of R(ej! ). The unaltered Sxx (ej! ) is used for Q(ej! )(12), the noninverted portion of the solution, ensuring zero-valued portions of the input spectrum remain so at the output. Invertibility of (27) also places some restrictions on HAC (ej! ), generally requiring unique timing offsets or distortions for the samplers. In the event of duplicate sampling, the redundant branches of the model should be removed, corresponding to the removal of linearly dependent columns from HAC (ej! ).
The solution to (10) is found by setting the derivative of See (ej! )) with respect to FAC (ej! ) equal to zero (the relevant operations can be found in [30, App. E]). This produces tr(
3 j! H j! T j! j! FAC ;opt (e )(HAC (e )Sxx (e )HAC (e )) =D
STxx (ej! )HAC (ej! )
from which (11) is obtained provided the inverse of H (ej! )ST (ej! )H (ej! ) exists. HAC AC xx To ensure the derivative’s zero-point is a global minimum, the convexity of tr(See (ej! )) is examined. Ordering the entries of FAC (ej! ) vectorially through the concatenation of its rows, the Levi matrix (or, complex Hessian) [31] of tr(See (ej! )) with respect to FAC (ej! ) is calculated. The result is a block diagonal matrix with identical size j! 2 C 2 C blocks equal to R(e )=D , a scaled version of (13). This structure means convexity is proven if and only if R(ej! ) is positive definite. Although a detailed examination of the necessary conditions for a positive definite R(ej! ) is not provided here, Section III-C discusses some conditions causing a noninvertible R(ej! ), resulting in a semidefinite matrix for which a unique minimum cannot be found.
REFERENCES [1] A. Papoulis, “Generalized sampling expansion,” IEEE Trans. Circuits Syst., vol. CAS-24, no. 11, pp. 652–654, Nov. 1977. [2] P. P. Vaidyanathan, Multirate Systems and Filter Banks. EnglewoodCliffs, NJ: Prentice-Hall, 1993. [3] M. K. Tsatsanis and G. B. Giannakis, “Principal component filter banks for optimal multiresolution analysis,” IEEE Trans. Signal Process., vol. 43, no. 8, pp. 1766–1777, Aug. 1995. [4] R. A. Haddad and K. Park, “Modeling, analysis, and optimum design of quantized M-band filter banks,” IEEE Trans. Signal Process., vol. 43, no. 11, pp. 2540–2549, Nov. 1995. [5] S. Ohno and H. Sakai, “Optimization of filter banks using cyclostationary spectral analysis,” IEEE Trans. Signal Process., vol. 44, no. 11, pp. 2718–2725, Nov. 1996.
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
3242
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
[6] K. Gosse and P. Duhamel, “Perfect reconstruction versus MMSE filter banks in source coding,” IEEE Trans. Signal Process., vol. 45, no. 9, pp. 2188–2202, Sep. 1997. [7] P. P. Vaidyanathan, “Theory of optimal orthonormal subband coders,” IEEE Trans. Signal Process., vol. 46, no. 6, pp. 1528–1543, Jun. 1998. [8] P. P. Vaidyanathan and A. Kirac, “Results on optimal biorthogonal filter banks,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 45, no. 8, pp. 932–947, Aug. 1998. [9] J. Tuqan and P. P. Vaidyanathan, “Oversampling PCM techniques and optimum noise shapers for quantizing a class of nonbandlimited signals,” IEEE Trans. Signal Process., vol. 47, no. 2, pp. 389–407, Feb. 1999. [10] P. P. Vaidyanathan and V. C. Liu, “Efficient reconstruction of band-limited sequences from nonuniformly decimated versions by use of polyphase filter banks,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 38, no. 11, pp. 1927–1936, Nov. 1990. [11] A. J. Coulson, “A generalization of nonuniform bandpass sampling,” IEEE Trans. Signal Process., vol. 43, no. 3, pp. 694–704, Mar. 1995. [12] P. P. Vaidyanathan and S.-M. Phoong, “Reconstruction of sequences from nonuniform samples,” in Proc. IEEE Int. Symp. Circuits Systems (ISCAS), Seattle, WA, 1995, vol. 1, pp. 601–604. [13] C. Herley and P. W. Wong, “Minimum rate sampling and reconstruction of signals with arbitrary frequency support,” IEEE Trans. Inf. Theory, vol. 45, no. 5, pp. 1555–1564, Jul. 1999. [14] Y. C. Eldar and A. V. Oppenheim, “Filterbank reconstruction of bandlimited signals from nonuniform and generalized samples,” IEEE Trans. Signal Process., vol. 48, no. 10, pp. 2864–2875, Oct. 2000. [15] P. P. Vaidyanathan, “Generalizations of the sampling theorem: Seven decades after Nyquist,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 48, pp. 1094–1109, Sep. 2001. [16] R. Venkataramani and Y. Bresler, “Optimal sub-Nyquist nonuniform sampling and reconstruction for multiband signals,” IEEE Trans. Signal Process., vol. 49, no. 10, pp. 2301–2313, Oct. 2001. [17] H. Johansson and P. Löwenborg, “Reconstruction of nonuniformly sampled bandlimited signals by means of digital fractional delay filters,” IEEE Trans. Signal Process., vol. 50, no. 11, pp. 2757–2767, Nov. 2002. [18] R. S. Prendergast, B. C. Levy, and P. J. Hurst, “Reconstruction of bandlimited periodic nonuniformly sampled signals through multirate filter banks,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 51, pp. 1612–1622, Aug. 2004. [19] A. S. Bopardikar, R. M. Rao, and B. S. Adiga, “Matched sampling systems, relation to wavelets and implementation using PRCC filter banks,” IEEE Trans. Signal Process., vol. 48, no. 8, pp. 2269–2278, Aug. 2000. [20] D. Seidner and M. Feder, “Vector sampling expansion,” IEEE Trans. Signal Process., vol. 48, no. 5, pp. 1401–1416, May 2000. [21] R. Venkataramani and Y. Bresler, “Filter design for MIMO sampling and reconstruction,” IEEE Trans. Signal Process., vol. 51, no. 12, pp. 3164–3176, Dec. 2003. [22] M. Unser and J. Zerubia, “Generalized sampling: Stability and performance analysis,” IEEE Trans. Signal Process., vol. 45, no. 12, pp. 2941–2950, Dec. 1997. [23] ——, “A generalized sampling theory without band-limiting constraints,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 45, no. 8, pp. 959–969, Aug. 1998. [24] M. Unser, “Sampling-50 years after Shannon,” Proc. IEEE, vol. 88, pp. 569–587, Apr. 2000. [25] T. I. Laakso, V. Välimäki, M. Karjalainen, and U. K. Laine, “Splitting the unit delay,” IEEE Signal Process. Mag., vol. 13, no. 1, pp. 30–60, Jan. 1996. [26] W. A. Gardner, Introduction to Random Processes with Applications to Signals and Systems. New York: McGraw-Hill, 1989. [27] R. S. Prendergast and T. Q. Nguyen, “Optimal filter bank reconstruction of periodically undersampled signals,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing (ICASSP), Philadelphia, PA, 2005, vol. 4, pp. 201–204.
[28] ——, “Optimal reconstruction of periodically sampled signals with probabilistic timing delays,” presented at the 38th Asilomar Conf. Signals, Systems, Computers Pacific Grove, CA, Nov. 2004. [29] ——, “Improving frequency domain super-resolution via undersampling model,” in IEEE Int. Conf. Image Processing, 2005, pp. 853–856. [30] T. K. Moon and W. C. Stirling, Mathematical Methods and Algorithms for Signal Processing. Upper Saddle River, NJ: Prentice-Hall, 2000. [31] S. G. Krantz, Function Theory of Several Complex Variables. Providence, RI: AMS Chelsea, 1992.
Optimum Blind Multichannel Equalization Using the Linear Prediction Algorithm Houcem Gazzah, Member, IEEE
Abstract—The linear prediction (LP) algorithm estimates a zero-forcing (ZF) equalizer from the channel output’s second-order statistics (SOS). We show that it allows, in fact, to determine a whole subspace of (but not all) ZF equalizers. The one obtained by maximizing the signal-to-noise ratio (SNR) at the equalizer output not only outperforms the original LP equalizer, but also attains the lowest achievable (by any no-delay ZF equalizer) mean-square error (MSE). The increase in complexity is only marginal, but equalization performance can be enhanced significantly, especially in the (practical) case of a small leading channel coefficient. Index Terms—Channel modeling, equalization, estimation.
I. INTRODUCTION Training sequences are becoming inadequate with the growing demand on bandwidth. This has motivated an intensive interest in blind techniques, and more particularly, those based on (the efficiently estimated) second-order statistics (SOS). This is made possible if the channel output is oversampled and/or received using multiple sensors. The linear prediction (LP) algorithm [1], [2] is among the most popular SOS-based approaches. The LP equalizer is obtained by linearly combining the columns of the predictor. Each, in fact, can be regarded as a zero-forcing (ZF) equalizer, and each has a different behavior in the presence of observation noise. This fact was not considered in the original LP algorithm, and hence, the computed equalizer is not optimal. In this correspondence, we optimize the LP algorithm using the signal-to-noise ratio (SNR) criterion [3], [4] (Section IV). We show the so-optimized LP equalizer to be nothing but the ZF minimum meansquare error (MMSE) zero-delay equalizer. The gain of the optimized solution is especially relevant in the practical case of a small leading channel coefficient. A performance analysis (Section V) and simulations (Section VI) support these claims. We denote by Diag(a1 ; . . .) the diagonal matrix whose ith diagonal entry is ai , by Ia the a 2 a identity matrix and by 0a;b the a 2 b zero matrix. Indexes are omitted whenever the dimensions can be inferred from the context. The ith row, j th column entry of the matrix Manuscript received February 15, 2005; revised October 9, 2005. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Helmut Boelcskei. The author is with the Institute for Digital Communications, University of Edinburgh, Edinburgh, EH9 3JL U.K. (e-mail: [email protected]). Digital Object Identifier 10.1109/TSP.2006.877669
1053-587X/$20.00 © 2006 IEEE
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
3237
On the other hand, we have shown that B (^b) is idempotent, that is, B (^b) = B 2 (^b), which means that it only has eigenvalues 0 and 1 and > 0. Combining the results, the zero-phase property hence [B (^b)]
Minimum Mean-Squared Error Reconstruction for Generalized Undersampling of Cyclostationary Processes
REFERENCES
Ryan S. Prendergast, Student Member, IEEE, and Truong Q. Nguyen, Fellow, IEEE
m;m
of the denominator of (13) is proved.
[1] S. M. Kay, Modern Spectral Estimation : Theory and Application. Englewood Cliffs, NJ: Prentice-Hall, 1988. [2] S. Kay and R. Nekovei, “An efficient two-dimensional frequency estimator,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 38, no. 10, pp. 1807–1809, Oct. 1990. [3] Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” IEEE Trans. Signal Process., vol. 40, no. 9, pp. 2267–2280, Sep. 1992. [4] M. P. Clark and L. L. Scharf, “Two-dimensional modal analysis based on maximum likelihood,” IEEE Trans. Signal Process., vol. 42, no. 6, pp. 1443–1452, Jun. 1994. [5] C. R. Rao, L. Zhao, and B. Zhou, “Maximum likelihood estimation of 2-D superimposed exponential signals,” IEEE Trans. Signal Process., vol. 42, no. 7, pp. 1795–1802, Jul. 1994. [6] A. Sourice, G. Plantier, and J.-L. Saumet, “Two-dimensional frequency estimation using autocorrelation phase fitting,” in Proc. Int. Conf. Acoustics, Speech, Signal Processing (ICASSP), Hong Kong, Apr. 2003, vol. 3, pp. 445–448. [7] P. Stoica and R. Moses, Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice-Hall, 2005. [8] S. A. Tretter, “Estimating the frequency of a noisy sinusoid by linear regression,” IEEE Trans. Inf. Theory, vol. 31, pp. 832–835, Nov. 1985. [9] G. W. Lank, I. S. Reed, and G. E. Pollon, “A semicoherent detection and Doppler estimation statistic,” IEEE Trans. Aerosp. Electron. Syst., vol. 9, pp. 151–165, Mar. 1973. [10] S. Kay, “A fast and accurate single frequency estimator,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 12, pp. 1987–1990, Dec. 1989. [11] H. Chien, “2-D estimation from AR models,” Ph.D. dissertation, Elect. Eng. Dept., Univ. of Rhode Island, Kingston, RI, 1981. [12] R. Kumaresan, L. L. Scharf, and A. K. Shaw, “An algorithm for polezero modeling and spectral analysis,” IEEE Trans. Acoust., Speech, Signal Process., vol. 34, no. 6, pp. 637–640, Jun. 1986. [13] Y. Bresler and A. Macovski, “Exact maximum likelihood parameter estimation of superimposed exponential signals in noise,” IEEE Trans. Acoust., Speech, Signal Process., vol. 34, no. 10, pp. 1081–1089, Oct. 1986. [14] G. C. Goodwin and R. L. Payne, Dynamic System Identification: Experiment Design and Data Analysis. New York: Academic, 1977, pp. 46–50. [15] H. C. So and K. W. Chan, “A generalized weighted linear predictor frequency estimation approach for a complex sinusoid,” IEEE Trans. Signal Processing, vol. 54, no. 4, pp. 1304–1315, Apr. 2006. [16] B. Völcker and P. Händel, “Frequency estimation from proper sets of correlations,” IEEE Trans. Signal Process., vol. 50, no. 4, pp. 791–802, Apr. 2002. [17] M. L. Fowler and J. A. Johnson, “Extending the threshold and frequency range for phase-based frequency estimation,” IEEE Trans. Signal Process., vol. 47, no. 10, pp. 2857–2863, Oct. 1999. [18] H. Lutkepohl, Handbook of Matrices. New York: Wiley, 1996.
Abstract—The generalized sampling problem is considered using a filter bank model with a stochastic process framework. A signal reconstruction solution minimizing the time-averaged mean-squared error is found using a filter bank optimization technique for cyclostationary signals. The use of a stochastic model approach to find a minimized error solution rather than a perfect reconstruction solution allows consideration of sampling problems that deterministic approaches cannot solve. For instance, a particular sampling density is not required, allowing reconstruction optimization in cases of undersampling. An extension is also presented for optimal reconstruction in the presence of additive noise and interference. The result is a generalized sampling model which can be used for an extensive variety of scenarios, combining aspects of undersampling, multiple-input multiple-output (MIMO) sampling, and periodic nonuniform sampling. Index Terms—Cyclostationary processes, generalized sampling, minimum mean-square error (MMSE) filtering, undersampling.
I. INTRODUCTION This correspondence examines a modified approach to the generalized sampling problem using a filter bank model with a discrete-time cyclostationary input. Traditional implementations of Papoulis’ generalized sampling expansion [1] have considered sampling of several linearly filtered versions of a bandlimited continuous-time signal at a collective Nyquist rate, from which perfect reconstruction (PR) can be achieved. Since the sampled signal is perfectly bandlimited, an equivalent discrete-time model is found using the multirate filter bank, and solutions can be determined based on established PR filter bank theory [2]. While previous generalized sampling approaches have mainly dealt with deterministic signals of known finite spectral support, this correspondence will consider stochastic signals with known spectral densities. The solution does not require a minimum sampling density to determine a signal reconstruction. Instead, the optimal reconstruction is found for a given sampling process. This provides a distinct advantage over the traditional approach, allowing for optimal reconstruction solutions in the event of nonuniform undersampling and in the presence of interference. As such, the result is directly applicable to undersampling problems such as image super-resolution, and may also be of interest for sampling rate reduction of analog-to-digital conversion and equalization systems. Much previous literature in the area of filter banks have shown that the spectral characteristics of stochastic signals can be taken advantage of to determine optimal filtering structures for a variety of applications [3]–[9]. The generalized sampling problem requires a different approach than many other filter bank optimization problems since it must incorporate a fixed model based on the sampling process. In other cases, typically both the analysis and synthesis filter banks (see Fig. 1) are jointly optimized. The filter bank model for a sampling problem Manuscript received March 4, 2005; revised August 28, 2005. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. H. Boelcskei. This research has been supported by a grant from Skyworks, Inc. and UC Dimi. The authors are with the Department of Electrical and Computer Engineering, University of California, San Diego, La Jolla, CA 92093 USA (e-mail: [email protected]; [email protected]). Digital Object Identifier 10.1109/TSP.2006.877649
1053-587X/$20.00 © 2006 IEEE
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
3238
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
x(nTM ) with a sufficiently short model period TM
Fig. 1. General multirate filter bank structure.
uses the analysis bank to represent a physical sampling process and any enhancement or processing is confined to the synthesis bank. As such, the general problem to be solved can be stated as follows. Given a particular generalized sampling process (represented through analysis filters H0 (ej! ); . . . ; HC 01 (ej! )), determine a linear processing technique (through synthesis filters F0 (ej! ); . . . ; FC 01 (ej! )) that is optimal for some criterion. The presented approach determines an optimal solution through minimization of the mean of je[n]j2 = jx[n] 0 y [n]j2 , the squared error. Since a stationary signal model and linear synthesis filters are used, the result is a Wiener filtering solution for the generalized sampling problem. The cyclostationary filter bank analysis techniques developed by Ohno and Sakai [5] are used to determine a frequency domain solution for the synthesis filters. For the purposes of this correspondence, the sampled signal will be considered using a discrete-time model. The limitations associated with the use of a discrete-time model for real-world signals can be dismissed since the resolution of this model can be arbitrarily increased. A considerable amount of work has been done using filter bank models of generalized sampling since Papoulis’ result. One problem of interest is a purely periodic sampling problem, which is considered through Fig. 1 when the analysis filters are simple delays. This has been considered for signals of finite spectral support in a number of different cases [10]–[18]. Primary interest is in periodic nonuniform sampling, since results for uniform sampling are well established. The use of periodic sampling patterns of sufficient average rate (equal to the Nyquist rate for lowpass signals) allows PR solutions to be found, although practical implementations can only approximate PR. Some modifications to the standard problem include undersampling [19], multi-input problems [20], [21], and reconstruction using alternative (to finite spectral support) bases [22]–[24]. This correspondence’s method is capable of dealing with a wide variety of sampling problems using a minimized mean-square error (MSE) approach. In some instances, PR solutions can be found, but optimal reconstruction solutions are still found when PR cannot be achieved. The organization of this correspondence is as follows. Section II provides background on the use of a discrete-time framework for generalized sampling problems and on filter bank analysis of cyclostationary signals. Section III presents the general optimal solution for the sampling problem with some related discussion and an extension allowing consideration of multiple inputs and additive noise. Concluding remarks are provided in Section IV, along with a brief mention of problem extensions and future work. II. PRELIMINARIES
enables representation without loss of information. In order to maintain the periodicity of the sampling process in a filter bank representation, it will be required that the sampling process period TP be divisible by TM , or DTM = TP , where D is a positive integer corresponding to the degree of decimation in the filter bank model. For a given TP , the decimation factor D is herein assumed to be sufficiently large so that x[n] = x(nTP =D) provides a lossless representation of x(t). Through this assumption, strictly discrete-time signal models can be considered. It should be noted that use of the discrete-time representation will generally not limit the modelling of a physical periodic sampling process. Referring to Fig. 1, each branch of the analysis filter bank models a physical uniform sampling with a unique delay and distortion function. The collection of these C samplers forms the periodic sampling process. In absence of distortion, each sampler is represented through the delay term Hi (ej! ) = e0j!d : Allowing the use of fractional delays enables consideration of any periodic sampling process. It is known that fractional delays can only be perfectly represented through filters of infinite length in the time domain [25]. However, as is done in this correspondence, ideal fractional delay filters are easily considered in the frequency domain over the limited range j! j < . B. Filter Bank Analysis of Cyclostationary Signals This subsection introduces definitions based on the analysis developed by Ohno and Sakai in [5] for the purpose of minimizing a filter bank’s time-averaged MSE (TAMSE) when a subband is lost. The approach considers a discrete-time cyclostationary input even though primary interest will be for wide-sense stationary (WSS) inputs. The reason for this is that even when the filter bank is provided a stationary input, it will generally produce an output that is cyclostationary with period D (corresponding to the order of decimation). Because of this, the error between the input and output signals is cyclostationary. For generalization purposes, the input is considered so as well (since stationarity is a subset of cyclostationarity). Specifically, x[n] is wide-sense cyclostationary with period D , or WSCS(D), meaning its first- and second-order expectations are periodically stationary with period D . This is expressed as
E [x[n]] = E [x[n + kD]]
and Rxx [n; m] = Rxx [n + kD; m + kD ]
(1)
for every integer k , where Rxx [n; m] = E [x[n]x 3 [m]]. Herein, a zeromean process is considered. The cyclic correlation function is defined as
[ u] = Rxx
D01
R [k + u; k]e0j 2 k D k=0 xx 1
(2)
for = n=D and integer n (the notation used in [5] is based on that of [26]). Since the exponential term is periodic with , (2) only needs to be determined for 0 (D 0 1)=D (or, n = 0; 1; . . . ; D 0 1) to be defined for all valid . The discrete-time Fourier transforms of the cyclic correlation functions produce the cyclic spectral densities (CSDs) through
A. Discrete-Time Framework for Sampling Model This section describes the discrete-time framework used to model a periodically sampled continuous-time signal. Provided this continuous-time signal x(t) is bandlimited, a discrete-time model x[n] =
(ej! ) = Sxx
1 u=01
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
[u]e0j!u : Rxx
(3)
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
Through (2) and (3), the time-averaged variance of a WSCS(D) process is defined as 1
x2 = Nlim !1
N 01
N n=0
= Rxx [0] = 0
E [x[n]x3 [n]] =
1
2
0
1
D01
R [n; n] D n=0 xx 2 0 j! Sxx (e )d!:
3239
the minimization of (9) with respect to the synthesis AC matrix can be stated as
F (4)
(
p;q)
=
(p0q )=D j! p Sxx (e W )
over
j!j D
(5)
p; q = 0; . . . ; D 0 1, using the standard definition of e0j (2=D) . This representation is only valid for j!j =D ; however, this is sufficient to completely represent all D CSDs of (3). [u] = 0 for all valid 6= 0 if and only if x[n] is WSS, it is Since Rxx easily verified that a diagonal Sxx (ej! ) is a necessary and sufficient for
W
111
H0 (ej! )
HAC (ej! ) =
.. .
..
H0 (ej! W (D01) )
HC 01 (ej! ) .. .
.
(6)
1 1 1 HC01 (ej! W (D01) ) for the analysis bank and an equivalent representation FAC (ej! ) for the synthesis bank. Like the CSD matrix, these AC representations are also only considered over j! j =D but completely and uniquely describe the analysis and synthesis filters in the frequency domain. Defining the matrix product
P(ej! ) = D1 FAC (ej! )HTAC (ej! )
)
tr(
See (ej! ))
See (ej! ) = Sxx (ej! ) 0 Syx (ej! ) 0 Sxy (ej! ) + Syy (ej! ) j! j! j! H j! H j! = Sxx (e ) 0 P(e )Sxx (e ) 0 Sxx (e )P (e ) j! j! H j! (8) + P(e )Sxx (e )P (e ) which represents the error as a function of the filter bank and the input signal statistics.
=
2
2 1
2
0
0 j! See (e )d! =
=D
0=D
tr(
1
2
See (ej! ))d!:
D01
=D
k=0
0=D
(10)
FAC;opt (ej! ) = DQH (ej! )R01 (ej! )
(11)
T (ej! )S (ej! ) Q(ej! ) = HAC xx
(12)
T (ej! )SH (ej! )H3 (ej! ): R(ej! ) = HAC xx AC
(13)
and
An expression for the minimized TAMSE is found by applying the solution (11) to (8), resulting in the optimal error CSD matrix
See;opt (ej! ) = Sxx (ej! ) 0 QH (ej! )R0H (ej! )Q(ej! ) H j! QH (ej! )FAC ;opt (e ) : j! = Sxx (e ) 0 D
(14)
As through (9), the optimal TAMSE e;opt can be determined from this expression. The remainder of this section provides further analysis and extensions on the basic solution. 2
As shown above, the time-averaged MSE is minimized through (11). However, since this error is WSCS(D), the expected squared error will generally vary periodically with the time index. This prompts the question: Is this expected squared error minimized at all times? By considering the errors of all D -fold decimated output sequences it will be shown that yes, the solution (11) does not periodically sacrifice instantaneous performance to optimize average performance. The left portion of Fig. 2 shows the block diagram producing one of these decimated sequences, defined through
ek [n] = e[nD + k] = x[nD + k] 0 y[nD + k]
2 e;k = lim N !1
1
(9)
0 j! (e ) is Thus, the TAMSE is a function of the trace of (8). Since See strictly nonnegative, the minimization of e2 is equivalent to the minimization of tr(See (ej! )). Since (8) is only defined for j! j =D ,
N 01
R [nD + k; nD + k]: N n=0 ee
(16)
The system on the left of Fig. 2 can be simplified to the equivalent system on the right, where Ei;k (ej! ) is the k th polyphase component of Fi (ej! ) based on the type-1 decomposition [2]
Fi (ej! ) =
0 j! k See (e W )d!
(15)
for k = 0; 1; . . . ; D 0 1. Each of these subsequences is WSS with an MSE determined by
III. SYNTHESIS FILTER OPTIMIZATION This section determines the frequency domain representation for the synthesis bank for minimizing the system’s MSE where the input signal statistics and sampling process (modelled through the analysis bank) are assumed known. The approach is especially of interest for the case of D > C , which generally results in a loss of information, thereby negating the possibility of a PR solution. The TAMSE is determined through an equivalent to (4) as 1
j!j D :
A. Optimality of Decimated Output Sequences (7)
CSD matrices relating the input x[n] and output y [n] are defined as Syy (ej! ) = P(ej! )Sxx (ej! )PH (ej! ), Syx (ej! ) = P(ej! )Sxx (ej! ), and Sxy (ej! ) = SHyx (ej! ) = SHxx (ej! )PH (ej! ). The superscript H denotes the conjugate transpose. The CSD matrix corresponding to e[n] = x[n] 0 y [n], the reconstruction error of the filter bank, is found to be
e2 =
for all
where
=
condition for a WSS input. The analysis and synthesis filters are represented using alias component (AC) matrices [2], defined as
e
(
The solution, detailed in the Appendix, is found to be
To determine a filter bank’s CSD input/output relationship, a matrix representation for the CSDs of (3) is defined. This D 2 D CSD matrix is expressed as Sxx (ej! ), and the (p; q )th element is defined as
Sxx (ej! )
min
D01 k=0
e0j!k Ei;k (ej!D ):
(17)
The mean-squared variance (15) can therefore be minimized through selection of an optimal set of Ei;k (ej! ). As will be proven below, combining the optimal Ei;k (ej! ) through (17) determines the same optimal synthesis bank of (11), showing FAC;opt (ej! ) minimizes the MSE for all time indexes.
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
3240
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
Fig. 3. Dual-analysis bank model used to consider the presence of an intermerging analfering source. The system can be generalized to consist of in the presence of distinct interfering sources. ysis banks to model
x[n]
Fig. 2. Equivalent reconstruction models used to determine the MSE of a decimated output sequence.
The k th set of polyphase components is represented through the vector E k (ej! ) = [E0;k (ej! ); . . . ; EC 01;k (ej! )]. Similarly to (11), the mean-square variance of ek [n] is minimized through selection of the optimal filters
E k;opt (ej!D ) = hTk (ej! )QH (ej! )R01 (ej! )
(18)
where hTk (ej! ) = ej!k [1; W k ; . . . ; W k(D01) ]. As before, (18) is only valid for j! j =D . The alias components of the polyphase representation (17) are found to be
Fi (ej! W p ) = =
D01 k=0 D01 k=0
e0j!k W 0pk Ei;k (ej!D W pD ) e0j!k W 0pk Ei;k (ej!D ):
These alias components make up the entries of plying (18) to (19), is written as
F
AC (e
j! ) =
F
AC (
(19)
ej! ) which, ap-
A(ej! )B(ej! )QH (ej! )R0 (ej! ) 1
where
1
A(ej! ) =
1
.. .
1 and
B(ej! ) = A
e0j! 0 e j! W 01 .. .
e0j! W 0(D01)
111 111
e0j!(D01) 0 j! (D 01) e W 0(D01)
.
e0j!(D01) W 0(D01)
..
111
hT(D01) (ej! )
B
:
2 e;k =
D =D Se e (ej! )d! = Se e (ej!D )d! (22) 2 0 2 0=D 1
where 1
Se e (ej!D ) = hTk (ej! ) D 1 Sxx(ej! ) 0 QH (ej! )R0H (ej! )Q(ej! ) h3k (ej! ): (23) Solving (22) for 0 k D 0 1 allows the expected squared error at any given time index to be found. An important implication of this result is that, since all D-fold decimated sequences are individually optimal, the output does not need to be implemented at the model rate TM . In fact, since selection of D is only bounded by the minimum value necessary for model accuracy, an output signal can be produced at nearly any rate (it must be an integer ratio of 1=TP ). This enables lower cost reconstruction implementations without reducing the modeling accuracy of x[n]. B. Additive Noise and Interference The dual analysis bank model of Fig. 3 is used to consider w[n], an additive noise or interference signal with known statistics. This input has the same requirements and assumptions as x[n]. The separate set of analysis filters G0 (ej! ); . . . ; GC 01 (ej! ) allows w[n] to be considered through either the same sampling process as x[n] (in which case j! ) = j! )) or a unique sampling process (e.g., to model AC (e AC (e multipath distortions of an interfering signal). One useful configuration is an advance chain, where Gl (ej! ) = ej!l . This considers the subband cross correlations using the WSCS(D) autocorrelation of w[n] through
H
E [w[nD + p]w3 [mD + q]] = E vp [n]vq3 [m]
.. .
hT0 (ej! ) .. .
tions. Additionally, through (16) the decimated output MSE is given by
G
(20)
K +1
K
(21)
Realizing (ej! ) = D 01 (ej! ), the result (20) simplifies to the original optimal solution (11), thereby showing equivalence of the solu-
(24)
for all n, m, and p; q = 0; 1; . . . ; C 0 1. For C p; q D 0 1, (24) can be set to zero. This configuration models noise that is introduced subsequent to the sampling process. To further generalize, K distinct interfering sources (w1 [n]; . . . ; wK [n]) can be considered using K + 1 merging banks of analysis filters. This framework is useful for considering generalized multiple-input multiple-output (MIMO) problems, since a separate synthesis bank can be determined for each output. For this general problem, a separate alias component matrix is constructed for each analysis bank. The matrix corresponding to the k th interfering source wk [n] is defined as through (6) and labeled k;AC (ej! ).
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
G
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
Using this generalization, the optimal reconstruction filter bank and error CSD matrix are still determined, respectively, through (11) and (14). However, the component matrices Q(ej! ) and R(ej! ) must be modified to account for the multiple source structure. The new definitions for multiple sources are
Q(ej! ) = HTAC (ej! )Sxx (ej! ) +
K k=1
GTk;AC (ej! )Sxw
(e
T (ej! )SH (ej! )H3 (ej! ) R(ej! ) = HAC xx AC K
+
+
+
l=1 K
T (ej! )SH HAC xw
(e
j! )
k=1 l=1
IV. CONCLUDING REMARKS AND PROBLEM EXTENSIONS
j! )G3 (ej! ) l;AC
GTk;AC (ej! )SwH w (ej! )Gl;3 AC (ej! ):
An alternate numerical solution was presented for this problem in the earlier work of Ohno and Sakai [5]. Their approach examined the error caused by missing subband of a biorthogonal filter bank in the correlation domain. Numerical optimization techniques were then used to find synthesis filters of predetermined length. Although their solution could not guarantee the globally optimal solution, it did have the advantage of allowing specific design constraints (e.g., linear phase synthesis filters).
(25)
3 (ej! ) GTk;AC (ej! )SwH x (ej! )HAC
k=1 K K
3241
(26)
It is worth noting that for a given set of analysis banks modeling, a MIMO scenario R(ej! ) does not change with respect to selection of the desired signal. This is useful for reducing the design complexity since each distinct synthesis bank is determined by applying a unique Q(ej! ) to a common R01 (ej! ).
A filter bank framework for generalized sampling without strict bandlimiting conditions was demonstrated. Optimization techniques were used to determine a set of filters to minimize the TAMSE of reconstruction for a given signal and sampling model. A generalized form of the problem was also presented for cases of additive noise and interfering sources. An example sampling scenario using this presented approach can be found in [27]. There also exist several extensions using the foundation presented in this correspondence. A similar method was used in [28] for signal reconstruction in the case of unknown random sampler timing offsets. A multidimensional approach is used in [29] for image super-resolution. Current research investigating the optimization of periodic sampling patterns is ongoing. APPENDIX CALCULATION OF (11)
C. Numerical Considerations While an analytical solution has been considered so far, a feasible solution requires a numerical approach. For this, the solution (11) is determined for varying frequencies over j! j =D , providing sampled versions of the synthesis filters’ frequency responses through piecewise combinations of the alias components. The sampled frequency responses should be of sufficiently fine resolution from which accurate filter realizations can be determined. A high-resolution model will 2 generally also allow an accurate calculation of e; opt through numerical integration of (9). For the solution (11) to exist at a particular ! , the matrix R(ej! ) of (13) must be invertible. As mentioned in Section II-B, Sxx (ej! ) is diagonal for a WSS x[n], with real nonnegative terms representing j! j! portions of the power spectrum. This means SH xx (e ) = Sxx (e ), allowing the (p; q)th entry of R(ej! ) to be found through
R(ej! )
D01 (
p;q)
=
k=0
Sxx (e
j! W k )H
p (e
j! W k )H 3 (ej! W k ) (27) q
where (p; q) = 0; . . . ; C 0 1. To be invertible, R(ej! ) must be of full rank C , requiring diagonal Sxx (ej! ) to have at least C nonzero entries. This is not true at frequencies for which there are fewer nonzero aliased components than samplers. Although R(ej! ) is not invertible in these cases, the reduced number of aliased components actually indicates oversampled frequencies which can be made free of reconstruction error. To circumvent this problem, the zero-valued portions of Sxx (ej! ) can be replaced with a small positive value (much smaller than the non-zero valued spectral components), thereby providing an invertible approximation of R(ej! ). The unaltered Sxx (ej! ) is used for Q(ej! )(12), the noninverted portion of the solution, ensuring zero-valued portions of the input spectrum remain so at the output. Invertibility of (27) also places some restrictions on HAC (ej! ), generally requiring unique timing offsets or distortions for the samplers. In the event of duplicate sampling, the redundant branches of the model should be removed, corresponding to the removal of linearly dependent columns from HAC (ej! ).
The solution to (10) is found by setting the derivative of See (ej! )) with respect to FAC (ej! ) equal to zero (the relevant operations can be found in [30, App. E]). This produces tr(
3 j! H j! T j! j! FAC ;opt (e )(HAC (e )Sxx (e )HAC (e )) =D
STxx (ej! )HAC (ej! )
from which (11) is obtained provided the inverse of H (ej! )ST (ej! )H (ej! ) exists. HAC AC xx To ensure the derivative’s zero-point is a global minimum, the convexity of tr(See (ej! )) is examined. Ordering the entries of FAC (ej! ) vectorially through the concatenation of its rows, the Levi matrix (or, complex Hessian) [31] of tr(See (ej! )) with respect to FAC (ej! ) is calculated. The result is a block diagonal matrix with identical size j! 2 C 2 C blocks equal to R(e )=D , a scaled version of (13). This structure means convexity is proven if and only if R(ej! ) is positive definite. Although a detailed examination of the necessary conditions for a positive definite R(ej! ) is not provided here, Section III-C discusses some conditions causing a noninvertible R(ej! ), resulting in a semidefinite matrix for which a unique minimum cannot be found.
REFERENCES [1] A. Papoulis, “Generalized sampling expansion,” IEEE Trans. Circuits Syst., vol. CAS-24, no. 11, pp. 652–654, Nov. 1977. [2] P. P. Vaidyanathan, Multirate Systems and Filter Banks. EnglewoodCliffs, NJ: Prentice-Hall, 1993. [3] M. K. Tsatsanis and G. B. Giannakis, “Principal component filter banks for optimal multiresolution analysis,” IEEE Trans. Signal Process., vol. 43, no. 8, pp. 1766–1777, Aug. 1995. [4] R. A. Haddad and K. Park, “Modeling, analysis, and optimum design of quantized M-band filter banks,” IEEE Trans. Signal Process., vol. 43, no. 11, pp. 2540–2549, Nov. 1995. [5] S. Ohno and H. Sakai, “Optimization of filter banks using cyclostationary spectral analysis,” IEEE Trans. Signal Process., vol. 44, no. 11, pp. 2718–2725, Nov. 1996.
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
3242
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 8, AUGUST 2006
[6] K. Gosse and P. Duhamel, “Perfect reconstruction versus MMSE filter banks in source coding,” IEEE Trans. Signal Process., vol. 45, no. 9, pp. 2188–2202, Sep. 1997. [7] P. P. Vaidyanathan, “Theory of optimal orthonormal subband coders,” IEEE Trans. Signal Process., vol. 46, no. 6, pp. 1528–1543, Jun. 1998. [8] P. P. Vaidyanathan and A. Kirac, “Results on optimal biorthogonal filter banks,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 45, no. 8, pp. 932–947, Aug. 1998. [9] J. Tuqan and P. P. Vaidyanathan, “Oversampling PCM techniques and optimum noise shapers for quantizing a class of nonbandlimited signals,” IEEE Trans. Signal Process., vol. 47, no. 2, pp. 389–407, Feb. 1999. [10] P. P. Vaidyanathan and V. C. Liu, “Efficient reconstruction of band-limited sequences from nonuniformly decimated versions by use of polyphase filter banks,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 38, no. 11, pp. 1927–1936, Nov. 1990. [11] A. J. Coulson, “A generalization of nonuniform bandpass sampling,” IEEE Trans. Signal Process., vol. 43, no. 3, pp. 694–704, Mar. 1995. [12] P. P. Vaidyanathan and S.-M. Phoong, “Reconstruction of sequences from nonuniform samples,” in Proc. IEEE Int. Symp. Circuits Systems (ISCAS), Seattle, WA, 1995, vol. 1, pp. 601–604. [13] C. Herley and P. W. Wong, “Minimum rate sampling and reconstruction of signals with arbitrary frequency support,” IEEE Trans. Inf. Theory, vol. 45, no. 5, pp. 1555–1564, Jul. 1999. [14] Y. C. Eldar and A. V. Oppenheim, “Filterbank reconstruction of bandlimited signals from nonuniform and generalized samples,” IEEE Trans. Signal Process., vol. 48, no. 10, pp. 2864–2875, Oct. 2000. [15] P. P. Vaidyanathan, “Generalizations of the sampling theorem: Seven decades after Nyquist,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 48, pp. 1094–1109, Sep. 2001. [16] R. Venkataramani and Y. Bresler, “Optimal sub-Nyquist nonuniform sampling and reconstruction for multiband signals,” IEEE Trans. Signal Process., vol. 49, no. 10, pp. 2301–2313, Oct. 2001. [17] H. Johansson and P. Löwenborg, “Reconstruction of nonuniformly sampled bandlimited signals by means of digital fractional delay filters,” IEEE Trans. Signal Process., vol. 50, no. 11, pp. 2757–2767, Nov. 2002. [18] R. S. Prendergast, B. C. Levy, and P. J. Hurst, “Reconstruction of bandlimited periodic nonuniformly sampled signals through multirate filter banks,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 51, pp. 1612–1622, Aug. 2004. [19] A. S. Bopardikar, R. M. Rao, and B. S. Adiga, “Matched sampling systems, relation to wavelets and implementation using PRCC filter banks,” IEEE Trans. Signal Process., vol. 48, no. 8, pp. 2269–2278, Aug. 2000. [20] D. Seidner and M. Feder, “Vector sampling expansion,” IEEE Trans. Signal Process., vol. 48, no. 5, pp. 1401–1416, May 2000. [21] R. Venkataramani and Y. Bresler, “Filter design for MIMO sampling and reconstruction,” IEEE Trans. Signal Process., vol. 51, no. 12, pp. 3164–3176, Dec. 2003. [22] M. Unser and J. Zerubia, “Generalized sampling: Stability and performance analysis,” IEEE Trans. Signal Process., vol. 45, no. 12, pp. 2941–2950, Dec. 1997. [23] ——, “A generalized sampling theory without band-limiting constraints,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 45, no. 8, pp. 959–969, Aug. 1998. [24] M. Unser, “Sampling-50 years after Shannon,” Proc. IEEE, vol. 88, pp. 569–587, Apr. 2000. [25] T. I. Laakso, V. Välimäki, M. Karjalainen, and U. K. Laine, “Splitting the unit delay,” IEEE Signal Process. Mag., vol. 13, no. 1, pp. 30–60, Jan. 1996. [26] W. A. Gardner, Introduction to Random Processes with Applications to Signals and Systems. New York: McGraw-Hill, 1989. [27] R. S. Prendergast and T. Q. Nguyen, “Optimal filter bank reconstruction of periodically undersampled signals,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing (ICASSP), Philadelphia, PA, 2005, vol. 4, pp. 201–204.
[28] ——, “Optimal reconstruction of periodically sampled signals with probabilistic timing delays,” presented at the 38th Asilomar Conf. Signals, Systems, Computers Pacific Grove, CA, Nov. 2004. [29] ——, “Improving frequency domain super-resolution via undersampling model,” in IEEE Int. Conf. Image Processing, 2005, pp. 853–856. [30] T. K. Moon and W. C. Stirling, Mathematical Methods and Algorithms for Signal Processing. Upper Saddle River, NJ: Prentice-Hall, 2000. [31] S. G. Krantz, Function Theory of Several Complex Variables. Providence, RI: AMS Chelsea, 1992.
Optimum Blind Multichannel Equalization Using the Linear Prediction Algorithm Houcem Gazzah, Member, IEEE
Abstract—The linear prediction (LP) algorithm estimates a zero-forcing (ZF) equalizer from the channel output’s second-order statistics (SOS). We show that it allows, in fact, to determine a whole subspace of (but not all) ZF equalizers. The one obtained by maximizing the signal-to-noise ratio (SNR) at the equalizer output not only outperforms the original LP equalizer, but also attains the lowest achievable (by any no-delay ZF equalizer) mean-square error (MSE). The increase in complexity is only marginal, but equalization performance can be enhanced significantly, especially in the (practical) case of a small leading channel coefficient. Index Terms—Channel modeling, equalization, estimation.
I. INTRODUCTION Training sequences are becoming inadequate with the growing demand on bandwidth. This has motivated an intensive interest in blind techniques, and more particularly, those based on (the efficiently estimated) second-order statistics (SOS). This is made possible if the channel output is oversampled and/or received using multiple sensors. The linear prediction (LP) algorithm [1], [2] is among the most popular SOS-based approaches. The LP equalizer is obtained by linearly combining the columns of the predictor. Each, in fact, can be regarded as a zero-forcing (ZF) equalizer, and each has a different behavior in the presence of observation noise. This fact was not considered in the original LP algorithm, and hence, the computed equalizer is not optimal. In this correspondence, we optimize the LP algorithm using the signal-to-noise ratio (SNR) criterion [3], [4] (Section IV). We show the so-optimized LP equalizer to be nothing but the ZF minimum meansquare error (MMSE) zero-delay equalizer. The gain of the optimized solution is especially relevant in the practical case of a small leading channel coefficient. A performance analysis (Section V) and simulations (Section VI) support these claims. We denote by Diag(a1 ; . . .) the diagonal matrix whose ith diagonal entry is ai , by Ia the a 2 a identity matrix and by 0a;b the a 2 b zero matrix. Indexes are omitted whenever the dimensions can be inferred from the context. The ith row, j th column entry of the matrix Manuscript received February 15, 2005; revised October 9, 2005. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Helmut Boelcskei. The author is with the Institute for Digital Communications, University of Edinburgh, Edinburgh, EH9 3JL U.K. (e-mail: [email protected]). Digital Object Identifier 10.1109/TSP.2006.877669
1053-587X/$20.00 © 2006 IEEE
Authorized licensed use limited to: IEEE Xplore. Downloaded on February 27, 2009 at 13:51 from IEEE Xplore. Restrictions apply.
