Analysis of Signals via Non-Maximally Decimated Non-Uniform Filter Banks

TL;DR

This paper develops a theoretical framework for designing optimal FIR synthesis filter banks for signal reconstruction from multirate observations, including non-uniform and non-maximally decimated cases, with conditions for perfect reconstruction.

Contribution

It introduces a novel theoretical approach with pseudocirculant conditions for perfect reconstruction in non-maximally decimated non-uniform filter banks.

Findings

Derived necessary and sufficient conditions for perfect reconstruction.

Provided a parametric form for all optimal solutions.

Verified results through simulations.

Abstract

This paper addresses the important problem of reconstructing a signal from multiple multirate observations. The observations are modeled as the output of an analysis bank, and time-domain analysis is carried out to design an optimal FIR synthesis bank. We pose this as a minimizing the mean-square problem and prove that at least one optimal solution is always possible. A parametric form for all optimal solutions is obtained for a non-maximally decimated filter bank. The necessary and sufficient conditions for an optimal solution, that results in perfect reconstruction (PR), are derived as time-domain pseudocirculant conditions. This represents a novel theoretical contribution in multirate filter bank theory. We explore PR in a more general setting. This results in the ability to design a synthesis bank with a particular delay in the reconstruction. Using this delay, one can achieve PR in cases where it might not have been possible otherwise. Further, we extend the design and analysis to non-uniform filter banks and carry out simulations to verify the derived results.

Tables

Table 1. TABLE I: Experiment 1: MSE for various delay values

$d$	$J_{\text{min}}^0$(dB)	$J_{\text{min}}^1$(dB)	$J_{\text{min}}$(dB)
0	-0.8673	-1.3434	1.9115
5	-33.1717	-43.7244	-32.8052
9	-62.8689	-68.3970	-61.7967
10	-68.3970	-70.7269	-66.3972
11	-70.7269	-61.6691	-61.1606
15	-57.3320	-45.5881	-45.3068
20	-29.5288	-33.7826	-28.1442

Table 2. TABLE II: Experiment 2: Variation in MSE with number of available subbands

$L$	Available analysis filters	$d$	$J_{\text{min}}^0$(dB)	$J_{\text{min}}^1$(dB)	$J_{\text{min}}^2$(dB)	$J_{\text{min}}^3$(dB)	$J_{\text{min}}$(dB)
1	$h_0\left(n\right)$	10	-12.9106	-12.9711	-15.0069	-14.9634	-7.8231
1	$h_1\left(n\right)$	10	-0.2141	-0.2360	-0.0490	-0.1161	5.8675
1	$h_2\left(n\right)$	10	-0.0826	-0.0577	-0.1747	-0.0823	5.9215
1	$h_3\left(n\right)$	10	-0.0865	-0.0499	-0.0612	-0.0241	5.9652
2	$h_0\left(n\right)$, $h_1\left(n\right)$	10	-16.7631	-19.0161	-16.7292	-19.1518	-11.7388
2	$h_0\left(n\right)$, $h_2\left(n\right)$	10	-14.3354	-13.1462	-18.6208	-15.8517	-9.0295
2	$h_0\left(n\right)$, $h_3\left(n\right)$	10	-13.7442	-13.8427	-15.6931	-15.6224	-8.6055
2	$h_1\left(n\right)$, $h_2\left(n\right)$	10	-0.2899	-0.2923	-0.2217	-0.2139	5.7663
2	$h_1\left(n\right)$, $h_3\left(n\right)$	10	-0.3044	-0.2881	-0.1109	-0.1410	5.8104
2	$h_2\left(n\right)$, $h_3\left(n\right)$	10	-0.1676	-0.1071	-0.2355	-0.1088	5.8662
3	$h_0\left(n\right)$, $h_1\left(n\right)$, $h_2\left(n\right)$	10	-20.5040	-20.3481	-23.1342	-23.1368	-15.5518
3	$h_0\left(n\right)$, $h_1\left(n\right)$, $h_3\left(n\right)$	10	-19.2065	-24.3326	-17.8184	-21.2263	-14.0058
3	$h_0\left(n\right)$, $h_2\left(n\right)$, $h_3\left(n\right)$	10	-15.5108	-14.0397	-20.4628	-16.7167	-10.0944
3	$h_1\left(n\right)$, $h_2\left(n\right)$, $h_3\left(n\right)$	10	-0.3786	-0.3432	-0.2834	-0.2416	5.7092
4	$h_0\left(n\right)$, $h_1\left(n\right)$, $h_2\left(n\right)$, $h_3\left(n\right)$	10	-45.1648	-32.0844	-145.2323	-143.7481	-31.8758
4	$h_0\left(n\right)$, $h_1\left(n\right)$, $h_2\left(n\right)$, $h_3\left(n\right)$	11	-32.0844	-145.2323	-143.7481	-143.5253	-32.0844
4	$h_0\left(n\right)$, $h_1\left(n\right)$, $h_2\left(n\right)$, $h_3\left(n\right)$	12	-145.2323	-143.7481	-143.5253	-145.3962	-138.3732
4	$h_0\left(n\right)$, $h_1\left(n\right)$, $h_2\left(n\right)$, $h_3\left(n\right)$	13	-143.7481	-143.5253	-145.3962	-30.2194	-30.2194

Equations

$$\mathbf{v}(n)={\left[\begin{array}{cccc}v(n-P+1)& v(n-P+2)& \dots & v(n)\end{array}\right]}^T.$$

$$\overline{\mathbf{v}}(n)={\mathbf{J}}_P\mathbf{v}(n),$$

$$\overline{\mathbf{B}}={\mathbf{J}}_M\mathbf{B}{\mathbf{J}}_N$$

$$\overline{\mathbf{AB}}=\overline{\mathbf{A}}\overline{\mathbf{B}}.$$

$$y(n)=x(Mn+l).$$

$$\mathbf{y}={\left[\begin{array}{cccc}y(0)& y(1)& \dots & y(P-1)\end{array}\right]}^T$$

