Page  00000001 Subband HTLS Tom Schouten, Bart De Moor Department ESAT/SISTA, KULeuven email: Abstract In this paper we present a reparametrization of an autoregressive modeling problem using a subband or interleaving representation. This approach increases the numerical condition of the problem in the case where the AR order is small and there are closely spaced poles or poles close to the unit circle. The side effect is the introduction of aliasing. We present a method that can undo the aliasing and classify detected poles as part of the noise or signal model. 1 Introduction 1.1 Resolution enhancement By using data interleaving on an oversampled signal one can enhance the accuracy and MA noise robustness of an AR modeling method, compared to the same algorithm applied to the undecimated signal. See (?) (?) for more information on this subject. Mere oversampling increases the data size and thus increases statistical accuracy but reduces the condition of the problem, hence introduces numerical inaccuracies. When interleaving of the oversampled signal is applied, the excess of data can be used to increase statistical precision without deteriorating the numerical condition. This idea can be applied to increase resolution for parameter estimation on signals which can be modeled as an AR process with groups of closely spaced poles. When one constructs a criticly sampled subband representation of the original signal (i.e. an interleaving scheme) and one uses all subbands in the calculation of the signal poles, the numerical condition of the problem improves because the subsampling operation moves the poles away from the unit circle and other poles. The only disadvantage is the introduction of aliasing due to the subsampling involved. So in short, instead of increasing statistical accuracy by using more data (oversampling) and rearranging it (interleaving), as described in the cited papers, we use the same amount of data but change the parameterization (critically sampled subband scheme) to increase the numerical condition of a certain class of problems. The method is analogous the the cited papers apart from the fact that the original signal is not oversampled and needs 'anti-aliasing' of the detected poles based on amplitude information in different subbands. In this report we present several approaches using different subband schemes. Section two dicusses a criticly sampled sliding FFT subband scheme. Section three discusses interleaving. Section four discusses other possible subband schemes. 1.2 The HTLS and derived methods Although the algorithm presented can be employed using other multi channel AR modeling algorithms, we use the Hankel Total Least Squares method (HTLS) and the multichannel extensions HTLSsum and HTLSstack. This method uses a subspace approach to extract signal poles (frequency, damping) and a linear least squares algorithm to estimate phase and amplitude. For more information on these algorithms we refer to (?) (?). 2 Sliding FFT subband method 2.1 Introduction First we will take a look at the representation of a single complex exponential in a DFT subband representation. For now we will limit ourselves to an equally weighted (square window), criticly sampled sliding FFT representation. To prevent as much aliasing as possible we might apply a non-rectangular window later on. We consider a signal s[n] = cpn with c e C, p = ed+jw E C and n e N. We represent this signal in a criticly sampled N-band representation with {Sm[k]m = 0,1,...,N - 1} N-i Sm[k] = e-j27rmn/Ns[n + kN] n=O (1) (2)

