Simultaneous blind signal separation and denoising

Simultaneous Blind Signal Separation and Denoising H. Hammam*, A. E. Abu El-Azm**, M. E. Elhalawany** and F. E. Abd El-Samie** *Telecom Egypt, Alexandria, Egypt, hossammha@yahoo.com E-mail: hossammha@yahoo.com ** Department of Electronics and Electrical Communications, Faculty of Electronic Engineering, Menoufia University, Menouf, Egypt. E-mails: halawany_m85@yahoo.com, abouelazm_atef@yahoo.com, fathi_sayed@yahoo.com p x1 (k ) = ∑ h11 (i ) s1 (k − i ) i=0 p + ∑ h12 (i ) s 2 (k − i ) + v1 (k ) i=0 p x 2 (k ) = ∑ h21 (i ) s1 (k − i ) i=0 p + ∑ h22 (i ) s 2 ( k − i ) + v 2 ( k ) i=0 Abstract -- This paper deals with the problem of blind separation of noisy signals. It proposes the application of the discrete wavelet transform as a preprocessing step of denoising prior to the blind separation algorithm. The separation of sinusoidal signals as well as speech signals in a noisy environment is studied with and without the use of the wavelet denoising step. The simulation results confirm the usefulness of this step. In all our simulation experiments, soft thresholding is used. I. INTRODUCTION Blind signal separation is an important branch of signal processing. It deals mainly with mixed signals which are frequently encountered in real life. Real life signals are frequently mixed with undesired signals. This fact has motivated the evolution of blind signal separation algorithms. The word blind means that there is no a priori information about the mixed signals and there sources. Several approaches have been proposed for blind signal separation [1]-[8]. Some of these approaches depend on independent component analysis [1]. Others depend on higher order statistics [2]. There is also a category of correlation based separation algorithms [3]. Adaptive signal separation algorithms also exist [4]. Generally, most of the separation algorithms deal with the problem in the presence of noise. The existence of noise causes some deterioration in the quality of the separated signals. Thus, there is a bad need for a pre-processing step of denoising to reduce the effect of noise in the separation process. Wavelet denoising is well known as an effective denoising approach [9],[10]. This denoising approach is based on the decomposition of the noisy signal into sub-bands and neglecting samples of less importance in the high frequency sub-bands. These samples may be due to noise. This process has the effect of reducing noise power on the signal. The signal decomposed into sub-bands is then reconstructed again using an inverse wavelet transform. This paper proposes the application of wavelet denoising on the noisy mixtures of signals and then the application of the blind signal separation algorithm. The rest of the paper is organized as follows. Sections II and III present the blind signal separation algorithm. Section IV gives the principles of wavelet denoising. Section V gives the simulation results. Finally, the concluding remarks are given in section VI. II. PROBLEM FORMULATION or in matrix form as follows: T T ⎞ s (k ) ⎛ ⎞ ⎛ v (k ) ⎞ ⎛ x1 (k ) ⎞ ⎛ h11 h12 ⎟⎜ 1 ⎟ + ⎜ 1 ⎟ ⎜⎜ ⎟⎟ = ⎜ T ⎟ ⎜ ⎟ ⎜ T ⎟ ⎜ ⎝ x 2 (k ) ⎠ ⎝ h 21 h 22 ⎠⎝ s 2 (k ) ⎠ ⎝ v 2 (k ) ⎠ (2) where [ h Tij = hij (0),…, hij ( p) ] s Ti (k ) = [si (k ),…, si (k − p)] (3) v1(k) and v2(k) are due to noise, hij is the impulse response from source j to sensor i and p is the order of the filter. For simplicity, the source signals are assumed to be of zero mean and statistically independent. From Eqs. 1 and 2, it is clear that the mixtures are convolutive sums of sources in the presence of noise. The problem is simplified by assuming that the signals arrive at the sensors unfiltered which is equivalent to setting h11=h22=1. Taking Z-transform of equation (2) and neglecting the effect of noise lead to ⎛ X 1 ( z ) ⎞ ⎛⎜ H 11 ( z ) H 12 ( z ) ⎞⎟⎛ S1 ( z ) ⎞ ⎜⎜ ⎟⎟ = ⎜ ⎟⎟ ⎜ ⎟⎜ ⎝ X 2 ( z ) ⎠ ⎝ H 21 ( z ) H 22 ( z ) ⎠⎝ S 2 ( z ) ⎠ (4) This model can be represented by the block diagram of Fig. 1. Simplifying Eq. 4 leads to [11]: Blind signal separation deals with mixtures of signals in the presence of noise. If there are two signal sources s1(k) and s2(k) and two observations x1(k) and x2(k), these observations can be represented as follows [11]: 978-1-4244-2116-9/08/$25.00 ©2008 IEEE (1) 107 ′ ( z ) ⎞⎛ S1′ ( z ) ⎞ ⎛ X 1 ( z ) ⎞ ⎛⎜ 1 H 12 ⎟⎜ ⎜⎜ ⎟⎟ = ⎟ ⎜ ′ ( z) 1 ⎟⎠⎜⎝ S 2′ ( z ) ⎟⎠ ⎝ X 2 ( z ) ⎠ ⎝ H 21 (5) s1(k) s2(k) + H11(z) x1(k) s1(k) + x1(k) + H21(z) H2(z) W2(z) H12(z) H1(z) W1(z) H22(z) + x2(k) s2(k) S1′ ( z ) = H11 ( z ) S1 ( z ) S 2′ ( z ) = H 22 ( z ) S 2 ( z ) ′ ( z) = H12 H12 ( z ) H 22 ( z ) ′ ( z) = H12 H12 ( z ) H 22 ( z ) + x2(k) y2(k) Fig. 2, Schematic diagram of the 2x2 separation algorithm. Fig. 1, A fully coupled 2 x 2 mixing system. where + y1(k) From Eq. 9, It is clear that the solution of the problem is to find a suitable W1(z) and W2(z) such that Y1(z) and Y2(z) each contains only S1(z) or S2(z). This is achieved only if either the diagonal or the anti-diagonal are zero. Assuming s1(k) and s2(k) are stationary, zero mean and (6) independent random signals, the cross-correlation of the two signals is equal to zero, that is: rs1s2 (l ) = E[ s1 (k ) s 2 (k + l )] = 0 ∀ l (11) For Hii(z)=1, which is the case of interest, Eq. 5 simplifies to: ⎛ X 1 ( z ) ⎞ ⎛⎜ 1 H 12 ( z ) ⎞⎟⎛ S1 ( z ) ⎞ ⎜ ⎟ ⎜⎜ ⎟⎟ = ⎜ 1 ⎟⎠⎜⎝ S 2 ( z ) ⎟⎠ ⎝ X 2 ( z ) ⎠ ⎝ H 21 ( z ) If y1(k) and y2(k) each contain components of s1(k) or s2(k) only, then the cross-correlation of y1(k) and y2(k) should also (7) be zero, then ry1 y2 (l ) = E[ y1 (k ) y 2 (k + l )] = 0 ∀ l (12) The objective of blind signal separation is to get the signals y1(k) and y2(k) which are as close as possible to x1(k) and x2(k). We can assume that: W1 ( z ) ⎞⎛ X 1 ( z ) ⎞ ⎛ Y1 ( z ) ⎞ ⎛ 1 ⎜⎜ ⎟⎟ = ⎜⎜ ⎟⎜ ⎟ 1 ⎟⎠⎜⎝ X 2 ( z ) ⎟⎠ ⎝ Y2 ( z ) ⎠ ⎝W2 ( z ) substituting Eq. 8 into Eq. 12 gives: ry1 y2 (l ) = E[( x1 (k ) + w 1T x 2 (k ))( x 2 (k + l ) + w T2 x1 (k + 1))] (8) If rxi x j (l ) = E[ xi (k ) x j (k + l )] ; Eq. (13) becomes [11]: where w Ti = [wi (0),… , wi (q)] x Ti (k ) = [xi (k ),… , xi (k − q)] ry1 y2 (l ) = (9) ⎛ Y1 ( z ) ⎞ ⎛1 + W1 ( z ) H 21 ( z ) W1 ( z ) + H12 ( z ) ⎞⎛ S1 ( z ) ⎞ ⎟⎟ = ⎜⎜ ⎟⎟⎜⎜ ⎟⎟ ⎜⎜ ⎝ Y2 ( z ) ⎠ ⎝ W2 ( z ) + H 21 ( z ) 1 + W2 ( z ) H 12 ( z ) ⎠⎝ S 2 ( z ) ⎠ (10) rx1x2 (l ) + w1T ⎛ rx x (l ) ⎞ 2 2 ⎜ ⎟ ⎜ ⎟ ⎜⎜ ⎟⎟ ⎝ rx2 x2 (l + q) ⎠ ⎛ rx x (l ) ⎞ ⎜ 11 ⎟ T ⎟ + w 1T R x2 x1 (l )w 2 + w2 ⎜ ⎜⎜ ⎟⎟ ⎝ rx1x1 (l + q) ⎠ Substituting Eq. 7 into Eq. 8 leads to [11]: III. BLOCK BASED SEPARATION ALGORITHM (13) (14) where R x2 x1 (l ) = E[x 2 (k )(x1 (k + l )) T ] is a (q+1) x (q+1) matrix which is a function of the cross-correlations of x1 and This section deals with a time domain iterative separation x2. The cost function C is defined as the sum of the squares of algorithm for the 2 x 2 convolutive system. The algorithm minimizes the output cross-correlations for an arbitrary the cross-correlations between the two inputs as follows [11]: number of lags with q+1 tap FIR filters. ∑r l2 C= l =l1 108 y1 y2 (l ) 2 (15) where l1 and l2 constitute a range of cross-correlation lags. C can be also expressed as: C = r y1 y2 T r y1 y2 (16) where r y1 y2 = [ry1 y2 (l1 ), , ry1 y2 (l 2 )]T (17) between adjacent filters is usually 2:1 leading to octave bandwidths and center frequencies that are one octave apart [9], [10]. The outputs of the filters are usually maximally decimated so that the number of DWT output samples equals the number of input samples and the transform is invertible as shown in Fig.4. A. IMPLEMENTATION OF THE DWT The art of finding a good wavelet lies in the design of the set of filters, H0, H1, G0 and G1 to achieve various tradeoffs r y1 y2 = rx1x2 + [Q +x2 x2 ]T w 1 + [Q −x1x1 ]T w 2 between spatial and frequency domain characteristics while (18) satisfying the perfect reconstruction (PR) condition [10]. In T + R x2 x1 A(w 2 )w 1 Fig. 4, the process of decimation and interpolation by 2:1 at or the output of H0 and H1 sets all odd samples of these signals T T + − to zero. For the low pass branch, this is equivalent to r y1 y2 = rx1x2 + [Q x2 x2 ] w 1 + [Q x1x1 ] w 2 1 (19) 1 + (− 1)n . Hence, X 0 ( z ) is multiplying x0 (n) by + R Tx1x2 A(w 1 )w 2 2 converted to {X 0 ( z ) + X 0 (− z )} . Similarly, X 1 ( z ) is converted where Q +x2 x2 and Q −x1x1 are (q+1) x (l2-l1+1) matrices, Rx2x1 is 1 to {X 1 ( z ) + X 1 (− z )} . 2 a (2q+1) x (l2-l1+1) matrix. These are all correlation matrices Thus, the expression for Y (z ) will be given by [10]: of x1 and x2 and are estimated using sample correlation estimates. A(w1) and A(w2) are (2q+1) x (q+1) matrices which contain w1 and w2, respectively. In order to find some 1 Y ( z ) = {X 0 ( z ) + X 0 (− z )}G0 ( z ) suitable w1 and w2, C is minimized such that: 2 1 ∂C + {X 1 ( z ) + X 1 (− z )}G1 ( z ) = [0, ,0]T , k = 1,2 (20) 2 ∂w k (25) ⎧ H 0 ( z )G0 ( z ) ⎫ 1 Let = X ( z )⎨ ⎬ + H 1 ( z )G1 ( z )⎭ 2 ⎩ ψ 1 = ([Q +x2 x2 ]T + R Tx2 x1 A (w 2 )) (21) ⎫ ⎧ H 0 (− z )G0 ( z ) 1 + X (− z )⎨ ψ 2 = ([Q −x1x1 ]T + R Tx1x2 A(w1 )) ⎬ ( ) ( ) + − H z G z 2 1 1 ⎭ ⎩ Thus: ( Substituting Eq. 21 into Eqs. 18 and 19 gives: r y1 y2 = rx1x2 + ψ 1w1 + [Q −x1x1 ]T w 2 or r y1 y2 = rx1x2 + ψ 2 w 2 + [Q +x2 x2 ]T w 1 The first PR condition requires aliasing cancellation and forces the above term in X (− z ) to be zero. (22) Hence, {H 0 (− z )G0 ( z ) + H1 (− z )G1 ( z )} = 0 , Which can be achieved if : (23) H1 ( z ) = z − k G0 (− z ) From Eq. 20, we obtain [11]: w2 = G1 ( z ) = z k H 0 (− z ) and w 1 = −(ψ 1Tψ 1 ) −1ψ 1T (rx1x2 + [Q −x1x1 ]T w 2 ) −(ψ 2Tψ 2 ) −1ψ 2T ) (rx1x2 + [Q +x2 x2 ]T (26) (24) where k must be odd (usually k = ±1). w1 ) w1 and w2 can be obtained by iterating between the two equations until convergence is achieved when the rate of change of parameter values is less than a pre-set threshold. By estimating w1 and w2, we then obtain a set of outputs y1(k) and y2(k). Each output contains s1(k) or s2(k) only. IV. THE DISCRETE WAVELET TRANSFORM X0(z) H0(z) (1/2){X0(z)+ X0(-z)} ↓2 ↑2 G0(z) + X(z) Wavelets have become a popular tool in most signal processing applications such as data fusion, denoising, compression and restoration. The conventional discrete wavelet transform (DWT) may be regarded as equivalent to filtering the input signal with a bank of bandpass filters, whose impulse responses are all approximately given by scaled versions of a mother wavelet. The scaling factor 109 H1(z) ↓2 X1(z) ↑2 G1(z) (1/2){X1(z)+ X1(-z)} Fig. 3, Perfect reconstruction wavelet filter bank. Y(z) ⎧ x ⎪2 x − T ⎪ f soft ( x) = ⎨ ⎪T + 2 x ⎪⎩ 0 The second PR condition is that the transfer function from X (z ) to Y (z ) should be unity: {H 0 ( z )G0 ( z ) + H1 ( z )G1 ( z )} = 2 (27) If we define a product filter P( z ) = H 0 ( z )G0 ( z ) and substitute from Eq. 26 into equation (27), then the PR condition becomes: H 0 ( z)G0 ( z) + H1 ( z)G1 ( z) = P( z) + P(− z) = 2 (28) This needs to be true for all z and, since the odd powers of z in P(z ) cancel with those in P(− z ) , it requires that p0 = 1 and that pn = 0 for all n even and non-zero. The polynomial P(z ) should be a zero phase polynomial to minimize distortion. In general, P(z ) is of the following form: P (z p1 z )= −1 + p 5 z 5 + p 3 z 3 + p1z + 1 + + p3z −3 + p5z −5 + (29) The design method for the PR filters can be summarized in the following steps: 1- Choose p1 , p3 , p5 , to give zero phase polynomial P(z ) with good characteristics. 2- Factorize P(z ) into H 0 ( z ) and G0 ( z ) with similar lowpass frequency response. 3- Calculate H1 ( z ) and G1 ( z ) from H 0 ( z ) and G0 ( z ) . To simplify this procedure, we can use the following relation: P(z ) = Pt (Z ) = 1 + pt ,1 Z + pt ,3 Z 3 + pt ,5 Z 5 + x ≥T T /2≤ x <T − T < x ≤ −T / 2 x <T /2 T denotes the threshold value and x represents the coefficients of the high frequency components of the DWT. V. RESULTS AND DISCUSSION In this section, two experiments are conducted to test the effect of wavelet denoising on blind signal separation. The first experiment is performed on two mixtures of two sinusoidal signals contaminated by white Gaussian noise. The blind signal separation algorithm is tested with and without the preprocessing denoising step. The results of this experiment shown in Fig. 4 reveal the effect of the wavelet denoising step. It is clear that a better quality of separated signals can be obtained with wavelet denoising. In the second experiment, a more realistic case is studied. This is the speech separation case. Two mixtures of noisy speech signals are used in this experiment and the results are shown in Fig. 5. These results also reveal the necessity of the wavelet denoising step prior to the blind separation. In both experiments, the Haar wavelet transform is used. The method of denoising used is the soft thresholding as it has a better performance than the hard thresholding. (30) (a) Noisy signal 1. where Z= ( 1 z + z −1 2 ) (31) B. WAVELET DENOISING Wavelet denoising is a simple operation which aims at reducing noise in a noisy signal. It is performed by choosing a threshold that is sufficiently large multiple of the standard deviation of the noise in the image. Most of the noise power is removed by thresholding of wavelet transform values. There are two types of thresholding; hard and soft thresholding [9], [10]. The equation of the hard thresholding is given by: ⎧x ⎪ f hard ( x) = ⎨ ⎪0 ⎩ (33) (b) Noisy signal 2. (c) Mixed signal 1. x ≥T (32) x <T (d) Mixed signal 2. On the other hand, that of soft thresholding is given by: 110 (e) Separated signal 1 without denoising. (d) Noisy mixture signal 2. (f) Separated signal 2 without denoising. (e) Separated speech signal 1 without denoising. (g) Separated signal 1 with denoising. (f) Separated speech signal 2 without denoising. (h) Separated signal 2 with denoising. Fig. 4. Separation of noisy sinusoidal signals with and without wavelet denoising (soft thresholding). (g) Separated speech signal 1 with denoising. (a) Noisy speech signal 1. (h) Separated speech signal 2 with denoising. Fig. 5. Separation of noisy speech signals with and without wavelet denoising (soft thresholding). VI. CONCLUSION (b) Noisy speech signal 2. In this paper, wavelet denoising is introduced to achieve a better performance of blind signal separation algorithms. In order to show the effectiveness of the proposed wavelet denoising step, two experiments have been presented. The experimental results show the effectiveness of the wavelet denoising step in the blind separation of noisy sinusoidal signal as well as speech signals. (c) Noisy mixture signal 1. 111 REFERENCES [1] [2] [3] [4] [5] [6] T. W. Lee, M. Girolami, and T. J. Sejnowski, “Independent Component Analysis Using an Extended J. F. Cardoso, “Blind Signal Separation: Statistical Infomax Algorithm for Mixed Sub-Gaussian and SuperPrinciples,” Proc. IEEE, vol. 86, pp. 2009-2025, Oct. Gaussian Sources,” Neural Computation, vol. 11, pp. 1998. 409-433, 1999. A. Belouchrani, K. Abed-Meraim, J. F. Cardoso, and E. [7] S. Amari, and J. F. Cardoso, “Blind Source SeparationMoulines, “A Blind Source Separation Technique Using Semiparametric Statistical Approach,” IEEE Trans. Second-Order Statistics,” IEEE Trans. Signal Signal Processing, vol. 45, pp. 2692-2700, Nov. 1997. Processing, vol. 45, pp. 434-444, Feb. 1997. [8] X. R. Cao, and R. W. Liu, “General Approach to Blind D.C.B. Chan, P.J.W. Rayner, and S.J. Godsill, "MultiSource Separation,” IEEE Trans. Signal Processing, vol. channel Signal Separation," Proc. ICASSP, vol. II, pp. 44, pp. 562-571, Mar. 1996. 649-652, Atlanta, GA, May 1996. [9] J. S. Walker, A. Primer, "Wavelets and Their Scientific J. F. Cardoso, and B. H. Laheld, “Equivariant Adaptive Applications," CRC Press LLC, 1999. Source Separation,” IEEE Trans. Signal Processing, vol. [10] A. Prochazka, J. Uhlir, P. J. W. Rayner and N. J. 44, pp. 3017-3030, Dec. 1996. Kingsbury, "Signal Analysis and Prediction," Birkhauser A. J. Bell and T. J. Sejnowski, “An InformationInc. , 1998. Maximization Approach to Blind Separation and Blind [11] D. C. Chan, “Blind Signal Separation” A PhD Deconvolution,” Neural Computation, vol. 7, pp. 1129dissertation, University of Cambridge, January, 1997. 1159, 1995. 112