$$\begin{array}{rl}\left[\begin{array}{c}y(0)\\ y(1)\\ ⋮\\ y(P-1)\end{array}\right]& =\left[\begin{array}{cccccc}1& 0& \dots & 0& \dots & 0\\ 0& 0& \dots & 1& \dots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \dots & 0& \dots & 1\end{array}\right]\left[\begin{array}{c}x(l)\\ x(l+1)\\ ⋮\\ x(M+l)\\ ⋮\\ x(M(P-1)+l)\end{array}\right],\\ \mathbf{y}& ={\mathbf{D}}_{P,M(P-1)+1}^0\mathbf{x},\end{array}$$

$$\mathbf{x}={\left[\begin{array}{cccc}x(0)& x(1)& \dots & x(P-1)\end{array}\right]}^T$$

$$\mathbf{y}={\mathbf{D}}_{\lceil \frac{P}{M}\rceil ,P}^l\mathbf{x},$$

$$\begin{array}{rl}{D}_{\lceil \frac{P}{M}\rceil ,P}^l(i,j)& =\{\begin{array}{ll}1,& \text{if }j=Mi+l\\ 0,& \text{otherwise}.\end{array}\end{array}$$

$$\mathbf{x}={\left[\begin{array}{cccc}x(0)& x(1)& \dots & x(P-1)\end{array}\right]}^T$$

$$\begin{array}{rl}\mathbf{y}& ={\mathbf{U}}_{M(P-1)+1,P}\mathbf{x},\end{array}$$

$${\mathbf{U}}_{M(P-1)+1,P}=({\mathbf{D}}_{P,M(P-1)+1}^0{)}^T.$$

$$\begin{array}{rl}{\mathbf{D}}_{P,M(P-1)+1}^0{\mathbf{U}}_{M(P-1)+1,P}& ={\mathbf{I}}_P,\\ {\mathbf{U}}_{M(P-1)+1,P}{\mathbf{D}}_{P,M(P-1)+1}^0& \ne {\mathbf{I}}_{M(P-1)+1}.\end{array}$$

$$y(n)=\sum _{k=0}^{Q-1}h(k)x(n-k).$$

$$\mathbf{y}(n)={\left[\begin{array}{cccc}y(n-P+1)& y(n-P+2)& \dots & y(n)\end{array}\right]}^T$$

$$\mathbf{y}(n)=\mathbf{Hx}(n),$$

$$\mathbf{x}(n)={\left[\begin{array}{cccc}x(n-P+2-Q)& x(n-P+3-Q)& \dots & x(n)\end{array}\right]}^T,$$

$$\mathbf{H}=\left[\begin{array}{cccccc}h(Q-1)& h(Q-2)& \dots & h(0)& \dots & 0\\ 0& h(Q-1)& \dots & h(1)& \dots & 0\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& 0& \dots & 0& \dots & h(0)\end{array}\right].$$

$$\mathbf{x}={\left[\begin{array}{cccc}x(0)& x(1)& \dots & x(P-1)\end{array}\right]}^T,$$

$$\left[\begin{array}{c}y(0)\\ y(1)\\ ⋮\\ y(Q-1)\\ ⋮\\ y(P+Q-2)\end{array}\right]=\left[\begin{array}{cccc}h(0)& 0& \dots & 0\\ h(1)& h(0)& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ h(Q-1)& h(Q-2)& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & h(Q-1)\end{array}\right]\left[\begin{array}{c}x(0)\\ x(1)\\ ⋮\\ x(P-1)\end{array}\right]$$

$$\mathbf{y}=\tilde{\mathbf{H}}\mathbf{x}.$$

$$\tilde{\mathbf{H}}=\overline{\mathbf{H}}{}^T,$$

$${\mathbf{r}}_{d_{i}v}={\mathbf{a}}_i^T{\mathbf{R}}_{vv},0\le i\le M-1,$$

$${\mathbf{v}}_{\text{s}}(n)={\left[\begin{array}{cccc}\overline{\mathbf{v}}{}_0^T(n)& \overline{\mathbf{v}}{}_1^T(n)& \dots & \overline{\mathbf{v}}{}_{L-1}^T(n)\end{array}\right]}^T,$$

$$\begin{array}{rl}{\mathbf{v}}_j(n)=& {\left[\begin{array}{cccc}{v}_j(n-P+1)& v_j(n-P+2)& \dots & v_j(n)\end{array}\right]}^T.\end{array}$$

$$\begin{array}{rl}{\mathbf{R}}_{vv}& =\mathbf{E}[{\mathbf{v}}_{\text{s}}(n){\mathbf{v}}_{\text{s}}^H(n)],\\ {\mathbf{r}}_{d_{i}v}& =\mathbf{E}[d_i(n){\mathbf{v}}_{\text{s}}^H(n)].\end{array}$$

$$\begin{array}{rl}{\mathbf{v}}_j(n)& ={\mathbf{D}}_{P,M(P-1)+1}^0{\mathbf{x}}_j(Mn),\end{array}$$

$${\mathbf{x}}_j(Mn)={\left[\begin{array}{cccc}{x}_j(M(n-P+1))& x_j(M(n-P+1)+1)& \dots & x_j(Mn)\end{array}\right]}^T.$$

$${\mathbf{x}}_j(Mn)={\mathbf{H}}_j\mathbf{u}(Mn),$$

Full text

I Introduction

Multirate filter banks have been used extensively for applications in the fields of image and video processing [1, 2, 3], communications [4, 5], subband adaptive filtering [6], radar systems [7], etc. The filter bank theory evolved from the need to decompose a signal spectrum into multiple frequency subbands in order to efficiently perform a number of operations like compression, coding, noise cancellation, etc. Initially, two channel quadrature mirror filter bank (QMF) was proposed in the literature [8], and later, it was extended to arbitrary $M$-channel case [9, 10]. To better exploit the signal characteristics in a number of applications like audio coding, spectral analysis, subband adaptive filtering, etc., non-uniform filter banks (NUFBs) having different channel bandwidth were introduced [11, 12, 13]. In [14], filter banks with different multirate factors were used, to realize the filter transfer functions for reducing the multiplicative complexity. NUFBs with rational sampling factors were considered in [12, 15]. However, they were found to be equivalent to integer decimated NUFBs [16]. A number of design methods like direct method [15], recombination method [11], modulation based method [17], etc. exist in the literature. Application specific design methods have also been proposed in the literature. For example, in [18], an NUFB was designed for acoustic echo-cancellation by maximizing the signal-to-alias ratio. All these methods design both the analysis and synthesis banks of an NUFB. The analysis and synthesis filters in a filter bank can be IIR or FIR. FIR filters have advantage over IIR in respect of stability and linear-phase property. In the literature, various design methods using transform like Lapped Orthogonal Transform [19], Generalized Lapped Orthogonal Transform [20], Extended Lapped Transform [21], etc. have been proposed to design a perfect reconstruction (PR) FIR maximally decimated uniform filter bank (UFB). A comprehensive analysis of FIR perfect reconstruction filter banks can be found in [22]. It has been shown that oversampled filter banks (OFBs) offer additional design freedom as compared to critically sampled filter banks [23, 24]. The inbuilt redundancy present in them has been exploited for a variety of applications [25, 26, 27, 28]. A frame theoretic analysis of OFBs can be found in [29]. Given an oversampled analysis bank, a parameterization for all synthesis banks resulting in PR was developed in [30]. Undersampled filter banks which provide underdetermined representation of signals have found application in compression of signals [31]. Despite the filter banks have been around for more than two decades, they continue to be applied to new and emerging areas. Graph signal processing is one such interesting area where filter banks techniques find applications in solving diverse problems [32].

The necessity in any scenario to estimate a signal or its properties from its subband components can be addressed well using the filter bank based processing. Expectedly, many research problems appear posing such requirement. Synthesis from subband components has been utilized in applications like subband coding [33, 34], analog/digital conversion [35], transmultiplexer [33], ultrawideband receiver [36], TV white space transceiver [37], channel coding [38], etc. To get better results, processing on subband components, rather than on full-band signal, have been carried out in problems like speech enhancement [39], ECG processing [40], subband adaptive filtering [41], seismic signals processing [42], etc. The processed subband signals are then combined using a synthesis bank to get the final output. In multirate sensors networks, sensor outputs are modeled as an analysis bank output, and these outputs are fused to get a unified measurement [43]. In [44], construction of a high-resolution signal from multiple low-resolution signals has been carried out. The problem of linear prediction from multirate observations has been examined in [45]. In many applications, rather than a signal, a property of interest is to be estimated from multirate observations. Time delay estimation from low rate samples is an example of that, and it finds application in speech localization, processing using microphone arrays, radar, underwater acoustics, wireless communications and other areas [46, 47]. Direction of arrival estimation is another example which finds application in wireless communications, medical imaging, seismology, etc. [48]. In [31], estimation of the spectrum of a signal from its low resolution signals was carried out. This is required when a signal can be observed only indirectly using low-rate sensors.

To address the problem of estimating a signal from its subbands, various solutions using Wiener filtering, Kalman filtering, $\mathcal{H}^{\infty}$ filtering, compressive sensing, adaptive filtering, etc. have been proposed in the literature. Wiener filtering based methods require joint stationarity of the subband input and desired signals along with the knowledge of second order statistics of them. Kuchler and Therrien [49] used Wiener filtering to estimate a signal from multiple multirate observations. Scrofani and Therrien [44] extended the analysis to two-dimensional signals. Vaidyanathan and Chen [50] designed a statistically optimal synthesis bank for a subband coder by using the biorthogonal solution followed by a matrix Wiener filter in the synthesis stage. In case of Kalman filtering, the statistics of the input and noise signals are required to be known even if that vary with time. Multirate Kalman filter based methods can be found in [51, 52, 53]. $\mathcal{H}^{\infty}$ optimization based methods, on the other hand, do not require any statistical assumption and perform worst-case design [54]. As $\mathcal{H}^{\infty}$ solution is generally not unique, mixed $\mathcal{H}^{2}/\mathcal{H}^{\infty}$ optimization methods [55, 56] have also been proposed in the literature. The complexity of $\mathcal{H}^{2}$ and $\mathcal{H}^{2}/\mathcal{H}^{\infty}$ solutions are very high since they are solved as linear matrix inequalities (LMIs) [56]. A least square based solution for the problem was presented in [57]. However, the solution presented was for two channels only. In [58], Li Chai et al. formulated the problem using structural similarity criterion and obtained closed-form solution under certain assumptions like the filter length is equal to the decimation factor, the mean of the source is zero, etc. Compresive sensing based method was used in [59] to estimate a signal from multirate observations. However, the signal has to be sparse. In [60, 61], adaptive filters were used to estimate a signal from two observation sequences. An adaptive least square lattice filter was used in [62]. However, the filter obtained was periodically time-varying. Statistical analysis of the estimation from the subbands has been performed by a number of authors. For example, in [63, 64, 65], power spectral density of a signal was estimated from multiple observations. In [66], a measure for information contained in a multirate observations was defined.

I-A Contribution

In this paper, the problem of synthesizing a signal from its multirate observations is addressed by performing time-domain analysis. This approach has led to some additional novel contributions to the field which are now briefly presented. Under time-domain analysis, the matrix representation for various operations like filtering, downsampling, upsampling, etc. are developed and used to parameterize synthesis bank solutions. This parameterization is consistent and at least one optimum synthesis filter bank is guaranteed. We design the synthesis bank for a more general set up in comparison to existing methods [57, 58, 59, 60, 61]. The synthesis filter length can be arbitrary but given. Further, we have used delay in the reconstruction as a parameter to have flexibility in the synthesis filter design, and using that, we obtain better reconstruction error than possible with a zero delay design. Another advantage of our analysis is that we obtain a linear time-invariant (LTI) optimal matrix filter unlike the linear periodically time-varying (LPTV) optimal filters obtained in earlier works [49, 62]. Our method always results in a solution in comparison to LMI based methods [55, 56] where a solution for an LMI may not be feasible [67]. The necessary and sufficient condition on a synthesis bank solution to result in PR is derived. Under this condition, the polyphase matrices for the analysis and synthesis banks satisfy the general pseudocirculant property [68, 69], unlike the limited identity matrix condition in the earlier works [23, 29]. The necessary and sufficient conditions on the analysis bank and the delay for a PR synthesis bank solution to exist are derived. This analysis forms an alternative to the frame theory based approach in previous works [29, 30]. Time-domain pseudocirculant result derived in this paper form an important theoretical contribution to the multirate theory. The derivation of a range of delay values for which PR is possible is another contribution. We have also extended the design methodology for any arbitrary NUFB. Under experimental results, a number of experiments are presented to highlight contributions and to show application of our work.

I-B Organization

This paper is organized as follows: Section II discusses the various matrix representations for time-domain signals and operations. In Section III, we introduce the framework used to design the synthesis filter for a given analysis bank. In Section IV, we derive the pseudocirculant condition in time-domain and present the necessary and sufficient conditions on an analysis bank and on the delay to result in PR. We provide the experimental results in Section V and finally, concluding remarks are given in Section VI.

I-C Notation

The $(i,j)$-th entry of a matrix $\mathbf{B}$ is given by $B_{i,j}$ or $B(i,j)$. The transpose and conjugate-transpose of a matrix or vector quantity are denoted by the superscripts $(.)^{T}$ and $(.)^{H}$, respectively. An identity matrix of size $P$ is represented by $\mathbf{I}$ when the size is implicit and by $\mathbf{I}_{P}$ otherwise. The pseudoinverse of a matrix $\mathbf{A}$ is denoted by $\mathbf{A}^{\dagger}$. The superscript $(.)^{*}$ stands for conjugation of a quantity. Convolution of two sequences is denoted by $\ast$. The ceiling function is represented by $\left\lceil.\right\rceil$.

II Matrix Representations

In this section, we develop the matrix representations for various multirate operations, which are required for an effective time-domain analysis of multirate filter banks.

II-A Time-domain Sequence Representation

A vector signal is developed from the given signal $v(n)$ by stacking the past $P$ samples as:

[TABLE]

It is termed as an observation vector for the signal $v(n)$ [70]. Its flipped version or reversal can be obtained as

[TABLE]

where $\mathbf{J}_{P}$ is a counter-identity or reverse-identity matrix. For a matrix $\mathbf{B}$, reversal operation results in reversing of order along both rows and columns. Mathematically, the same can be expressed as

[TABLE]

if the size of the matrix is $M\times N$. For two matrices $\mathbf{A}$ and $\mathbf{B}$ whose product is defined, we have the following relation:

[TABLE]

II-B Decimation

Consider a system, as shown in Fig. 1, with a $M$-fold decimator preceded by a time-advance by $l$. The input and output of it are related as

[TABLE]

Such a system can be used to obtain polyphase components of a signal when $l$ is in the range $[0,M-1]$. Two matrix representations are possible for a system. If an observation vector

[TABLE]

for the output is to be evaluated, then the same can be done as follows:

[TABLE]

where $\mathbf{D}_{P,M(P-1)+1}^{0}$ is known as a decimation matrix. On the other hand, if an input sample vector

[TABLE]

is given, then the corresponding output samples can be obtained as

[TABLE]

where the elements of the decimation matrix $\mathbf{D}_{\lceil P/M\rceil,P}^{l}$ are given by

[TABLE]

Here $l$ is in the range $0\leq l\leq M-1$. For $l\geq M$, the same matrix representations can be extended.

II-C Expansion

If an input observation vector

[TABLE]

is available, then the output of a $M$-fold expander can be obtained as

[TABLE]

where $\mathbf{U}_{M(P-1)+1,P}$ is an expansion matrix. This matrix is related to a decimation matrix as

[TABLE]

It can be easily verified that

[TABLE]

II-D Convolution

Consider a causal linear time-invariant (LTI) filter with impulse response $\{h(0),\allowbreak h(1),\allowbreak\dots,\allowbreak h(Q-1)\}$. Its input and output are related as

[TABLE]

We present the matrix representations for the above relation for two use cases. Under first case, consider an observation vector

[TABLE]

for the output which needs to be evaluated. The same can be achieved using

[TABLE]

where

[TABLE]

and $\mathbf{H}$ is a convolution matrix given by

[TABLE]

The second representation is obtained for a scenario where an input observation vector is given, and the corresponding output samples need to be evaluated. Consider an input observation vector

[TABLE]

then the complete output due to it can be expressed as

[TABLE]

or

[TABLE]

It is implicit in the above representation that $x(n)=0$ for $n\notin[0,P-1]$. The utility of the representation is in evaluation of convolution of two filters. The two convolution matrices, so far obtained, are related with each other as:

[TABLE]

where $\mkern 1.5mu\overline{\mkern-2.0mu\mathbf{H}\mkern-2.0mu}\mkern 1.5mu$ denotes the reversal of the matrix.

In the following section, we address the design of a synthesis filter bank using the time-domain representations discussed so far.

III Design Framework

We have multirate observations, derived from various sub bands, and we aim to reconstruct or estimate the original high resolution signal. Initially, we consider observations to be at the same low rate and model them as output of a uniform analysis bank. The case when the observations are at different rates is addressed later in this paper. The subbands are applied to a matrix synthesis filter, as shown in Fig. 2, where reconstruction happens at a low rate to reduce computational cost. The output of the synthesis filter is then unblocked [71] to get the final high resolution signal. We model our problem of designing the matrix synthesis filter bank in the framework of Wiener filtering. The Wiener filter is an optimum filter to minimize the mean-square error (MSE). For a multi-input multi-output system, the solution is known as a matrix Wiener filter. The details of the design setup are explained next.

Consider the filter bank shown in Fig. 2. A wide-sense stationary (WSS) input signal $u(n)$ passes through a given analysis bank consisting of the filters $\left\{h_{i}(n)\right\}_{i=0}^{L-1}$. The analysis filters are assumed to be FIR and causal. The output signals $\left\{x_{i}(n)\right\}_{i=0}^{L-1}$ of the analysis bank are uniformly decimated by the factor $M$ to obtain the subband signals $\left\{v_{i}(n)\right\}_{i=0}^{L-1}$ which are then applied to a matrix synthesis filter $\mathbf{A}(n)$. We choose the synthesis filter to be FIR and causal, and formulate a Wiener filtering problem for a given length of the synthesis filter. The desired signal for the problem is derived from the input signal $u(n)$ by adding a given delay $d$ and then performing blocking operation [71]. The delay is added to compensate the time lag that appears in the input signal when it passes through the filter bank. The blocking operation ensures that the desired signal $\mathbf{d}(n)$ is jointly stationary with the subband signal $\mathbf{v}(n)$ [72]—a necessary condition required for application of Wiener filtering [73].

III-A FIR Matrix Synthesis Filter Expression

In Appendix VI, we obtain the following expression for an FIR matrix Wiener filter with length $P$:

[TABLE]

where $\mathbf{R}_{vv}$ is the auto-correlation matrix of the input signal $\mathbf{v}(n)$ to the filter, $\mathbf{r}_{d_{i}v}$ is the cross-correlation vector between the desired signal component $d_{i}(n)$ and the signal $\mathbf{v}(n)$, and $\mathbf{a}_{i}$ is the $i$-th row of the filter. This expression can be simplified in the context of the UFB given in Fig. 2. In order to do that, a vector $\mathbf{v}_{\text{s}}(n)$ is created by stacking the observation vectors of the subband signals as

[TABLE]

where $\mathbf{v}_{j}(n)$ is an observation vector representing the $j$-th input to the Wiener filter at time $n$ and is given by

[TABLE]

The quantities $\mathbf{R}_{vv}$ and $\mathbf{r}_{d_{i}v}$ can be expressed in terms of $\mathbf{v}_{\text{s}}(n)$ as

[TABLE]

To make use of the Wiener filter in (24), we derive an expression of the vector $\mathbf{v}_{\text{s}}(n)$ in terms of the input to the filter bank. We can note from Fig. 2 that $v_{j}(n)=x_{j}(Mn)$ and this relationship can also be expressed in matrix form as

[TABLE]

where

[TABLE]

Further, the vector $\mathbf{x}_{j}(Mn)$ can be expressed in terms of the input $u(n)$. We start with an assumption that all the analysis filters are causal and have equal length $Q$. It can be assured by appropriate zero padding. Now we can obtain $\mathbf{x}_{j}(Mn)$ as

[TABLE]

where $\mathbf{H}_{j}$ is a convolution matrix. The entire input to the Wiener filter can be expressed as follows:

[TABLE]

where $\mathbf{K}$ is a matrix of size $LP\times[M(P-1)+Q]$. Using this result, the auto-correlation matrix $\mathbf{R}_{vv}$ can be written as

[TABLE]

Similarly the cross-correlation vector $\mathbf{r}_{d_{i}v}$ becomes

[TABLE]

where $\mathbf{r}_{uu}^{i+d}$ is a row vector given by

[TABLE]

By substituting these results in (24), we obtain

[TABLE]

The resulting equation is of form

[TABLE]

where

[TABLE]

If $\mathbf{A}$ is invertible, then we have a unique solution. Otherwise, the existence of a solution depends on the nature of the system of equations. If the same is consistent, then all solutions can be expressed as

[TABLE]

where $\mathbf{w}$ is an arbitrary vector. In the next theorem, we prove that the equation always results in a solution.

Theorem 1**.**

For an FIR uniform analysis bank, there will be at least one optimal synthesis bank solution that minimizes the mean-square error.

Proof.

We have to prove that the solution given by (38) is always an exact solution of the equation. Substituting the solution in the equation, we obtain

[TABLE]

The matrix $\mathbf{A}$ can be expressed as

[TABLE]

by performing eigendecomposition of the auto-correlation matrix $\mathbf{R}_{uu}$. With this,

[TABLE]

Using the pseudoinverse identities [74], the equation reduces to

[TABLE]

The auto-correlation matrix $\mathbf{R}_{uu}$ is in general positive definite but always positive semi-definite. We consider the solution of equation under these two cases:

Case I: Positive Definite Correlation Matrix

In this case, the auto-correlation matrix is invertible, so is the diagonal matrix $\mathbf{\Lambda}$. Hence we can process the equation as follows:

[TABLE] 2. 2.

Case II: Positive Semi-Definite Correlation Matrix

The diagonal matrix $\mathbf{\Lambda}$ is not invertible under this case. However we exploit the relation between the vector $\mathbf{r}_{uu}^{i+d}$ and the matrix $\mathbf{R}_{uu}$. For $(i+d)<(M(P-1)+Q)$, the vector $\mathbf{r}^{i+d}_{uu}$ is one of the row of the correlation matrix. If $(i+d)\geq(M(P-1)+Q)$, then the vector $\mathbf{r}^{i+d}_{uu}$ can still be expressed as a linear combination of the rows of the auto-correlation matrix [75]. Thus, we can always express the vector $\mathbf{r}^{i+d}_{uu}$ as

[TABLE]

where $\mathbf{R}_{uu}(j,:)$ is the $j$-th row of the correlation matrix. With this,

[TABLE]

We can observe that we have

[TABLE]

under both the cases and hence the solution given by (38) is exact. ∎

From the above theorem, we can conclude that if the matrix $\mathbf{A}$ is non-invertible, we have infinite number of optimum solutions possible. This gives an opportunity to optimize the synthesis bank for a desired property.

III-A1 Computational Complexity

The computational complexity can be evaluated from various algebraic manipulations required for the solution in (38). To start with, the matrix $\mathbf{A}$ evaluation needs O\big{(}LP(M(P-1)+Q)^{2}\big{)}+O\big{(}(LP)^{2}\big{)} calculations and its pseudoinverse requires O\big{(}(LP)^{3}\big{)} calculations. The term $(\mathbf{I}-\mathbf{A}^{\dagger}\mathbf{A})$ can then be evaluated with O\big{(}(LP)^{3}\big{)} cost. As the vector $\mathbf{b}_{i}$ needs to be evaluated for each $i\in[0,M-1]$, the total cost of obtaining $\mathbf{b}_{i}$’s is O\big{(}MLP(M(P-1)+Q)\big{)}. All the products $\mathbf{A}^{\dagger}\mathbf{b}_{i}$ and $(\mathbf{I}-\mathbf{A}^{\dagger}\mathbf{A})\mathbf{w}$ can then be obtained with complexity O\big{(}M(LP)^{2}\big{)}. Therefore, the complexity is quadratic with respect to $M$ and $Q$, and cubic with respect to $L$ and $P$.

Our synthesis filter is multiple-input multiple-output (MIMO) FIR system, thus enabling the operations at the lowest sampling rate. The MIMO implementation brings down the complexity drastically, thereby for each reconstructed output sample, the number of multiplication and additions required are only $LP$ and $(LP-1)$, respectively.

III-B Extension to NUFB

The results developed so far for a UFB can be extended for an NUFB case. Before applying Wiener filtering on an NUFB, we convert it into an equivalent UFB with the help of blocking operation. This conversion is needed as Wiener filtering requires joint stationarity of the input and the desired signals [72]. Each subband of the NUFB with a decimation factor $M_{i}$ is blocked with a factor

[TABLE]

where $L$ is the number of subbands and $M$ is the LCM of all decimation factors. The process is shown in Fig. 3a. An equivalent representation of the process, as shown in Fig. 3b, can be derived using the noble identities [72, 68]. Thus with blocking, the NUFB behaves like a UFB with decimation factor equal to $M$ and with analysis filters given by

[TABLE]

Now the Wiener synthesis filter for the NUFB can be obtained using (38).

III-C Minimum Mean-square Error

We now determine the MSE for a given length of the Wiener filter. The result can be used to choose the length of the synthesis filter based on the desired level of the error. The MSE is defined as

[TABLE]

Its minimum value is achieved using the Wiener filter and is denoted by $J_{\text{min}}$. Further, $J_{\text{min}}$ is the sum of channel-wise MSEs which are given by

[TABLE]

Using the principle of orthogonality [76], the above equation reduces to

[TABLE]

Substituting the expressions for $\mathbf{r}_{d_{i}v}$ and $\mathbf{a}_{i}$ from (33) and (38), respectively, we obtain

[TABLE]

The term with the arbitrary vector is zero due to (46) and the equation reduces to

[TABLE]

This proves that all the solutions for the synthesis bank results in the same MSE, which is intuitive.

III-C1 Optimal Delay

The channel MSE $J_{\text{min}}^{i}$ depends on the delay $d$. The same is also observed in our simulation results. The exact form of the relation can be obtained from (37) as

[TABLE]

$J_{\text{min}}^{i}$ can be observed as a quadratic form with the matrix $\mathbf{B}$ being a Hermitian matrix. If we are interested to obtain an optimal delay value, we can use the existing results [77, 78] on the quadratic forms.

In the following section, we derive the necessary and sufficient conditions on a filter bank for the MSE to be zero.

IV FIR Perfect Reconstruction Filter Bank

In the following analysis, we prove that the designed synthesis bank can achieve PR provided that the analysis bank satisfies certain condition. To justify PR property for any arbitrary input signal on a strong analytical background, we verify the pseudocirculant property [68, 69] for the solution.

IV-A Pseudocirculant Property

We evaluate the product $\mathbf{P}(z)=\mathbf{A}(z)\mathbf{E}(z)$ in time-domain where $\mathbf{A}(z)$ and $\mathbf{E}(z)$ are the polyphase component matrices for the synthesis bank and the analysis bank, respectively [69]. The matrix $\textbf{P}(z)$ is designed as an identity matrix in previous works [23, 29]. Our design approach is more general where $\textbf{P}(z)$ satisfies the pseudocirculant property. In time-domain, we obtain the matrix $\textbf{P}(z)$ from convolution of the two matrix filters as

[TABLE]

In the above equation, the filter $e_{r,j}(n)$ can be obtained from the analysis filters as [68]

[TABLE]

As $Q$ is the length of the analysis filters, the length of the filter $e_{r,j}(n)$ will be at max

[TABLE]

By zero padding, we ensure that all the constituent filters of the polyphase filter $\mathbf{E}(n)$ have the same length $T$. With that, the convolution given in (55) can be expressed in matrix form using (22) as

[TABLE]

where $\tilde{\mathbf{E}}_{r,j}$ is a convolution matrix of size $(P+T-1)\times P$. We like to express this convolution matrix in terms of the analysis filter $h_{r}(n)$. For that, we use a polyphase identity, as shown in Fig. 4, obtained from the existing result [68] by finding the [math]-th polyphase component of a time-advanced filter $h_{r}(n+j)$.

Using this identity, we obtain the following relation

[TABLE]

It can be seen that Fig. 4b contains a combination of upsampling, filtering, time-advance and downsampling. Even then the matrix representations developed in this paper for various multirate operations have been successfully applied to deal with it. Continuing forward, we use the above relation in (58) and obtain

[TABLE]

Next we derive a necessary and sufficient condition, in the form of a theorem, for the analysis and synthesis banks of a filter bank to result in PR.

Theorem 2**.**

An FIR uniform filter bank is PR if and only if

[TABLE]

where $c$ is a non-zero constant and the column vector $\mathbf{e}_{i+d}$ is defined as follows:

[TABLE]

Proof.

The proof follows the Vaidyanathan’s well-known pseudocirculant property. A PR filter bank must satisfies the generalized pseudocirculant condition given by [68]

[TABLE]

We have evaluated $\mathbf{P}(z)$ in time-domain and next we find conditions for it to be pseudocirculant. Let $\mathbf{K}^{T}\mathbf{a}_{i}=\mathbf{s}_{i}$, then (60) becomes

[TABLE]

If the matrix filter $\mathbf{P}(z)$ is pseudocirculant, then the following properties are true in time-domain:

Property 1: In a row, only one $\mathbf{p}_{i,j}$ is non-zero i.e. only one out of $\mathbf{p}_{i,0},\allowbreak\mathbf{p}_{i,1},\allowbreak\dots,\allowbreak\mathbf{p}_{i,M-1}$ is non-zero. 2. 2.

Property 2: In a column, only one $\mathbf{p}_{i,j}$ is non-zero i.e. only one out of $\mathbf{p}_{0,j},\allowbreak\mathbf{p}_{1,j},\allowbreak\dots,\allowbreak\mathbf{p}_{M-1,j}$ is non-zero. 3. 3.

Property 3: For a vector $\mathbf{p}_{i,j}$, at most one entry is non-zero. 4. 4.

Property 4: All the entries on upper-left to lower-right diagonals are the same i.e.

[TABLE] 5. 5.

Property 5: The first column will have

[TABLE]

The above properties will result in the following constraints on the vector $\mathbf{s}_{i}$:

Property 1 implies that only one vector

[TABLE]

is non-zero in a row. Further, Property 3 means that only one element of the non-zero vector is non-zero. Together these two conditions imply that only one element of the vector $\mathbf{s}_{i}$ is non-zero. Assume that $\mathbf{p}_{i,j}\neq 0$ for $j=j_{i}$ and ${p}_{i,j_{i}}(l_{i})=s_{i}(j_{i}+Ml_{i})\neq 0$, then we can express

[TABLE] 2. 2.

Property 4 results in

[TABLE]

From the above equation, all the elements of the vector $\mathbf{s}_{i+1}$, except at indices that are multiples of $M$, can be obtained from the vector $\mathbf{s}_{i}$. 3. 3.

Property 5 implies

[TABLE]

Another implication is that

[TABLE]

Combining the two results (69) and (71), we get

[TABLE]

Thus the vector $\mathbf{s}_{i+1}$ is just a shifted version of the vector $\mathbf{s}_{i}$.

The last constraint implies that all the vectors $\mathbf{s}_{i}$ can be obtained from a single vector $\mathbf{s}_{0}$. If the vector $\mathbf{s}_{0}$ has form

[TABLE]

then all the vectors $\mathbf{s}_{i}$ can be expressed as

[TABLE]

Due to the above equation, we have

[TABLE]

In z-domain, this means that

[TABLE]

and the reconstructed output is $\hat{u}(n)=cu(n-d_{0}-M+1)$. One can easily see that the variable $d_{0}$ is the same as the delay $d$ in the desired signal path. ∎

It is to be noted that the above theorem puts restriction on both the analysis and synthesis banks. It may happen that a solution may not exist for the synthesis bank to result in PR. We will first obtain a parameterization of all PR solutions for the synthesis bank and later derive necessary & sufficient conditions for the same to exist.

IV-B Parametric Form for PR Solutions

In the previous subsection, we obtained the following condition on the analysis and synthesis banks to have PR:

[TABLE]

The above equation is of form

[TABLE]

and thus represents a system of equations. If the same is consistent, then all the PR solutions for the synthesis bank are given by

[TABLE]

where $\mathbf{w}$ is an arbitrary vector. We can observe that the solutions have turned out to be independent of the input signal.

IV-C Necessary and Sufficient Conditions on Analysis Bank and Delay

In order for a PR solution to exist, the analysis bank has to satisfy certain condition. We present the same in the form of a theorem about the matrix $\mathbf{K}$ which represents the analysis bank.

Theorem 3**.**

For an FIR uniform analysis bank, at least one PR synthesis bank solution is possible if the matrix $\mathbf{K}$ satisfies the relation

[TABLE]

where $\mathbf{0}_{r\times r}$ is the largest zero sub-block possible with $r\geq M$, and $\mathbf{B}_{p\times p}$, $\mathbf{D}_{(q-p-r)\times(q-p-r)}$ and $\mathbf{C}_{p\times(q-p-r)}$ are arbitrary matrices with the first two being Hermitian.

Proof.

We require the system of linear equations given by

[TABLE]

to be consistent to have an exact solution. In order to obtain the condition for that, we substitute a PR solution in the above equation and obtain

[TABLE]

Using the fact that $\mathbf{K}\mathbf{K}^{\dagger}\mathbf{K}=\mathbf{K}$, the equation simplifies to

[TABLE]

If the matrix $\mathbf{K}^{\dagger}\mathbf{K}$ is of dimension $q\times q$, then in order for the above equation to be satisfied, we require

[TABLE]

where $p\leq d$, $r\geq M$, and $\mathbf{X}_{q\times p}$ and $\mathbf{Y}_{q\times(q-p-r)}$ are arbitrary matrices of dimensions $q\times p$ and $q\times(q-p-r)$, respectively. As $\mathbf{K}^{\dagger}\mathbf{K}$ is a Hermitian matrix, we have

[TABLE]

Thus the two matrices $\mathbf{X}_{q\times p}$ and $\mathbf{Y}_{q\times(q-p-r)}$ have following form:

[TABLE]

and

[TABLE]

where matrices $\mathbf{B}_{p\times p}$ and $\mathbf{D}_{(q-p-r)\times(q-p-r)}$ are Hermitian matrices. Hence the condition on the matrix $\mathbf{K}$ becomes

[TABLE]

∎

The PR requirement imposes the following condition on the desired signal delay $d$.

Lemma 1**.**

With an analysis bank matrix $\mathbf{K}$ satisfying the above condition, PR is possible only when the delay $d$ is in the range:

[TABLE]

Proof.

The proof easily follows from the equations (83) and (84). ∎

The above lemma together with Theorem 3 constitute the necessary and sufficient conditions for a PR synthesis solution to exist. The lemma points that with suitable delay, PR is achievable for an analysis bank which is otherwise not possible. This way we have attempted to expand the space of PR FIR filter banks. Further, we can have PR for a number of delay values which gives an opportunity to use delay as an optimization parameter, based on various design criterion like Hermitian symmetry, time or frequency localization [23], coding gain [79], etc. For this, the number of channels required is $L>M$ and the channels should satisfy the PR condition given by Theorem 3. For a fixed delay, it is an unconstrained problem with respect to the arbitrary vector $\mathbf{w}$ [23, 30]. The required optimized solution can be selected among the solutions for various delays by imposing an additional constraint. Design specific constraints and their solutions can be considered for future work.

It can be seen that the matrix $\mathbf{K}$ easily satisfy the PR condition if it has full column rank. The PR condition induces certain restriction on its null space which is given by the following lemma.

Lemma 2**.**

For a PR UFB, the analysis bank matrix $\mathbf{K}$ has a null space spanned by the following set of vectors:

[TABLE]

where $\mathbf{f}_{1\times p}$ and $\mathbf{g}_{1\times(q-p-r)}$ are arbitrary vectors.

Proof.

The matrix $\mathbf{K}$ represents the analysis bank of a filter bank. Hence, when an input signal $u(n)$ is applied to the analysis bank, the resulting subband vector can be obtained as

[TABLE]

Consider a case where the vector $\bar{\mathbf{u}}(Mn)$ lies in the null space of the matrix $\mathbf{K}$. As a result, the vector $\mathbf{v}_{s}(n)$ as well as the output of the synthesis bank, both are zero. The reconstructed signal is then

[TABLE]

For PR, we require

[TABLE]

Combining the two, we get

[TABLE]

Thus the vector $\bar{\mathbf{u}}(Mn)$ should be of form

[TABLE]

A basis vector for the null space of the matrix $\mathbf{K}$ should satisfy the above condition for all possible values of $d$ and hence should have the form

[TABLE]

where $\mathbf{f}_{1\times p}$ and $\mathbf{g}_{1\times(q-p-r)}$ are arbitrary vectors. ∎

In this section, we presented various necessary and sufficient conditions for a PR synthesis bank solution to exist. In the next section, we present some experiments which validates our results.

V Experimental Results

In this section, we verify the derived results with a number of experiments. These experiments are chosen to demonstrate the wide scope of application of our design. In Experiment 1, we obtain the MSE for various delay values in the desired signal and compare the results. We present the improvement in reconstruction, in terms of the MSE, with addition of a new subband in Experiment 2. We also examine the impact of delay on PR. In Experiment 3, a PR non-uniform filter bank obtained through our design is verified with a non-stationary input signal.

V-A Experiment 1

We consider a 2-channel maximally decimated uniform analysis bank designed using fir1 routine of Matlab with the following parameters: the analysis filter $h_{0}(n)$ is low-pass filter with cut-off 0.6 and length 9, and $h_{1}(n)$ is high-pass filter with cut-off 0.4 and length 10. The input to the filter bank is a zero mean, unit variance second-order autoregressive(AR(2)) process with coefficients $0.7$ and $0.1$. We choose the length of the constituent filters of the matrix synthesis filter to be 11 and design the filter for various delay values using (38). We obtain ensembled-averages of the channel-wise squared errors, shown in Figs. 5a and  5b, for these synthesis filters by performing 200 independent runs. The corresponding theoretical values, obtained using (53), are provided in Table I. It is to be mentioned here that the analysis bank does not satisfy Theorem 3. We can see that the observed error values are closer to the theoretical values for both channels. Further, we can notice that the best result among the selected delay values, in terms of the total MSE, is with delay 10. This shows that a better reconstruction can be obtained if some delay in the reconstructed output is acceptable.

V-B Experiment 2

In this experiment, we examine the effect of number of subbands on the MSE along with the role of delay in achieving PR. We consider an Extended Lapped Transform (ELT) filter bank with decimation factor $M=4$ and overlapping factor $K^{\prime}=2$. The prototype filter for ELT is designed using the optimized butterflies angles mentioned in [21]. The input signal to the filter bank is a realization of an AR(1) process with coefficient $0.95$, zero mean and unit variance. The length of the synthesis filter is $4$ and the delay is $10$. We vary the number of subbands ($L$) available to the synthesis filter and obtain the MSE for each channel. On analyzing the obtained results, as given in Table II, we can conclude that for number of subbands $L<4$, there is a variation in the MSE with subbands chosen. However independent of the choice of subbands, if we add additional subband to the system, then we always see an improvement in the total MSE as well as in the channel-wise MSEs. This is possible since we have exploited the correlation between the input signal components while designing the matrix synthesis filter. With the number of subbands $L=4$, the analysis bank satisfies the PR condition given by Theorem 3. We vary the delay and find out that PR is possible only with $d=12$. This is consistent with the range of delay values provided by (89).

V-C Experiment 3

In this experiment, we design a PR synthesis bank for a non-uniform analysis bank and test it with a non-stationary input signal. We have the following setup:

[TABLE]

The analysis filters are designed using Parks-McClellan algorithm. We convert the NUFB into an equivalent UFB with the help of blocking operation. The resulting UFB has decimation factor equal to 6 and the following analysis filters:

[TABLE]

We choose the length of the synthesis filter to be 7 and delay $d=0$. The resulting matrix $\mathbf{K}$ satisfies the PR condition (80) and using (79), we obtain a PR solution with $\mathbf{w}=\mathbf{0}$. We apply a non-stationary signal generated using $\texttt{randn}(n)\sin(0.1n^{2})$ to the filter bank and evaluate the reconstruction error as (refer Fig. 2):

[TABLE]

Fig. 6 shows the resulting error versus time curve. The small error which we obtain is due to finite precision in MATLAB.

VI Conclusion

We have addressed an important problem of designing an optimal FIR synthesis bank for a given non-maximally decimated analysis bank. To design such a synthesis bank, a framework was developed to apply matrix Wiener filtering to the problem and solution is obtained for a given length of synthesis filter and delay in the reconstructed output. We experimentally showed that better MSE can be obtained in some cases by adding appropriate delay in the reconstructed output. We found that the reconstruction error keeps decreasing as more and more multirate observations are added for a signal. The matrix representation developed in this paper can be utilized to carry out time-domain analysis of any multirate circuit. The infinite number of PR or non-PR synthesis bank solutions possible with our design can be optimized for a desired property in a future work.

[FIR Matrix Wiener Filter Derivation]

Consider an optimum matrix filtering problem as shown in Fig. 7 where the vector signal $\mathbf{v}(n)$ is the input and the signal $\mathbf{d}(n)$ is an estimate of the desired signal. In order to obtain the Wiener solution, the input and desired signals are assumed to be jointly stationary [73]. The Wiener solution can be derived for either FIR or IIR case. Here we derive the FIR Wiener solution for a matrix filter whose constituent filters are causal and of length $P$.

We begin by applying the principle of orthogonality [76] to the setup and obtain

[TABLE]

The above equation can be written as

[TABLE]

As the input and desired signals are jointly stationary, we can write

[TABLE]

The above equation can be expressed in matrix form as

[TABLE]

where

[TABLE]

A further compact form is

[TABLE]

where

[TABLE]

The obtained result is a system of equations containing $LP$ equations in terms of $LP$ unknowns. These are the Wiener-Hopf equations for the matrix Wiener filter. Using them, we can determine the $i$-th row of the matrix Wiener filter.