Page  00000002 and k e N. After some algebraic manipulations we find the following expression for the subband signals Sm[k] = Am(d, W)CpkN (3) where 1 pN SPm pm = ed+j(w-27rm/N) (4) 1. Construct an M-channel AR model for M subbands using HTLSStack or HTLSSum. This produces the (aliased) poles qk and the subband amplitudes Cm,k corresponding to qk. 2. Obtain the winding factor(s) and the amplitude(s) for each signal pole pkl = q/N W1 using the information present in the detected subband amplitudes Cm,k for each detected aliased pole qk -= P 2.2 Illustration: recovery of a single pole in the abscence of noise To illustrate the problem we look at the following toy problem. The signal consists of 1 pole p with amplitude c. (5) yields the following relationship between the aliased pole q = pN, the integer winding factor 1 defined by p = q/INW1, the detected subband amplitudes Cm for pole q and the signal pole amplitude c represents the amplitude and phase distortion of the original complex exponential. When we write this in terms of p and use W = ej2,/N we get 1-p Am(p) = 1-p -. (5) 1 - pW-m (3) is again a complex exponential. The amplitude and phase are modified by (5) which depends on signal parameters d, w and subband scheme parameters m, N. The pole of the signal is, aside from the decimation power N, unaffected by the subband scheme. If we want to use subband scheme (1) for detection of the original pole p, we have to undo the aliasing and the amplitude effect. How can we manage this? Fist let us inspect the complex amplitude factor (5). Non-zero values of m introduce a non-zero shift. Increasing Idl introduces spread and thus increases aliasing into neighboring bands. The maximum is reached for p = Wm and is equal to N. The function is zero for p = W1 with (1 - m) mod N 0. Aliasing introduces ambiguity with respect to the exact location of a detected pole. Next we will describe how to undo the aliasing by using the information present in the amplitudes of different realizations of the same signal pole in neighboring subbands. If we would use mere subsampling, there would only be a phase difference between the different realizations. The filtering present in the subband scheme introduces a phase and amplitude difference. In theory, the order of each subband signal is equal to the order of the original signal. Exeptions are poles that collide with one of the base functions of the subband scheme p = W1 with 1 e Z. They are present in only one subband and are removed from the others, reducing their order compared to the original signal. Since this case is very rare, and overestimation of the order is generally not a problem, we can assume the order of each subband to be equal to the order of the signal. Next we explore the possibility to perform anti-aliasing in the exact case by using the information given by the M subbands of the subband representation. The procedure is as follows. Cm C= c 1-q 1 - ql/NW-m (6) The unknowns in these equations are c and 1. This equation is nonlinear in 1. It is possible to write the equations linear in c and W1. One needs only two subband amplitudes to obtain c and 1. However, numerical experiments have indicated that it is better to expand this equation to a multi-mode version in the case of noise on the Cm and to be able to resolve collision problems due to aliasing. 2.3 Moving to multiple poles: some remarks 2.3.1 Maximal order of the original signal There is a maximal model order limited by the size of the HTLSStack matrix, which is in turn defined by the length of the subband signals. If the signal order exceeds this limit, the method breaks down as with ordinary HTLS methods. The maximal order is limited by half of the subband signal length. 2.3.2 Collision problem A problem can arise when the signal has multiple poles. It is possible they collide after subsampling. Suppose there are 2 poles pi, p2 for which we have pj = p=, with pi # p2. A possible fix for this problem is the use of a different decimation factor. We will investigate this case using 2 subband schemes with decimation factors N and N', with N / N'. Suppose pi = pWki and p2 = pWk2, with ki,k2 E Z, kI / k2 and im(logp) E [0, 2r/N]. Is it possible to

