A Steady-State Kalman Predictor-Based Filtering Strategy for Non-Overlapping Sub-Band Spectral Estimation

Abstract

This paper focuses on suppressing spectral overlap for sub-band spectral estimation, with which we can greatly decrease the computational complexity of existing spectral estimation algorithms, such as nonlinear least squares spectral analysis and non-quadratic regularized sparse representation. Firstly, our study shows that the nominal ability of the high-order analysis filter to suppress spectral overlap is greatly weakened when filtering a finite-length sequence, because many meaningless zeros are used as samples in convolution operations. Next, an extrapolation-based filtering strategy is proposed to produce a series of estimates as the substitutions of the zeros and to recover the suppression ability. Meanwhile, a steady-state Kalman predictor is applied to perform a linearly-optimal extrapolation. Finally, several typical methods for spectral analysis are applied to demonstrate the effectiveness of the proposed strategy.

1. Introduction

As one of the most important tools, spectral estimation [1] has been extensively applied in radar, sonar and control systems, in the economics, meteorology and astronomy fields, speech, audio, seismic and biomedical signal processing, and so on. In particular, sparse representation [2–4] opens an exciting new vision for spectral analysis. However, such methods are usually accompanied by high computational complexity, which makes their availability somewhat limited.

Sub-band decomposition-based spectral estimation (SDSE) [5] is an important research direction in spectral estimation, because it has several advantageous features, e.g., computational complexity decrease, model order reduction, spectral density whiteness, reduction of linear prediction error for autoregressive (AR) estimation and the increment of both frequency spacing and local signal-to-noise ratio (SNR) [6]. These features have been theoretically demonstrated under the hypothesis of the ideal infinitely-sharp bandpass filter bank [7]. Subsequent studies [8–10] indicate that these benefits aid complex frequency estimation in sub-bands, thereby enabling better estimation performance than that achieved in full-band. In addition, the computational complexity of most algorithms for spectral analysis has a superlinear relationship with the data size, and sub-band decomposition can considerably speed up these algorithms. Independently handling each sub-band enables parallel processing, which can further improve the computational efficiency. Both advantages are crucial for reducing the computational burden, especially when analyzing multi-dimensional big data, such as polarimetric and/or interferometric synthetic aperture radar images of large scenes.

Unfortunately, the ideal infinitely-sharp bandpass filter cannot be physically realized, and non-ideal (realizable) filters introduce energy leakage and/or frequency aliasing phenomena [11]. Due to these non-ideal frequency characteristics of analysis filters, spectral overlap between any two contiguous sub-bands occurs during the sub-band decomposition. Then, the performance of SDSE severely degrades.

In the relevant literature, several methods have been proposed to mitigate spectral overlap. We classify these methods into three categories. The first category is defined as ideal frequency domain filtering with a strict box-like spectrum, such as “ideal” Hilbert transform-based half-band filters [9] and harmonic wavelet transform-based filters [12,13]. Theoretically, sub-band decomposition with these filters is immune to spectral overlap. However, discrete Fourier transform will inevitably induce spectral energy leakage, which can likewise distort sub-band decomposition. The second category is known as convolution filtering with wavelet packet filters [8], Kaiser window-based prototype cosine modulated filters, discrete cosine transform (DCT) IV filters [10] and Comb filters [6,14]. It seems that increasing the filter order can improve the filtering performance and also the spectral overlap suppression capability. However, in the context of involving a finite-length sequence and performing convolution filtering, the nominal improvement of performance will lead to spectral energy leakage and inferior filtering accuracy [10]. Considering the compromise between suppressing spectral overlap and reducing spectral energy leakage, we have to restrict the filter order. The third category is frequency-selective filtering, and a representative method is SELF-SVD (singular value decomposition-based method in a selected frequency band) [15]. Essentially, SELF-SVD attempts to attenuate the interferences of the out-of-band components by the post-multiplication with an orthogonal projection matrix. Unfortunately, the attenuation is often insufficient when the out-of-band components are much stronger than the in-band components or the SNR is relatively low. In this case, the estimation of the in-band frequencies is seriously affected.

In this paper, a new filtering strategy is proposed to suppress spectral overlap for sub-band spectral estimation. First, we discuss the formation mechanism of spectral overlap. Nominally, a high-order finite impulse response (FIR) filter usually has a powerful ability in spectral overlap suppression. However, once we perform such a filter on a finite-length sequence with the convolution operation, the non-given samples at the forward and backward sampling periods of the sequence are assumed to be zeros. A certain filtering error therefore occurs and conversely disrupts the decomposed sub-bands. As a result, sub-band spectral analysis severely suffers from the mutual overlap of adjacent sub-band spectra. Second, we propose a filtering strategy to eliminate the filtering error and recover the suppression ability. This strategy intuitively takes the place of the artificial zeros with some extrapolated samples. Toward the problem of data extrapolation, many algorithms have been proposed based on various theories, such as linear prediction [16], Gerchberg–Papoulis [17], Slepian series [18], linear canonical transform [19] and sparse representation [20]. To establish an efficient method for the extrapolation in context and to evaluate the effectiveness of the proposed strategy, we preliminarily develop a linearly-optimal extrapolation based on the classical AR model identification and the Kalman prediction [21–23]. Third, we derive the formulas to estimate the residual filtering error and adapt two common information criteria with adaptive penalty terms for AR order determination. Moreover, equiripple FIR filters are applied as analysis filters in coordination with the proposed filtering, because of their advantageous features [11]. Finally, the entire algorithm and the computational complexity are summarized. Some details, such as the sub-band spectrum mosaicking procedure and parameter selection, are discussed in practice.

The remainder of the paper is organized as follows. In Section 2, the formation mechanism of spectral overlap is discussed. Based on this, a steady-state Kalman predictor-based filtering strategy is developed to suppress the overlapped spectra. In Section 3, the proposed filtering strategy is discussed for SDSE. In Section 4, experimental results with several typical algorithms for spectral analysis demonstrate the effectiveness of the proposed strategy. Finally, Section 5 concludes this paper.

2. Signal Filtering Based on AR Model Identification and Kalman Prediction

