AN IMPROVED METHOD FOR NONSTATIONARY SIGNALS COMPONENTS EXTRACTION BASED ON THE ICI RULE Jonatan Lerga, Victor Sucic∗ Faculty of Engineering, University of Rijeka Vukovarska 58, HR-51000 Rijeka, Croatia jlerga@riteh.hr; vsucic@riteh.hr Boualem Boashash College of Engineering, Qatar University PO Box: 2713, Doha, Qatar boualem@qu.edu.qa ABSTRACT This paper proposes an improved adaptive algorithm for com- ponents localization and extraction from a noisy multicompo- nent signal time-frequency distribution (TFD). The algorithm, based on the intersection of confidence intervals (ICI) rule, does not require any a priori knowledge of signal components and their mixture. Its efficiency is significantly enhanced by using high resolution and reduced cross-terms TFDs. The ob- tained results are compared for different signal-to-noise ra- tios (SNRs) and various time and lag window types used in the modified B-distribution (MBD) calculation, proving the method to be a valuable tool in noisy multicomponent signals components extraction in the time-frequency (TF) domain. 1. INTRODUCTION Real-life signals, such as those found in telecommunications, radar, sonar, biomedical measurements or seismology, are often nonstationary, meaning that their frequency contest exhibits time-varying properties. In practice, those signals are often also corrupted by additive noise so that signal pro- cessing algorithms performances normally depend on the signal-to-noise ratio (SNR). Since the late 1980’s, the time-frequency (TF) signal analysis has been intensively studied for the purpose of finding ade- quate tools for analyzing non-stationary signals. A special attention was dedicated to multicomponent signals TFDs due to the interactions among signal components (cross-terms) which carry redundant information and may obscure signal features of interest. An important TFD property is its TF resolution. Good time and frequency resolution seem to be conflicting require- ments, hence research efforts are still being deployed to modify existing TFDs or design new TFDs with reduced cross-terms while keeping high TF resolution [1]. Conse- quently, various reduced interference distributions (RIDs) have been proposed, one of the recent ones being the modi- fied B-distribution (MBD) [2]. ∗This work is a part of the research project ”Optimization and Design of Time-Frequency Distributions”, No. 069-0362214-1575, which is financially supported by the Ministry of Science, Education and Sports of the Republic of Croatia. Furthermore, when dealing with multicomponent signals, the components extraction often precedes other signal processing procedures, such as for example the components instanta- neous frequency (IF) estimation [3]. Numerous components extraction procedures have been proposed, often termed as blind source separation (BSS) techniques [4]. Unlike in [5] where the fixed size bandwidths based BSS algorithm was used, in the here proposed method we have in- troduced adaptivity in the bandwidth size selection. This has resulted in improved components localization and extraction from a noisy multicomponent signal TFD. 2. COMPONENTS EXTRACTION PROCEDURE In this section, we present an improved components selection and extraction algorithm from a noisy multicomponent signal based on the adaptive bandwidth size calculated for each time instant of the signal TFD. 2.1. The algorithm The algorithm (see the flowchart in Fig. 1) has several major steps: • Step 1: TFD ρ(t, f) calculation • Step 2: Finding the ρ(t, f) maximum location (t0, f0) and the TFD slice ρ(t0, f) selection • Step 3: Adaptive bandwidth detection (for each time instant t0 and each component apart) • Step 4: The selected bandwidth (defined by its bound- aries [f0(t0) − fl(t0), f0(t0) + fr(t0)]) is used to lo- calize the considered component for the time instant t0, (where f0(t0) is the frequency for which the TFD slice ρ(t0 − 1, f) reaches its maximum). Note that ρ(t0 − 1, f) is used in order to insure the IF continuity, fl(t0) and fr(t0) are the adaptive varying bandwidth lengths calculated using the intersection of confidence (ICI) rule for each time instant apart [3]. • Step 5: Once the component is localized for the time instant t0, it is extracted from the TFD so that ρ(t0, f) is set to zero for f ∈ [f0(t0)− fl(t0), f0(to) + fr(t0)]. 2011 7th International Workshop on Systems, Signal Processing and their Applications (WOSSPA) 978-1-4577-0690-5/11/$26.00 ©2011 IEEE 307 The above procedure is then repeated for each time instant of ρ(t, f) as long as the ρ(t0, f) maximum in f0(t0) vicin- ity [f0(t0) − fl(t0), f0(to) + fr(t0)] is larger than the preset threshold value c (chosen so that only a component is ex- tracted, while avoiding selecting noise as the component) de- fined as a percentage of the TFD maximum. This stage of the algorithm results in one extracted component. In order to ex- tract the remaining components, the algorithm steps 2-5 are recursively repeated, providing the remaining TFD energy is larger than the preset threshold d (defined as a percentage of the total TFD energy). 2.2. TFD choice Due to the presence of a undesirable interference terms, the choice of an appropriate TFD for a given multicomponent signal representation in the (t, f) domain is of crucial im- portance for the components extraction efficiency. Further- more, the TFD is expected to have high TF resolution. In gen- eral, there exists a tradeoff between those two TFD features [2] which led to various RIDs, one being the MBD, which was shown to outperform other fixed-kernel TFDs in terms of cross-terms reduction and resolution enhancement [2]: MBDz(t, f) =  +∞ −∞  +∞ −∞ cosh−2β(t− u)∫ +∞ −∞ cosh−2β ξdξ · ·z ( u + τ 2 ) · z∗ ( u− τ 2 ) · e−j2pifτdudτ, (1) where the parameter β controls the TF resolution and cross term suppression (0 < β ≤ 1). Although the MBD was used in this paper due to its promi- nent properties in the multicomponent signal TF representa- tion, other high TF resolution and reduced cross-terms TFDs can be also used with the signal components extraction algo- rithm proposed in this paper. 2.3. Adaptive bandwidth size selection method For the adaptive bandwidth size calculation we have used a method based on the ICI rule [6]. The method’s first step (see the flowchart in Fig. 2) is to find f0 for each TFD slice ρ(t0, f) such that the slice ρ(t0 − 1, f) reaches its maximum in (t0 − 1, f0). Once f0 is found, the method calculates the subbandwidths f+l (t0) and f+r (t0) which form the over- all bandwidth for the considered time instant t0, defined as [f0(t0)− f+l (t0), f0(t0) + f+r (t0)]. In order to get the subbandwidths fl(t0) and fr(t0), the ICI rule introduces two sequences of growing subbandwidths lengths to the left hand and the right hand side of f0, denoted as Kl and Kr, where Hl(t0, kl) = {fl(t0, 1) < fl(t0, 2) < · · · < fl(t0, Kl)} and Hr(t0, kr) = {fr(t0, 1) < fr(t0, 2) < · · · < fr(t0, Kr)}. The accompanying confidence intervals are: Dl(t0, kl) = [Ll(t0, kl), Ul(t0, kl)], 1 ≤ kl ≤ Kl, (2) Dr(t0, kr) = [Lr(t0, kr), Ur(t0, kr)], 1 ≤ kr ≤ Kr, (3) TFD ρ(t,f) calculation Finding the maximum location (t0-1,f0(t0)) Multicomponent signal All components are extracted? No Estimated components IFs Yes Adaptive bandwidth calculation (f0(t0)-fl(t0),f0(t0)+fr(t0)) Component extraction for time instant t0 t0→ t0-1 t0→ t0+1 max(ρ(t0,f)) in f0(t0) vicinity>εc Yes No max(ρ(t0,f)) in f0(t0) vicinity>εc Yes No ρ(t,f) → ρ(t,f) without extracted component TFD slice ρ(t0,f) selection Fig. 1. Components extraction flowchart. where the lower limits Ll(t0, kl) and Lr(t0, kr) and the upper limits Ul(t0, kl) and Ur(t0, kr) (the indices l and r stand for the left hand and the right hand side of f0), are respectively defined as: Ll(t0, kl) = ρ(t0, f0)− Γσ(t0, kl), (4) Lr(t0, kr) = ρ(t0, f0)− Γσ(t0, kr), (5) Ul(t0, kl) = ρ(t0, f0) + Γσ(t0, kl), (6) Ur(t0, kr) = ρ(t0, f0) + Γσ(t0, kr). (7) The parameter Γ regulates the confidence interval width (where the probability P (Γ) → 1 as Γ grows), and σ(t0, k) is the standard deviation of the estimated ρˆ(t0, f0) value [6]. The ICI rule then tracks the values of the largest lower and the smallest upper confidence intervals limits for each sequence of confidence intervals Dl(t0, kl) and Dr(t0, kr), defined as: Ll(t0, kl) = maxi=1,...,klLi(t0, kl), (8) U l(t0, kl) = mini=1,...,klUi(t0, kl), (9) Lr(t0, kr) = maxi=1,...,krLi(t0, kr), (10) Ur(t0, kr) = mini=1,...,krUi(t0, kr), (11) 308 Finding the ρ(t0-1,f) maximum location (t0,f0) TFD ρ(t0,f) slice Sequence of the smallest upper intervals limits calculation Ul(t0,kl)=mini=1,…,kl Ui(t0,kl) Overall bandwidth width for time instant t0 calculation (f0(t0)-fl(t0),f0(t0)+fr(t0)) Sequence of Kl confidence intervals calculation Dl(t0,kl)=[Ll(t0,kl),Ul(t0,kl)] Sequence of Kr confidence intervals calculation Dr=[Lr(t0,kr),Ur(t0,kr)] Sequence of the largest lower intervals limits calculation Ll(t0,kl)=maxi=1,…,kl Li(t0,kl) Setting fl(t0) as the largest bandwidth width for which L(t0,kl)≤U(t0,kl) Sequence of the smallest upper intervals limits calculation Ur(t0,kr)=mini=1,…,kr Ui(t0,kr) Sequence of the largest lower intervals limits calculation Lr(t0,kr)=maxi=1,…,kr Li(t0,kr) Setting fr(t0) as the largest bandwidth width for which L(t0,kr)≤U(t0,kr) Extracted component for time instant t0 Fig. 2. Adaptive bandwidths selection flowchart. giving the subbandwidth size f+l (t0) as the largest one for which it is still true that: Ll(t0, kl) ≤ U l(t0, kl), (12) and the subbandwidth size f+r (t0) as the largest one for which it is true that: Lr(t0, kr) ≤ U r(t0, kr). (13) The f+l (t0) and f+r (t0) were shown to be the subbandwidth lengths closest to the optimal lengths f∗l (t0) and f∗r (t0), re- spectively, giving the optimal bias to variance tradeoff and minimizing the estimation mean square error (MSE) [3]. 3. SIMULATION RESULTS AND DISCUSSION The proposed method for components localization and ex- traction was applied to a noisy three component signal (with the length N = 128) of the form x(n) = z1(n) + z2(n) + z3(n) + (n), where zm(n) = Am exp(jφm(n)) (Am = 1) is a signal component, and (n) is complex-valued additive Gaussian noise. The signal contains two sinusoidal frequency modulated (FM) components and one linear FM component with different time supports (which partially overlap). The IF laws of the components are: ω1(n) = 0.35+0.05 cos(2pi(n− N1/2)/N1− pi/2), ω2(n) = 0.05+0.3(n−1)/(N2−1), and ω3(n) = 0.075 + 0.025 cos(2pi(n − N3/2)/N3 − pi/2) (the component lengths are N1 = 96, N2 = 48, and N3 = 48). (a) (b) (c) (d) Fig. 3. (a) Noisy signal components extraction (MBD calcu- lated using rectangular time and lag windows of length h = 9, SNR=10, Γ = 25, c = 0.25, d = 0.25, N = 128, β = 0.1 as in [2]). (b) First component extracted. (c) Second compo- nent extracted. (d) Third component extracted. 20 40 60 80 100 120 20 40 60 80 100 TFD time and frequency smoothing windows lengths Remaining TFD energy [% ] Rectangular Hamming Hanning Fig. 4. Remaining MBD energy after components extraction as a function of the MBD time and lag window lengths (Γ = 25, SNR=10, c = 0.25, N = 128, m = 100, β = 0.1 as in [2]). The method’s component extraction efficiency has been an- alyzed for various time and lag window types and lengths h used in the MBD calculation, as well as for different SNRs. The results of the components extraction using the proposed adaptive method are given in Fig. 3, showing that each com- ponent was precisely localized and extracted. However, the method performance is effected by the TFD cross-terms sup- pression and TF resolution capabilities which depend on the MBD time and lag window lengths [2]. The remaining MBD energy after the components extraction (averaged over m = 309 5 10 15 20 25 30 35 40 45 50 0 20 40 60 80 ? Remaining TFD energy [% ] Rectangular Hamming Hanning Fig. 5. Remaining MBD energy after components extraction as a function of Γ (MBD calculated using the time and lag windows of length h = 9, SNR=10, c = 0.25, m = 100, β = 0.1 as in [2]). 0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 10 Frequency ?(t 0,f ) ?(t0,f) ?(t0,f0) Extr. comp. at t0 (?=35) Extr. comp. at t0 (?=25) Extr. comp. at t0 (?=15) Fig. 6. Slice ρ(t0, f ) and the component extraction for the time instant t0 and different Γ values. Table 1. Remaining TFD energy after components extraction for different SNRs and Γ values (MBD calculated using the rectangular time and lag windows of length h = 9, c = 0.25, N = 128, m = 100, β = 0.1 as in [2]). Γ SNR [%] 2 5 10 15 5 72.49 65.67 61.94 70.96 10 56.37 44.79 42.73 54.57 15 44.79 35.26 30.30 41.62 20 38.86 26.46 23.23 32.52 25 35.34 24.90 18.82 27.31 50 25.54 22.34 9.70 13.80 100 Monte Carlo simulation runs) as a function of the MBD time and lag window lengths and the ICI parameter Γ for the rectangular, Hamming, and Hanning window types is shown in Figs. 4 and 5, respectively. Table 1 gives the results of the remaining TFD energy after the components extraction for different SNRs and Γ values. The choice of Γ affects the time-varying bandwidth size se- lection, and thus the component extraction quality, as shown in Fig. 6. Too small Γ values give undersized bandwidth lengths which result in omitting some parts of the considered component, hence increasing the remaining TFD energy af- ter the components extraction. On the other hand, too large Γ values result in oversized bandwidth lengths which may lead to inaccurate component selection (e.g. selecting seg- ments of neighboring components), although the remaining MBD energy after the extraction would be smaller than the one obtained using a smaller Γ value. Thus, when evaluating the method performance, beside the remaining TFD energy after the component extraction, one should also take into con- sideration whether the components are correctly localized, as ensured by this proposed algorithm. 4. CONCLUSION This paper presents an automatic adaptive method to localize and extract signal components from a noisy multicomponent signal TFD. The method is based on the asymmetrical adap- tive bandwidth selection for each time instant and each com- ponent apart. In order to get the proper bandwidth size, the ICI rule was used. The method’s performance was analyzed for different SNRs and various time and lag window lengths and types used in the MBD calculation. The presented results show that the method is an efficient tool for automatic com- ponents extraction of noisy multicomponent signals in the TF domain. 5. REFERENCES [1] B. Boashash, Ed., Time Frequency Signal Analysis and Processing: A Comprehensive Reference, Elsevier, Ox- ford, UK, 2003. [2] Z. M. Hussain and B. Boashash, “Adaptive instantaneous frequency estimation of multicomponent FM signals us- ing quadratic time-frequency distributions,” Signal Pro- cessing, IEEE Transactions on, vol. 50, no. 8, pp. 1866– 1876, Aug 2002. [3] J. Lerga and V. Sucic, “Nonlinear IF estimation based on the pseudo WVD adapted using the improved sliding pairwise ICI rule,” Signal Processing Letters, IEEE, vol. 16, no. 11, pp. 953–956, Nov 2009. [4] B. Barkat and K. Abed-Meraim, “Algorithms for blind components separation and extraction from the time- frequency distribution of their mixture,” EURASIP J. Appl. Signal Process., vol. 2004, pp. 2025–2033, 2004. [5] J Lerga, V. Sucic, and B. Boashash, “An efficient algo- rithm for instantaneous frequency estimation of nonsta- tionary multicomponent signals in low SNR,” EURASIP J. Adv. Signal Process., vol. 2011, 2011. [6] J. Lerga, M. Vrankic, and V. Sucic, “A signal denoising method based on the improved ICI rule,” Signal Process- ing Letters, IEEE, vol. 15, pp. 601–604, 2008. 310