Page  00000003 have p;' = pj' in the second subband scheme? Suppose W' = ej2r/N', pi = p'W'kl and P2 = p'W'k2, with k[,ki e Z, k[ $ ki and im(logp') E [0,27r/N]. We then have pi/Pj = Wki-k2 = Wikl-k2 or N'(k - k2) = N(k' - ks). This is only possible when N and N' are not coprime. When they are, two poles that will collide in one subband scheme, will not collide in the other. We can use this to avoid the collision problem, but how can we detect it? Experiments show that the detection of a collision is more difficult than the modeling of all possible collisions. We will introduce an exact method we dubbed linear multimodal fit (LMF) that performs this operation. In short the LMF approach needs only one subband scheme. We can construct a linear multimodal model for all the possible parents of q, being ql/N W. We can solve for the amplitudes ci in an exact manner. The data dependencies are as follows. First we obtain q, then for each q we obtain the Cm, the subband amplitude corresponding to pole q in band m. Out of the Cm we compute an amplitude factor ci for each possible parent of q. Then we need to analyse the obtained spectrum to separate signal from noise components. 2.4 Linear multimodal fit (LMF) We have the following relationship between the subband amplitudes Cm corresponding to the aliased pole q if we include the effect of all possible parents of q being ql/NW1. N-1l m ck l/N -m (7) k=0 1 If we define tl = 1-ql-w- (7) transforms to tk = v/Nq(-k)modN (11) 2.5 Signal/Noise separation The LMF method described above produces a spectrum of amplitudes ci of pole pi = ql/NWI, corresponding to one of the roots of q = pN. Because this is an exact fit, noise on the Cm will effect the c1. We still need a method to determine if c1 can be seen as signal or noise. We do this using a threshold method. Some remarks. If q is a noise pole (or ghost pole), the Cm and ck will be fairly uncorrelated because q is then an artifact of the HTLSstack/HTLSsum procedure and will have uncorrelated amplitudes in all the subsampled sequences Sm. If q is a signal pole, the reverse is true. How can we use this information to remove both HTLSstack/HTLSsum artifacts (spurious noise poles) and anti-aliasing procedure artifacts (poles that collide after subsampling)? A possible approach is the use of variance estimation and thresholding. We use the universal threshold, being a/2 log N, proposed in wavelet denoising literature. It has several properties such that it can be considered as optimal. We use a variance estimator that is robust with respect to outliers a = MAD/.6745, with MAD the average of the absolute values of the vector. Experimental results show that this kind of thresholding rule performs pretty well. 2.6 The algorithm 1. Construct subband scheme. Divide the original signal s into N subbands Sm. 2. Obtain aliased poles. Calculate K (aliased) poles qk using HTLSstack or HTLSsum. 3. Obtain pseudoinverse. Construct a pseuodoinverse for the LS problem Sm[n] + A Sm[n] = --o Cm,kqk. 4. Obtain subband amplitudes. For each subband m, use the pseudoinverse to calculate Cm,k, the subband amplitudes of the qk. 5. Calculate LMF model. For each qk, obtain ck,l solving (7). 6. Perform anti-aliasing and signal/noise separation. For k = 1: K - 1, calculate ak and retain all Ck,l: iCk,l] > ak2logN. N-1 Cm = Cktm-k k=O (8) or C =Tc (9) in matrix notation. T is circulant because we have the property tl = tl+N. (9) can thus be easily solved using an FFT algorithm. If we define t = (to, t,..., tN--1) to be the first column of T and F the normalized DFT matrix, we have T = Ftdiag(Ft)F. The solution of (9) is then given by c = Ftdiag-l(Ft)FC. (10) We can express diag(Ft)-1 as a function of q without having to compute t and Ft. Suppose t' = Ft. The elements of t' are given by

Page  00000004 3 Interleaving subband scheme Instead of using a subband scheme, we can use interleaving as an alternative. It should still be possible to undo the aliasing process using the information present in the phase shift due to interleaving, as opposed to amplitude and phase altering in the subband scheme. 3.1 Introduction First we present a method to find the winding factor when there is no noise present and there are no colliding poles. If we write the analog of (2) without the filtering we get the following expression. 3.3 Noise robustness For the separation of noise and signal poles we use the same method as presented in the previous section. 4 Windowed FFT The possibility to construct a windowed FFT version to further reduce aliasing instead of just applying the aliasing recovery operation described above was examined. The filterbank representation with a window w[n] applied looks like N-i Sm[k] = e-j27mn/Nw[n]s[n + kN] n=O (19) Sm [k] = s[m+ kN] = cpmpkN The analog for (5) is Am(p) = pm. Given the complex amplitudes Cm, c recovered out of the equations (12) and p ca Cm = cpm As we did in the previous sections, we can expres equations in the variables c and 1, using q as a kr parameter. (14) then becomes Cm = c(ql/NWl)m, which is the analog of (6). The next paragraphs discuss the influence of noise and colliding modes. 3.2 LMF Changing the filterbank characteristic will of course alter (5). The question is, will a proper choice for w[n] still produce a tractable formula to recover the wind(13) ing factor and original amplitude? And what about the S numerical condition of the problem? n be It is still possible to use the techniques presented above. This is because the equations stay linear if we (14) impose a model that includes all the possible aliased modes. The set matrix of this equation is still circus the lar, so the procedure will be fast. We will show this in aown the following paragraphs. Suppose a sliding FFT subband scheme with weighting window w[n]. The amplitude response will be, with(15) out loss of generality, Am(p), a (nonlinear) function of will p. The multimodal equation then becomes 3 V I111 N-i Cm = CkAm(qI/NWk). k=0 (20) Again, as in (7), we write Cm as a sum of contributions from all possible parents of q being zk = ql/NWk. This set is square and can be solved in an exact manner. In short It is easy to veryfy that Am(p) has the following property N--1 Cm = ck(ql/NWk)m. k=O If we multiply (16) by q-m/N we get N-1 Cq-m/N = CkWkm. k=0 The Ck are then given by 1Nck = NT Cmq-m W-km m=O This can be computed using the FFT. (16) Am(p) = Am+i(pW) = Ao(pW -m) (20) then becomes N-i Cm =I ckAo(ql/NWk-m), k=0 (21) (22) or C = Tc, which is clearly a circulant set. If we (17) set ti = Ao(ql/NW-l) we get again (8). If we define t = (to, ti,.., tN-1)T to be the first column of T and F the normalized DFT matrix, we have T = Ftdiag(Ft)F. The solution of (9) is then given by (18) c = Ftdiag-l(Ft)FC. (23)

Page  00000005 This is generally applicable to all windowing setups for which we can compute Ao(p). To speed up the procedure, it would help to have a closed form expression for Ao (p), or even better for Ft, as we did above. In general we have N-i Ao(p) = Spkw[k]. k=0 The expression for t1 becomes (24) N-i ti = I (q/NW-l)kw[k]. k=0 (25) If we define wq[k] = qk/Nww[k], a modulated version of the window, we get t = Fwq. Then we have Ft = FFw, which is a mirrored version of wq. More precisely (Ft)[k] = wq[-k mod N]. This general expression for Ft allows us to draw some conclusions about (23). Because of the inverse operation on diag(Ft), coefficients near 0 will be dominant in the equation. This can happen when q is relatively close to the origin or when w[k] approaches zero for some k. What will be the effect of these two influences on the sparseness of c? Experiments indicate that sparseness breaks down when q is close to the origin in the presence of noise. We will investigate this. Supose C is corrupted by zero mean gaussian noise AC in (23). The expression for Ac then becomes Ac = Ftdiag- (FFwq)FAC. (26) Because (FFwq)[0] = 1 the DC part is unnaffected. We then have E Ac = Z AC. If q is on the unit circle, there is only a phase effect. If q| < 1 or I q > 1 this will change the colour of the noise. This phenomenon troubles the threshold based separation of signal and noise and can produce spurious estimates. That is why it is best to limit the distance to the unit circle for the qi to a certain safety margin. 5 Examples The results of Figure 1 were obtained by using an interleaving scheme to estimate the parameters of a 65536 -sample signal consisting of a sum of 10 slightly damped, low frequency complex exponentials and added gaussian noise. The SNR in this example is -1dB. The signal was interleaved to give 256 bands of 256-sample signals. The figure shows only the real part. Matlab scripts are available upon request. x104 Figure 1: The top graph shows the original signal, the middle graph shows the signal with noise added and the bottom graph shows the reconstructed signal. 6 Conclusions This paper describes a method to tackle the problem of estimating a low order AR model of a signal with closely spaced poles or poles close to the unit circle. It does so by rearranging the signal in a subband or interleaving representation, applying a multichannel AR estimation step and combining these results to obtain the disired AR parameters. The benefits are better numerical condition, significant speedup and the ability to distinguish between signal and noise contributions. References [1] B. Halder, T. Kailath, Efficient estimation of closely spaced sinusoidal frequencies using subspace-based methods, IEEE Signal Processing Letters, Vol. 4, No. 2, Feb 1997, 49-51. [2] P. Stoica, A. E. Nordsji, Subspace-based frequency estimation in the presence of moving-average noise using decimation, Signal Processing 63 (1997), 211-220. [3] L. Vanhamme, Advanced time-domain methods for nuclear magnetic resonance spectroscopy data analysis, PhD thesis, Faculty of Engineering, K.U.Leuven (Leuven, Belgium), K.U.Leuven, Nov. 1999, 210 p. [4] L. Vanhamme, S. Van Huffel, Multichannel quantification of biomedical magnetic resonance spectroscopy signals, in Advanced Signal Processing Algorithms, Architectures, and Implementations VIII, (Luk F.T., ed.), Proc. of SPIE's 43rd Annual Meeting (SPIE), San Diego, California, Jul. 1998, vol. 3461 of SPIE Proceedings series, SPIE-the International Society for Optical Engineering (Washington, USA), 1998, pp. 237-248.