This section focuses on signal filtering. To reduce the filtering error induced by convolution filtering, we propose an extrapolation-based filtering strategy and apply a steady-state Kalman predictor for extrapolation. Two criteria with adaptive penalty terms for order determination are developed based on the estimation of the residual filtering error.

2.1. Problem Statement of Signal Filtering

FIR filters are typical linear time-invariant (LTI) systems. According to the linear system theory, the filter can be mathematically expressed as the convolution of its impulse response with the input. Suppose that {x_n} is an input sequence and {h_n} is the impulse response of a causal FIR filter; the filtered sequence {y_n} can be derived as [11]:

$$y_n=h_n*x_n=\sum _{k=0}^{N_f-1}{h}_k{x}_{n-k}$$

(1)

where * denotes the convolution operator and N_f is the filter length (i.e., the length of the impulse response; the relationship between N_f and the filter order N_o can be written as N_f = N_o + 1). Alternatively, taking discrete-time Fourier transforms (DTFT), we can represent Equation (1) in the frequency domain as:

$$Y\left(e^{j\omega }\right)=H\left(e^{j\omega }\right)X\left(e^{j\omega }\right)$$

(2)

In addition, the filtered sequence length L, the input sequence length N and the filter length N_f satisfy:

$$L=N+N_f-1=N+N_o$$

(3)

Theoretically, given a large enough stop-band attenuation, spectral overlap can be thoroughly suppressed. Moreover, the spectral estimation error in sub-bands can be neglected, as long as the width of the transition band and the ripple of the passband are sufficiently small. Nonetheless, the pursuit of excellent filtering performance substantially increases both the filter order and the length of the filtered sequence (refer to Equation (3)). Such a high order is more likely to create error in part or even all of the filtered samples. This result is contrary to our original objective, and the resultant filtering quality is undesirable.

From the perspective of a discrete-time system, the output sequence of the convolution operation is equivalent to the zero-state response of the filter system, because the initial state of every delay cell is zero prior to the excitation of the input sequence. We take the example of the direct-type FIR system [24]. The value of the output sample at any time depends on all or part of the input samples and the system state at that time. The first N_f − 1 output samples suffer from biases, because a part of the delay cells do not yet become input-driven states; analogously, the last N_f − 1 output samples are invalid, because a part of the delay cells restore the initial zero-states. Thus, the length of the valid part of the output sequence, defined as L_υ, satisfies:

$$L_{\upsilon }=L-2\left(N_f-1\right)=N-N_f+1=N-N_o$$

(4)

Actually, if we rewrite Equation (1) in the following matrix form:

$$\underset{\mathbf{\text{y}}}{\underbrace{\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_{L-1}\end{array}\right]}}=\underset{\mathbf{\text{X}}}{\underbrace{\left[\begin{array}{cccc}{x}_0& 0& \cdots & 0\\ x_1& x_0& \ddots & 0\\ ⋮& ⋮& & ⋮\\ x_{N-1}& x_{N-2}& \ddots & x_0\\ 0& x_{N-1}& & x_1\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & x_{N-1}\end{array}\right]}}\underset{\mathbf{\text{h}}}{\underbrace{\left[\begin{array}{c}{h}_0\\ h_1\\ ⋮\\ h_{N_o}\end{array}\right]}}$$

(5)

then we can find that the matrix X possesses many zero elements, which probably makes the outputs y₀, y₁, …, y_No−1; y_L−No, y_L−No+1, _…, y_L−1 not ideal. For example, y₀ = x₀h₀, while the ideal output should be ỹ₀ = x₀h₀ + x₋₁h₁ + x₋₂h₂ + ⋯ + x_−Noh_No. This means that the unknown samples x_−No, x_−No+1, _…, x₋₁ are assumed to be zeros. The filtering error of y₀ is y₀ − ỹ₀ = −(x₋₁h₁ + x₋₂h₂ + ⋯ + x_−Noh_No). Likewise, the outputs y₁, y₂, …, y_No−1; y_L−No, y_L−No+1, _…, y_L−1 all suffer from errors under the zero assumption. Thus, we can conclude that the meaningless zeros are the error sources of the filtering.

Referring to Equation (4), we note that, if the filter order is not less than the length of the input sequence, the output samples are all invalid. Thus, improving filtering performance by means of unlimited increasing filter order is meaningless.

In the next subsection, we will identify an efficient way to resolve this problem.

2.2. Filtering Procedure Based on Signal Extrapolation

The desired output of the filtering process should have two characteristics:

• The original and filtered sequences should be of equal length;

• During the filtering process, the states of the delay cells in the filter system should always maintain input-driven states, i.e., there are no artificial zeros, but authentic samples in X.

As shown in Equation (5), the convolution filtering assumes the unknown samples x_−No, x_−No+1, _…, x₋₁; x_N, x_N+1, _…, x_N+No−1 to be zeros, which leads to the filtering error. Thus, an intuitive thought is to extrapolate the sequence ${\left\{x_n\right\}}_{\left(n=0\right)}^{\left(N-1\right)}$ along two sides to provide a series of estimates for the unknown samples. Taking place of the zeros in the matrix X with these estimates can mitigate the filtering error. The input sequence is extrapolated along both sides, yielding two extrapolated sequences, called Part A and Part B (see Figure 1). Suppose that L_A and L_B are the lengths of Part A and Part B, respectively; then, those L_A + L_B extrapolated samples are used to replace zeros in X. According to Equation (3), the length of the associated output sequence is L_A + L_B + N + N_o. From Equation (4), the effective length of the output can be given by L_A + L_B + N − N_o. To satisfy the requirement that the original and filtered sequence are of equal length, the extrapolated length can be derived as:

$$L_A+L_B+N-N_o=N\Rightarrow L_A+L_B=N_o$$

(6)

Now, we can conclude that the extrapolated length should be equal to the filter order. We define L_G as the constant group delay of the filter. Between time N_o and time N + N_o, the output samples are valid. The output sample at time N_o + n (n = 1, 2, ⋯ , N) corresponds to the input sample at time N_o − L_G + n (n = 1, 2, ⋯ , N), because of the group delay. Consequently, the input sample before time N_o − L_G is merely used as a training sequence of the system state. Thus, we can obtain the relationships:

$$\{\begin{array}{l}{L}_A=N_o-L_G\hfill \\ L_B=L_G\hfill \end{array}$$

(7)

Let x̂_n and ŷ_n be the extrapolated sequence and associated filtered result, respectively. Then, they satisfy:

$$\{\begin{array}{cc}\widehat{x_n}:& L_G-N_o⩽n⩽N+L_G-1\\ h_n:& 0⩽n⩽N_o\\ \widehat{y_n}:& L_G⩽n⩽L_G+N-1\end{array}$$

(8)

$${\widehat{x}}_n=x_n\left(0⩽n⩽N-1\right)$$

(9)

The filtering process can be rewritten in matrix form as:

$$\mathbf{ŷ}=\widehat{\mathbf{\text{X}}}\mathbf{\text{h}}$$

(10)

where:

$$\mathbf{ŷ}={\left[{ŷ}_{L_G},{ŷ}_{L_G+1},\dots ,{ŷ}_{L_G+N-1}\right]}_{\left(N\times 1\right)}^{\text{T}}$$

(11)

and:

$$\widehat{\mathbf{\text{X}}}={\left[\begin{array}{cccc}{\widehat{x}}_{L_G}& {\widehat{x}}_{L_G-1}& \dots & {\widehat{x}}_{L_G-N_o}\\ {\widehat{x}}_{L_G+1}& {\widehat{x}}_{L_G}& \dots & {\widehat{x}}_{L_G-N_o+1}\\ ⋮& ⋮& \ddots & ⋮\\ {\widehat{x}}_{L_G+N-1}& {\widehat{x}}_{L_G+N-2}& \dots & {\widehat{x}}_{L_G+N-N_o-1}\end{array}\right]}_{\left(N\times N_f\right)}$$

(12)

2.3. Signal Extrapolation Based on AR Identification and Kalman Prediction

According to the linear prediction theory [25], the AR model is an all-pole model, whose output variable only linearly depends on its own previous values, that is,

$$\{\begin{array}{l}\Phi \left(q^{-1}\right)x_n={\mathit{\varepsilon }}_n\hfill \\ \Phi \left(q^{-1}\right)=\text{∑}_{l=0}^p{\varphi }_l{q}^{-l}\hfill \end{array},n=0,1,\cdots ,N-1$$

(13)

where q⁻¹ denotes the unit delay, p is the model order, ϕ₀, ϕ₁, …, ϕ_p denote the coefficients of the model and ϕ₀ = 1. The sequence ${\left\{{\mathit{\varepsilon }}_n\right\}}_{n=-\infty }^{\infty }$ is a white noise process, which satisfies:

$$\{\begin{array}{ll}E\left({\mathit{\varepsilon }}_n\right)=0,\hfill & \forall n\hfill \\ E\left({\mathit{\varepsilon }}_n^2\right)={\sigma }^2,\hfill & \forall n\hfill \\ E\left({\mathit{\varepsilon }}_n{\mathit{\varepsilon }}_{n^{\prime }}\right)=0,\hfill & n\ne n^{\prime }\hfill \end{array}$$

(14)

where E (·) denotes the expectation operator.

We choose the forward-backward approach [1] as the coefficient estimator for AR model, for its precision and robustness. Both criteria, including the Akaike information criterion (AIC) and Bayesian information criterion (BIC) [26] can be applied to determine the model order; whereas both criteria sometimes suffer from overfitting. An alternative method of order determination will be discussed in Section 2.4.

A linearly-optimal prediction for AR sequences is derived in [21–23] under the minimum mean square error (MMSE) criterion. However, the prediction formula involves a polynomial long division and a coefficient polynomial recursion [23], making the calculation of the prediction somewhat inconvenient. Alternatively, the following steady-state Kalman predictor [27] provides an equivalent prediction with the MMSE predictor, while offering a simpler formula to facilitate the computation.

The AR model is regarded as a dynamic system. A specific state-space representation for a univariate AR(p) process can be written as [25]:

$$\{\begin{array}{c}{\xi }_{n+1}=\mathbf{\text{F}}{\xi }_n+\mathbf{\Gamma }{\mathit{\varepsilon }}_n\\ x_n=\mathbf{\text{H}}{\xi }_n+{\mathit{\varepsilon }}_n\end{array}$$

(15)

where:

$$\mathbf{\text{F}}={\left[\begin{array}{ccccc}-{\varphi }_1& 1& 0& \cdots & 0\\ -{\varphi }_2& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \cdots & 0\\ -{\varphi }_{p-1}& 0& 0& \cdots & 1\\ -{\varphi }_p& 0& 0& \cdots & 0\end{array}\right]}_{\left(p\times p\right)}$$

(16)

$$\mathbf{\Gamma }={\left[-{\varphi }_1-{\varphi }_2\cdots -{\varphi }_{p-1}-{\varphi }_p\right]}_{\left(p\times 1\right)}^{\text{T}}$$

(17)

and:

$$\mathbf{\text{H}}={\left[10\cdots 0\right]}_{\left(1\times p\right)}$$

(18)

The coefficient polynomials of x_n and ε_n are Φ(q⁻¹) and one, respectively. Since they are relative prime polynomials (or coprime), i.e., the transfer function is irreducible, the system of the AR model is a joint controllable and observable discrete linear stochastic system [28]. Thus, there exists a steady-state Kalman predictor:

$$\{\begin{array}{c}{\widehat{\xi }}_{n+1|n}=F{\widehat{\xi }}_{n|n-1}+\mathbf{\text{K}}{e}_n\\ x_n=\mathbf{\text{H}}{\widehat{\xi }}_{n|n-1}+e_n\end{array}$$

(19)

Since both ε_n and e_n are the innovation processes of x_n, they are equal [27]:

$$e_n={\mathit{\varepsilon }}_n$$

(20)

By comparing Equation (15) with Equation (19), we have:

$$\{\begin{array}{ll}{\xi }_n\hfill & ={\xi }_{n|n-1}\hfill \\ \mathbf{\text{K}}\hfill & =\mathbf{\Gamma }\hfill \end{array}$$

(21)

Therefore, the one-step steady-state Kalman predictor can be derived as [28]:

$$\begin{array}{ll}{\widehat{x}}_{n+1|n}\hfill & =H{\widehat{\xi }}_{n+1|n}=x_{n+1}-{\mathit{\varepsilon }}_{n+1}\hfill \\ \hfill & =x_{n+1}-\Phi \left(q^{-1}\right)x_{n+1}=-\text{∑}_{l=1}^p{\varphi }_l{x}_{n-l+1}\hfill \end{array}$$

(22)

Similarly, the k-step steady-state Kalman predictor can be presented as:

$$\{\begin{array}{l}{\widehat{\xi }}_{n+k|n}=F{\widehat{\xi }}_{n+k-1|n-1}+F^{k-1}\mathbf{\Gamma }{\mathit{\varepsilon }}_n\hfill \\ {\widehat{x}}_{n+k|n}=\mathbf{\text{H}}{\widehat{\xi }}_{n+k|n}\hfill \end{array}$$

(23)

Here, define:

$$\{\begin{array}{c}{\mathbf{\text{F}}}^{k-1}\mathbf{\Gamma }={\left[g_{k_0},g_{k1},\cdots ,g_{k_{p-1}}\right]}^{\mathbf{\text{T}}}\\ G_k\left(q^{-1}\right)=g_{k_0}+g_{k_1}{q}^{-1}+\cdots +g_{k_{p-1}}{q}^{-\left(p-1\right)}\end{array}$$

(24)

Then, we can obtain:

$$\Phi \left(q^{-1}\right){\widehat{x}}_{n+k|n}=G_k\left(q^{-1}\right){\mathit{\varepsilon }}_n$$

(25)

and that is,

$${\widehat{x}}_{n+k|n}=G_k\left(q^{-1}\right)x_n=\text{∑}_{l=0}^{p-1}{g}_{k_l}{x}_{n-l}$$

(26)

Analogously, the k-step backward extrapolation formula can be given by:

$${\widehat{x}}_{n-k|n}=\sum _{l=0}^{p-1}{g}_{k_l}^{*}{x}_{n+l}$$

(27)

where the superscript “*” denotes the complex conjugate operator. To guarantee reasonable and effective extrapolations, the step-size k should satisfy:

$$\{\begin{array}{l}{L}_G-N_o+k⩾0\hfill \\ N+L_G-1-k⩽N-1\hfill \end{array}\Rightarrow L_G⩽k⩽N_o-L_G$$

(28)

In order to evaluate the residual filtering error of the proposed filtering strategy, we derive the mean square error (MSE) in Appendix A1.

2.4. Adaptive Information Criteria for AR Order Determination

Given the impulse response of an analysis filter and AR coefficients, we can directly calculate MSE by Equations (A2) and (A10). The precision of AR coefficient estimation is concerned with AR order. Consequently, the filtering error at different AR orders can be evaluated with the preceding formulas; conversely, the calculation of MSE can be used for order determination.

AIC and BIC are two common information criteria, whose purpose is to find a model with sufficient goodness of fit and a minimum number of free parameters. In terms of the maximum likelihood estimate ${\widehat{\sigma }}_p^2$ we can denote AIC and BIC as [26]:

$$\underset{P}{min}\text{AIC}=log\left({\widehat{\sigma }}_p^2\right)+\frac{2\left(p+1\right)}{N}$$

(29)

$$\underset{P}{min}\text{BIC}=log\left({\widehat{\sigma }}_p^2\right)+\frac{plogN}{N}$$

(30)

As explained in [29], due to the lack of samples, both criteria encounter the risk of overfitting, where the selected order will be larger than the truth order. In particular, AIC has the nonzero overfitting probability as the sample number tends to infinity. Theoretically, both criteria consist of two terms: the first term involves MSE, and it decreases with the increment of the order p; the other term is a penalty that is an increasing function of p. The preferred model order is the one with the lowest AIC or BIC value. As shown in Figure 2a, the objective function curve $\widetilde{{\text{S}}_1{\text{P}}_1{\text{E}}_1}$ reaches its minimum value at the point P₁, which gives the correct order p. However, sometimes, both criteria may fail to determine available orders, and those failures are often related to inadequate penalties. Figure 2b illustrates a representative case. Since the change of the objective function instantly slows down as the order exceeds p, the point P₂ is the preferred point for order determination. However, the penalty strength is insufficient, so that the objective function is still falling after P₂. To handle this situation, we propose an adaptive mechanism to adaptively adjust the penalty strength. A geometric interpretation is depicted in Figure 2b. We assume that the order interval for computation consists of the correct order. Then, the ray $\overrightarrow{{\text{S}}_2{\text{E}}_2}$ forms the X₂ axis, while the ray $\overrightarrow{{\text{O}}_2{\text{Y}}_2}$ forms the Y₂ axis perpendicular to the ray $\overrightarrow{{\text{S}}_2{\text{E}}_2}$ throughout the intersection O₂ of the ray $\overrightarrow{{\text{S}}_2{\text{E}}_2}$ and the objective function axis. Under the new coordination system X₂O₂Y₂, the minimum point P₂ of the curve $\widetilde{{\text{S}}_2{\text{P}}_2{\text{E}}_2}$ can help to determine the correct order. Meanwhile, this modification has no impact on the case that the criterion works well (see Figure 2a).

Therefore the adaptive AIC (AAIC) based on MSE of the residual filtering error can be derived as:

$$\{\begin{array}{l}\underset{p}{min}\text{AAIC}=log\left({\widehat{\sigma }}_p^2\right)+\frac{2}{N}\left(p^{\alpha }+1\right)\hfill \\ s.t.\text{AAIC}\left(p_s\right)=\text{AAIC}\left(p_e\right)\hfill \end{array}$$

(31)

where p_s and p_e denote the start point and the end point of the computing order interval, respectively. If p_s= 1 the adaptive parameter α can be given by

$$\alpha ={log}_{p_e}\left[\frac{N}{2}log\left(\frac{{\widehat{\sigma }}_1^2}{{\widehat{\sigma }}_{p_e}^2}\right)+1\right]$$

(32)

Analogously, the adaptive BIC (ABIC) can be represented as:

$$\{\begin{array}{l}\underset{p}{min}\text{ABIC}=log\left({\widehat{\sigma }}_p^2\right)+p^{\beta \frac{logN}{N}}\hfill \\ \beta ={log}_{p_e}\left[\frac{N}{logN}log\left(\frac{{\widehat{\sigma }}_1^2}{{\widehat{\sigma }}_{p_e}^2}\right)+1\right]\hfill \end{array}$$

(33)

3. Implementation of SDSE

In this section, we discuss the implementation details of SDSE based on the proposed filtering strategy. In particular, equiripple FIR filters are used as the analysis filters for their advantageous features. To suppress spectral overlap and improve spectral precision in practice, we introduce a mosaicking operation for sub-band spectra and discuss the compensation of the residual error of the composite spectrum. After that, we summarize the entire algorithm and analyze the computational complexity.

3.1. Properties and Design of Equiripple FIR Filters

Besides the advantages of FIR filters, i.e., exact linear phase response and inherent stabilization, equiripple FIR filters have an explicitly specified transition width and passband/stop-band ripples (see Figure 3). As analysis filters, equiripple FIR filters can bring some important benefits, such as stop-band attenuation with a fixed maximum, the explicitly specified width of the invalid part of the sub-band spectrum (which corresponds to the transition-band spectrum) and a limited maximum deviation of the valid part of the sub-band spectrum (which corresponds to the passband spectrum). As shown in Figure 3, the specifications of a typical equiripple FIR filter consist of the passband edge ω_p, stop-band edge ω_s and maximum error in passband and stop-band δ_p, δ_s, respectively. The approximate relationship between the optimal filter length and other parameters developed by Kaiser [11] is:

$$N_f\approx \frac{-20{log}_{10}\left(\sqrt{{\delta }_p{\delta }_s}\right)-13}{14.6\mathrm{\Delta }f}+1$$

(34)

where Δf denotes the width of the transition-band,

$$\mathrm{\Delta }f=\frac{{\omega }_s-{\omega }_p}{2\pi }$$

(35)

The maximum passband variation and the minimum stop-band attenuation in decibels are defined as:

$$A_p=20{log}_{10}\left(\frac{1+{\delta }_p}{1-{\delta }_p}\right)\text{dB}$$

(36)

and:

$$A_s=-20{log}_{10}\left({\delta }_s\right)\text{dB}$$

(37)

respectively.

When the specification of a filter is explicitly specified, we can complete the design with the Parks–McClellan (PM) algorithm [30], since it is optimal with respect to the Chebyshev norm and results in about 5 dB more attenuation than the windowed design algorithm [11].

3.2. Practical Consideration of Equiripple FIR Filters

Firstly, the equiripple low-pass FIR filter is combined with a preprocessing step—complex frequency modulation—to form a passband filter for sub-band decomposition (see Figure 4).

The magnitude response of the analysis filter is shown in Figure 5, where ω_H and ω_L denote the high and the low edge of the stop-band, respectively. They satisfy:

$${\omega }_H-{\omega }_L=2{\omega }_s$$

(38)

As long as A_s is large enough and the downsample rate M meets the condition:

$$M⩾\frac{2\pi }{{\omega }_H-{\omega }_L}\Rightarrow M⩾\frac{\pi }{{\omega }_s}$$

(39)

frequency aliasing can be practically suppressed.

Secondly, due to the existence of the transition-band of each analysis filter, each sub-band spectrum contains two invalid parts. The spectral estimations of these invalid parts lead to errors. Consequently, according to [31], when mosaicking these sub-band spectral estimations into full-band, we should omit these invalid parts of spectral estimations. This procedure is illustrated in Figure 6. Thus, the composite full-band spectral estimation is practically immune to the spectral overlap.

Thirdly, due to the existence of passband ripples in equiripple FIR filters, there theoretically remains a small error in sub-band spectral estimations. Generally, by adjusting the maximum passband variation, we can limit the error to an allowable range. More precise spectral estimation necessitates compensation for the residual error. Since the ripple curve for any given equiripple FIR filter can be accurately measured, the compensation can be performed by weighting sub-band spectral estimations with the measured ripple curve.

Finally, we focus on selecting appropriate filter parameters in SDSE, which can improve the performance and reduce computational cost. The filter order should at least meet:

$$N⩾\text{max}\left(L_G,N_o-L_G\right)$$

(40)

The maximum stop-band attenuation should exceed the dynamic range of the signal to be analyzed. Once the aforementioned conditions are satisfied, the shortest transition-width can be chosen by Equation (34). Moreover, specific requirement will help to set the maximum passband variation.

3.3. Computational Complexity of SDSE

Algorithm 1 Non-overlapping sub-band spectral estimation with the steady-state Kalman predictor-based filtering strategy.

Input: The sequence ${\left\{x_n\right\}}_{n=0}^{N-1}$.

Parameters: The maximum passband variation Ap and the minimum stop-band attenuation As in decibels; the sub-band number M.

Filter Design:

Set the stop-band edge ωs by Equation (39);

Set the passband edge ωp by Equations (34), (35), (36), (37) and (40);

Design the equiripple FIR filter by Parks–McClellan (PM) algorithm [30], and then, compute the impulse response ${\left\{h_n\right\}}_{n=0}^{N_o}$ and the group delay LG.

AR Identification and Order Selection:

for pi = ps to pe (usually set ps = 1, pe ⩽ N/2 − 1) do

Estimate coefficients $\left(N{p}_i^2\right)$ of AR model by the forward and backward estimator, with flops;

Estimate the MSE $\left(\frac{N_o^4}{48}+\frac{5}{24}{N}_o^3+\frac{N_o^2}{6}\right)$ by (A10) and (A11), with flops;

end for

Select an order p by Equation (31) or Equation (33), with (N) flops.

Sequence Extrapolation:

Set the step-size k by Equation (28);

Calculate ${\left\{g_{k_l}\right\}}_{l=0}^{p-1}$ by Equations (16), (17) and (24), with ((2p)k−1) flops;

Implement forward and backward extrapolations by Equations (26) and (27), and obtain ${\left\{{\widehat{x}}_n\right\}}_{n=L_G-N_o}^{L_G+N-1}$ , with (No p) flops.

Sub-Band Spectral Estimation:

Set a rational factor $M_0=〚\frac{\pi }{{\omega }_s}〛$, where [[]] denotes a rational approximation;

for i = 1 to M do

Compute ωH and ωL by Equation (38) and ωH + ωL = (2i − 1) π/M;

Perform pre-modulation and filtering for ${\left\{{\widehat{x}}_n\right\}}_{n=L_G-N_o}^{L_G+N-1}$ by Figure 4 and Equation (10), and the computational complexity is in (2 (N + No) log (N + No)) flops;

Decimate the sequence ${\left\{{\widehat{x}}_n^{(i)}\right\}}_{n=0}^{\lceil \frac{N}{M_0}\rceil -1}$, by a factor of M0, and obtain the sub-band sequence , where ⌈ ⌉ denotes the ceiling function;

Perform spectral analysis for the sub-band sequence ${\left\{{\widehat{x}}_n^{(i)}\right\}}_{n=0}^{\lceil \frac{N}{M_0}\rceil -1}$, and denote the length of the sub-band spectrum as Ls;

Compute the length of overlapped spectrum by $\lceil \left(\frac{1}{M_0}-\frac{{\omega }_p}{\pi }\right)\frac{L_s}{2}\rceil$ and omit the overlapped parts at both the left and the right side of the sub-band spectrum.

end for

Mosaic the residual sub-band spectrums into an entire spectrum.

Output: The entire spectrum.

As shown in Algorithm 1, we summarize SDSE with the proposed filtering strategy and give the computational complexity of the major steps. First, the proposed strategy can greatly reduce the computational burden. We take the commonly-used amplitude and phase estimation (APES) [32] algorithm as an example. The full-band APES needs (N² log N) flops [33], while the computation requirement is decreased to $\left(M{\left(\lceil \frac{N}{M_0}\rceil \right)}^{2}log\left(\lceil \frac{N}{M_0}\rceil \right)\right)$ flops by SDSE with the proposed strategy. Second, except the sub-band spectral estimation, the main computation requirement is induced by the AR identification and the order selection. The computational complexity of this step is generally much lower than that of the sub-band spectral estimation. In particular, if a proper order or a small enough order interval is preselected before the AR identification, the computation of this step can be negligible.

4. Simulations and Analysis

In this section, both the feasibility and the effectiveness of the proposed strategy are evaluated by typical numerical simulations, including FIR filtering and line spectral analysis of 1D or 2D sequences.

4.1. Filtering Analysis Using the Proposed Strategy

Suppose that the input sequence {x_n} is a mixed complex exponential sequence:

$$\{\begin{array}{ll}\hfill x_n& =s_n^{\left(1\right)}+s_n^{\left(2\right)}+{\upsilon }_n\hfill \\ \hfill s_n^{\left(1\right)}& =exp\left(0.45j\pi n\right)\hfill \\ \hfill s_n^{\left(2\right)}& =\text{∑}_{l=0}^{25}100exp\left\{\left(0.55+0.035l\right)j\pi n\right\}\hfill \\ \hfill n& =0,1,\dots ,127\hfill \end{array}$$

(41)

where the measurement noise {υ_n} is a complex Gaussian process. All real parts and imaginary parts of {υ_n} are independent and identically distributed (i.i.d.) zero-mean Gaussian distributions with variance σ², i.e., Re (υ_n), Im (υ_n) ∼ (0, σ²). Our purpose is to non-distortedly extract the weak component $s_n^{\left(1\right)}$ from x_n or completely eliminate the strong component $s_n^{\left(2\right)}$.

The equiripple half-band low-pass FIR filter is chosen for the extraction. The specifications of the filter are:

$$A_p=1.4295\times {10}^{-3}\text{dB},A_{\text{s}}=81.6852\text{dB},\mathrm{\Delta }f=0.08$$

(42)

The length of the designed filter based on the given specifications is 119.

As shown in Figure 7a, the decreasing trend of the estimated residual error by the proposed strategy is consistent with the real error. When the order exceeds 57, the decrease of the estimated filtering error instantly slows down. Hence, the preferred order is 57. By comparison, due to the deficiency of the penalty strength, none of AIC and BIC can provide the right order; whereas, based on the adaptive penalty terms, both AAIC and ABIC get the right order 57 (see Figure 7b).

As shown in Figure 8b, the weak component $s_n^{\left(1\right)}$ is completely covered by the sidelobe of the out-of-band strong component $s_n^{\left(2\right)}$ thus, recognizing the existence of the weak component from the mixed spectrum is completely impossible. From the view of the magnitude response (see Figure 8a), the filter has the nominal ability of eliminating the interference of the out-of-band strong components for the in-band weak component. Due to the existence of the convolution filtering error, we still cannot find out the weak component from the convolution spectrum, as shown in Figure 8b. By contrast, once the samples contaminated by the filtering error are omitted by Equation (4) from the filtered sequence, the weak component reappears in the spectrum of the remaining samples (refer to the truncated spectrum in Figure 8b). However, the truncated spectrum has a much wider main lobe than the original spectrum, which means the spectral resolution suffers from a severe decrease. In order to simultaneously maintain the resolution and filter out the interference, we apply the proposed filtering strategy to handle the case. As shown in Figure 8c, based on the proposed strategy, the restored spectrum for the noiseless sequence closely coincides with the truth weak spectrum in shape, especially retaining the spectral resolution. In addition, even when the signal-to-noise (SNR) of $s_n^{\left(1\right)}$ is low to −3 dB (when σ² = 1), the recovery is still effective (see the magnified details of Figure 8c).

4.2. Line Spectral Analysis Using 1-D Signals

A complex exponential model can be mathematically represented as:

$$\{\begin{array}{ll}\hfill x_n& =s_n^{\left(1\right)}+s_n^{\left(2\right)}+{\upsilon }_n\hfill \\ \hfill s_n^{\left(1\right)}& =\text{∑}_{k=1}^5{\alpha }_{k}exp\left[j\left({\omega }_{k}n+{\varphi }_k\right)\right]\hfill \\ \hfill s_n^{\left(2\right)}& =\text{∑}_{i=0}^{16}100\left\{exp\left[j\left({\omega }_i^{(+)}n+{\varphi }_i^{(+)}\right)\right]+exp\left[j\left({\omega }_i^{(-)}n+{\varphi }_i^{(-)}\right)\right]\right\}\hfill \end{array}$$

(43)

(43) where:

$$\begin{array}{c}{\alpha }_1={\alpha }_3={\alpha }_5=5,{\alpha }_2={\alpha }_4=1\\ {\omega }_1=-0.075\pi ,{\omega }_2=-0.03125\pi ,{\omega }_3=0.0125\pi {\omega }_4=0.05625\pi \\ {\omega }_5=0.1\pi ,{\omega }_i^{(+)}=\left(0.15+0.05i\right)\pi ,{\omega }_i^{(-)}=-\left(0.15+0.05i\right)\pi \end{array}$$

and:

$$n=0,1,\cdots ,N-1;N=128$$

{υ_n} is a real-value sequence of i.i.d. zero-mean Gaussian random variables with variance σ² = 1.5811, i.e., u_n ∼ (0, σ²). ϕ_k, ${\varphi }_i^{+}$ and ${\varphi }_i^{-}$ are i.i.d. uniform random variables on the interval from zero to 2π, i.e., ϕ_k, ${\varphi }_i^{(+)}$, ${\varphi }_i^{(-)}$ ∼ [0, 2π).

In this case, we can get each component's SNR of $s_n^{\left(1\right)}$:

$${\text{SNR}}_1={\text{SNR}}_3={\text{SNR}}_5{=5\text{dB},\text{SNR}}_2={\text{SNR}}_4=-2\text{dB}$$

We decompose the mixed-signal x_n into four sub-bands using the proposed method with the filter parameter set as:

$$A_p=0.01\text{dB},A_s=60\text{dB},\mathrm{\Delta }f=0.05$$

(44)

The sub-band, whose radian frequency is within [−0.125π, +0.125π), is used for frequency estimation. Furthermore, we estimate the frequencies of complex sinusoids of $s_n^{\left(1\right)}$ that are contained in both mixed-signal x_n and the decomposed sub-band signal, via MUSIC, ESPRIT [34,35] and SELF-SVD [15] algorithms (see

Table 1. Comparison of frequency estimation in full-band and sub-bands.

Means and MSEs of Estimating Frequencies	Full-band			Sub-band
	MUSIC	ESPRIT	SELF-SVD	MUSIC	ESPRIT	SELF-SVD
ω̄1 (ω1 = –0.07500π)	–0.0750π	–0.0750π	–0.0790π	–0.0749π	–0.0749π	–0.0750π
ω̄2 (ω2 = –0.03125π)	–0.0309π	–0.0309π	∼	–0.0308π	–0.0313π	–0.0308π
ω̄3 (ω3 = 0.01250π)	0.0125π	0.0125π	0.0175π	0.0123π	0.0123π	0.0126π
ω̄4 (ω4 = 0.05625π)	0.0572π	0.0572π	∼	0.0555π	0.0557π	0.0555π
ω̄5 (ω5 = 0.10000π)	0.1000π	0.1000π	0.1027π	0.1000π	0.1000π	0.1000π
${\widehat{\sigma }}_{{\overline{\omega }}_1}^2\times {10}^5$	0.0072	0.0072	1.8126	0.0057	0.0058	0.0044
${\widehat{\sigma }}_{{\overline{\omega }}_2}^2\times {10}^5$	1.6462	1.6462	∼	0.1662	0.1259	0.1816
${\widehat{\sigma }}_{{\overline{\omega }}_3}^2\times {10}^5$	0.0125	0.0125	2.8847	0.0081	0.0074	0.0060
${\widehat{\sigma }}_{{\overline{\omega }}_4}^2\times {10}^5$	3.0443	3.0443	∼	0.1729	0.1199	0.1775
${\widehat{\sigma }}_{{\overline{\omega }}_5}^2\times {10}^5$	0.0360	0.0360	0.7379	0.0049	0.0041	0.0046

“∼” denotes meaningless estimates.(Monte Carlo analysis: 100 runs). SELF-SVD, singular value decomposition-based method in a selected frequency band.

). As shown in Table 1, we analyze the performance based on the Monte Carlo method. Compared with ESPRIT, SELF-SVD in full-band spectral estimation suffers from obvious performance degradations or even failures. Although SELF-SVD can theoretically attenuate the out-of-band components for the in-band frequency estimations, the ability of attenuation is not always sufficient, especially when the power of the out-of-band components are much stronger than that of the in-band components or the SNR is relatively low. Instead of performing the SVD method in the entire frequency domain as ESPRIT, SELF-SVD just performs it in the frequency interval of interest. Obviously, the remaining out-of-band interferences will be treated as in-band components, so that the frequency estimation with SELF-SVD sometimes fails. In the experiment, the power ratio of the out-of-band components to the in-band components ω₂ and ω₄ is up to 10,000 times. As a result, the corresponding frequency estimation with SELF-SVD fails to work. When we eliminate the out-of-band interferences with our method, the estimation of SELF-SVD for the residual signal exhibits similar performance as ESPRIT. In addition, MSEs of MUSIC and ESPRIT indicate that the frequency estimation in the sub-band is much more accurate than that in the full-band.

4.3. Line Spectral Analysis Using 2D signals

Let C_k (k = 1, 2, …, K) be a series of random integers with unique values generated from a uniform discrete distribution on [1,1024 × 1024]. We define two sets of nonnegative integers as:

$$\{\begin{array}{l}{p}_k=\lfloor \frac{C_k}{1024}\rfloor \hfill \\ q_k=mod\left(C_k,1024\right)\hfill \end{array}k=1,2,\cdots ,K$$

(45)

where ⌊·⌋ rounds a number to the nearest integer toward zero, and mod (·) is the modulo operator.

The 2D signal model can be expressed as:

$$\{\begin{array}{ll}{x}_{n_1,n_2}\hfill & =s_{n_1,n_2}^{\left(1\right)}+s_{n_1,n_2}^{\left(2\right)}+{\upsilon }_{n_1,n_2}\hfill \\ s_{n_1,n_2}^{\left(1\right)}\hfill & =\text{∑}_{k=1}^{K}exp\left[j2\pi \left(\frac{n_1{p}_k}{4N_1}+\frac{n_2{q}_k}{4N_2}\right)+j{\varphi }_k\right]\hfill \\ s_{n_1,n_2}^{\left(2\right)}\hfill & =\text{∑}_{l_1=0}^{16}1,000\left\{\begin{array}{c}exp\left[j2\pi \left(\frac{118.5n_1}{N_1}+\frac{\left(95.5+4l_1\right)n_2}{N_2}\right)+j{\varphi }_{l_1}^{\left(1\right)}\right]\\ +exp\left[j2\pi \left(\frac{54.5n_1}{N_1}+\frac{\left(95.5+4l_1\right)n_2}{N_2}\right)+j{\varphi }_{l_1}^{\left(2\right)}\right]\end{array}\right\}\hfill \\ \hfill & +\text{∑}_{l_2=1}^{15}1,000\left\{\begin{array}{c}exp\left[j2\pi \left(\frac{\left(54.5+4l_2\right)n_1}{N_1}+\frac{95.5n_2}{N_2}\right)+j{\varphi }_{l_2}^{\left(1\right)}\right]\\ +exp\left[j2\pi \left(\frac{\left(54.5+4l_2\right)n_1}{N_1}+\frac{159.5n_2}{N_2}\right)+j{\varphi }_{l_2}^{\left(2\right)}\right]\end{array}\right\}\hfill \end{array}$$

(46)

where:

$$n_1=0,1,\cdots ,N_1-1;n_2=0,1,\cdots ,N_2-1$$

and:

$$N_1=N_2=256,K=8,192$$

{υ_n1,n2} is a real-value sequence following υ_n1,n2 ∼ (0, σ²) with σ² = 0.005. ${\varphi }_{k,}{\varphi }_{l_1}^{\left(1\right)},{\varphi }_{l_1}^{\left(2\right)},{\varphi }_{l_2}^{\left(1\right)},{\varphi }_{l_2}^{\left(2\right)}$are uniform random variables on the interval from zero to 2π, i.e., ${\varphi }_{k,}{\varphi }_{l_1}^{\left(1\right)},{\varphi }_{l_1}^{\left(2\right)},{\varphi }_{l_2}^{\left(1\right)},{\varphi }_{l_2}^{\left(2\right)}\sim 𝒰[0,2\pi )$.

The spectrum of this 2D sequence is shown in Figure 9. Since the magnitude of $s_{n_1,n_2}^{\left(2\right)}$ is 60 dB greater than that of $s_{n_1,n_2}^{\left(1\right)}$, the sidelobe of the former significantly affects the spectral estimation of the latter. This affect is especially more severe for the components around $s_{n_1,n_2}^{\left(1\right)}$. The region inside the red pane is used to verify the performance of the proposed method.

The parameters of the analysis filter are selected as:

$$A_p=0.2\text{dB},A_s=80\text{dB},\mathrm{\Delta }f=0.05$$

(47)

The comparison of Figure 10a and 10c indicates that the Fourier spectrum of $s_{n_1,n_2}^{\left(1\right)}$ is severely affected by $s_{n_1,n_2}^{\left(2\right)}$. By contrast, the result shown in Figure 10b seems to be almost exactly the same as the desired result shown in Figure 10c. This decomposition result verifies the effectiveness of the proposed method.

To further testify the performance of our method, we select the APES [32] and the iterative adaptive approach (IAA) [36,37] for spectral estimation. Since the ideal frequency domain filters suffer from energy leakage and/or frequency aliasing problems, the APES result shown in Figure 10c is somewhat blurred. By contrast, the APES result of the decomposed sub-band based on the proposed strategy (see Figure 10d) is quite similar to the actual spectrum (see Figure 10e). Theoretically, the IAA is superior to the APES. However, as shown in Figure 10g, it is even more likely than the APES to suffer from out-of-band interferences. From the view of the sub-band IAA spectrum (see Figure 10h), most of the interferences are eliminated, while the remaining filtering error still has impacts on the spectrum. Thus, the spectral estimation experiment reveals that the sub-band decomposition based on the proposed method can provide relatively ideal performance; whereas the developed method for extrapolation is imperfect, so it can affect the performance of the IAA algorithm.

In addition, a simulated single-polarized SAR image of an airplane based on the physical and optical model is processed via the APES. The computation time of full-band APES (refer to Figure 11a) is 26.85 h, while the time of sub-band APES (refer to Figure 11b) is just 0.84 h. Obviously, the two imaging results only have tiny differences, which are hardly recognized.

5. Conclusion

This paper has investigated the problem of suppressing spectral overlap in sub-band spectral estimation. The spectral overlap phenomenon is originated from the non-ideal behavior of the analysis filtering, i.e., the filtering error. The error formation in convolution filtering was therefore discussed, based on which an extrapolation-based filtering strategy was proposed to greatly suppress spectral overlap. Several classical theories, including AR identification, Kalman prediction and the equiripple FIR filtering technique, are integrated into the strategy for linearly-optimal extrapolation. To resolve the “overfitting” in order determination with AIC and BIC, we modified the penalty terms for both criteria. The improved criteria adaptively adjust the penalty strength and avoid “overfitting” to some extent. Both 1D and 2D complex exponential signals are utilized to validate the performance of the proposed method. Moreover, we employed SAR image formation for a single-polarized SAR data, simulated based on electromagnetic theory, to testify the efficiency of our method. Future research will focus on developing more sophisticated methods for the problem of extrapolation, with which we can avoid model order determination and further improve the extrapolation precision.