Statistical Multirate High-Resolution Signal Reconstruction using the EMD-IT Based Denoising Approach

The reconstruction problem of a high-resolution (HR) signal from a set of its noise-corrupted low-resolution (LR) versions is considered. As a part of this problem, a hybrid method that consists of four operation units is proposed. The first unit applies noise reduction based on the empirical mode decomposition interval-thresholding to the noisy LR observations. In the second unit, estimates of zero-interpolated HR signals are obtained by performing up-sampling and then time shifting on each noise reduced LR signal. The third unit combines the zero-interpolated HR signals for attaining one HR signal. To eliminate the ripple effect, finally, median filtering is applied to the resulting reconstructed signal. As compared to the work that employs linear periodically time-varying Wiener filters, the proposed method does not require any correlation information about desired signal and LR observations. The validity of the proposed method is demonstrated by several simulation examples.


Introduction
Multirate statistical signal processing have found great interest over the last two decades for coping with the problems such as spectrum estimation, prediction, sensor fusion, time-delay estimation, and reconstruction or estimation of a stationary random process from multiple observations measured at different sampling rates [1].For the solution of these problems, different methods are suggested by several researchers.Because the focus of our study is concerned with high-resolution (HR) signal reconstruction, in this paper, we shall concentrate on the methods related to this topic.
Using the principle of maximum entropy (ME), in [2], an optimal least-mean-squares estimator is developed for estimating the samples of a wide-sense stationary (WSS) random signal from multiple observations at low sampling rates.Within the scope of the reconstruction or estimation of a WSS random signal through its noisy LR measurements observed at different sampling rates, studies that employ the optimal linear filtering referred as Wiener filtering are carried out under the guidance of Therrien [3][4][5][6][7][8][9].The problem of estimating a WSS random signal from two sequences observed separately at full-rate and half-rate, and both of which are affected by measurement noise is investigated in [3].Using the least-squares (LS) approach, in [4], three optimal filter structures are proposed for the estimation of the HR random process from two observations, one of which is measured at full-rate with a low signal-to-noise ratio (SNR) and the other is measured at arbitrary low-rate with a high SNR.A generalized version of [3] is considered in [5] that is a preliminary work of dissertation presented in [6] where optimal multirate filtering procedure is formulated taking the second-order statistics into consideration for the goal of estimating an HR signal from LR measurements observed at different sampling rates.In [7], derivations in [5] and [6] are adapted to reconstruct an HR signal from a set of its noisy and distorted LR versions.Eventually, on the basis of the LS approach, formulations of multirate optimal filtering derived in [5][6][7] are extended to the two-dimensional case [8] (refer to [9] for detailed information).It is remarked that both ME and Wiener filter based methods require the knowledge of second-order statistics related to the desired HR signal and the noisy LR observations.In addition to these mentioned two main groups of methods, recently, adaptive filtering and compressive sensing based approaches dealing with the statistical multirate signal estimation are suggested [10][11][12].
This paper aims at providing the reconstruction of an HR signal from a set of its noisy LR observations, without requiring the knowledge of any second-order statistics.In this context, a hybrid method consisting of four operation units is proposed.The paper is organized as follows.The statement of the problem and its relation to prior work are addressed in Sec. 2. Section 3 is devoted to the presentation of the proposed HR signal reconstruction method.To demonstrate the validity of the proposed method, two simulation examples are considered in Sec. 4. Paper ends up with concluding remarks.

Problem Statement and Relation to Prior Work
The problem studied here is concerned with the reconstruction of an HR unobservable random signal d[n] from a set of its LR measurements {x i [m i ]|i = 0,1,…,M 1} corrupted by additive white Gaussian noise (AWGN).In other words, these measurements are noisy, time shifted and down-sampled versions of d [n] and are generated by the model shown in Fig. 1, where AWGN, time shifting and down-sampling factor are indexed to the particular signal of interest [9].Thus, the ith LR observation signal   One of the existing methods to solve the problem emphasized above employs linear periodically time-varying (LPTV) Wiener filters that are optimal in mean-squares sense [7].As shown in Fig. 2a, this method requires knowing down-sampling rate L, LPTV Wiener filter order P, cross-correlations between desired HR signal and LR signals, and cross-correlations among LR signals.However, the proposed method does not require any information related to the second order statistics of d[n] and LR signals {x i [m]}.To know down-sampling rate L alone is enough for performing signal reconstruction task (see Fig. 2b).

Proposed Method
Under the assumption of down-sampling rate L is known, the method proposed for the aim of reconstructing a desired HR signal from a set of its noisy LR observations consists of four signal processing units.
The first signal processing unit serves the purpose of reducing the noise effect on noisy LR observations.For this aim, empirical mode decomposition (EMD) based denoising approach is used.As a part of Hilbert-Huang transform, EMD is a method of decomposing any complicated signal into a monotonic residual component and a series of intrinsic mode functions (IMFs) with degree of frequencies in descending order [13], [14].Note that the first IMF is related to the highest-frequency sequence of the ith noisy LR signal x i [m], [14].Taking account of the fact that the frequency components of AWGN spread out HR signal reconstruction based on the method of [7].over the entire frequency range, the first IMF of each noisy LR observation signal can be considered as it is mainly associated with the AWGN [15].Thus, a noise reduced version of the signal can be obtained by extracting its first IMF from the noisy signal.Such an attempt has recently been made in [16] for the purpose of reducing the AWGN effect on each of the noisy LR observation signals before using them with the LPTV Wiener filter structures for carrying out the HR signal reconstruction task.

] [n d
Nevertheless, there are still effects of noise on other IMFs.That is, separating the first IMF from noisy signal alone may not usually be sufficient for filtering noise in an effective manner.Taking this fact into consideration, in the first subunit of the proposed method, the noise reduced versions of the noisy LR observations are obtained on the basis of the EMD interval thresholding (EMD-IT) algorithm [17].In the framework of this algorithm, each noisy LR signal is decomposed into a monotonic residual sequence r (Ki) [m] and K i empirical modes termed IMFs {c i (k) [m]|k = 1,2,…,K i } by applying the EMD technique to it.The ith LR observation signal can then be expressed as Let z j (k) be the instant of the jth zero-crossing in the kth IMF.For isolated IMF samples that are located in the closed interval z j , it is possible to decide whether this interval is noise-or signaldominant by considering the single extrema c i (k) [p j (k) ] in this interval.Note that p j (k) denotes the instant of the extrema for the kth IMF.If the signal is weak as compared with the noise, the absolute value of this extrema will lie below the predefined threshold.Otherwise, in the presence of strong signal, it can be expected that the absolute value of this extrema will exceed the threshold [17].
Let T i (k) be the threshold value used for denoising the kth IMF of the ith LR signal.Application of the interval thresholding to the kth IMF of the ith LR signal is performed as follows [17]: for j = 1,2,… ,N (k) , where N (k) indicates the number of zerocrossings in the kth IMF.In (2), c i (k) [z j (k) ] denotes the samples of the kth IMF belonging to the interval z j (k) .The threshold value required in (2) is selected as [17] ) where C is a multiplication constant and M (i) indicates the sample number of the ith LR observation signal.The IMF energies denoted by {E i (k) } in (3) can be computed as in the following [17]: corresponds to the energy of the first IMF and can be defined as the variance estimate of the AWGN, as well.The parameters β and ρ in (4) depend on the number of sifting iterations performed as a part of the EMD algorithm.Thus, at the end of the EMD-IT based denoising procedure, each noise reduced LR signal can be obtained by It appears from (6) that, after the EMD-IT algorithm [17], the noise reduced version of the ith LR signal is acquired by extracting the first M 1 1 low-order noisy IMFs from (1) and by replacing the high-order less noisy IMFs between M 1 and M 2 in (1) with their thresholded counterparts.
In the second signal processing unit of the proposed method, each denoised LR signal produced by ( 6) is firstly upsampled by L i = L and then each resulting signal is time shifted with z i in order to get estimates of zero interpolated HR signals d ĩ[n] for i = 0,1,…, L 1.
In the third signal processing unit, the nonzero valued samples of the zero interpolated HR signals are combined in terms of the following rule: where n L stands for the common residue of n modulo L.
Finally, in the last stage of the proposed method, the HR signal obtained by ( 7) is applied to a median filter for eliminating the ripples appearing on it.Eventually, an estimate of the desired HR signal, ] [ ˆn d , is acquired.A combination of all signal processing units constituting the proposed method for achieving an estimate of the desired HR signal d[n] from its measured L noisy LR signals is shown by the block diagram in Fig. 3.

Simulation Results
In order to inspect the validity and the performance of the proposed method, we try to reconstruct two different signals from its three LR versions that are corrupted by AWGN.For a reliable evaluation, the results obtained by the proposed method are also compared to those of method [7].Since the LR signals are maximally decimated (downsampled) forms of the considered HR signals, the decimation factor is L i = L = 3 and the number of Wiener filter coefficients required in [7] is taken as P = 8.In the context of the proposed method, the IMFs of each LR signals are obtained by eight sifting iterations and the values of β and ρ constants are respectively fixed to 0.719 and 2.01 during the computation of IMF energies [17].Considering the minimum normalized mean square error (NMSE) related to the original and estimated HR signals, the multiplication constant C in (3) is determined via grid searching in the closed interval [0.1, 1.4] with 0.1 increments such that it meets the minimum NMSE.In order to avoid the ripple effect on the samples of signals reconstructed by the proposed method, a median filter with a window size of 3 is used.Furthermore, the simulation results provided by the method of [7] and the proposed method are compared according to the squared errors e 2 [n] and the NMSE criteria defined by where e[n] expresses the estimation error equaling to the difference between desired {d[n]} and reconstructed HR signals { ] [ ˆn d }.Note that, since the first PL 1 samples of the desired signal cannot be reconstructed with the method of [7], for making one-by-one comparison between the proposed method and the method of [7], the NMSE values are calculated over the last N PL 1 samples.On the other hand, all samples of the desired signal can be estimated by the proposed signal reconstruction method.The desired signal d[n] and its three noisy LR forms are plotted in Fig. 4.
As a part of the EMD-IT based denoising procedure [17] used in the proposed method, the multiplication  constant C in (3) is determined as 0.3 by grid searching and the values of M 1 and M 2 in ( 6) are taken as 6 and K i  2, respectively.The signal reconstruction results and the squared errors provided by the compared two methods are depicted in Fig. 5.It is observed from these figures that the proposed method exhibits better performance than that of method [7].The NMSE values calculated for the proposed method and the method of [7] are 0.031 and 0.109, respectively, support this thought as well.On the other hand, it is worthwhile to note that performances of the EMD-based methods are influenced by the frequency of signal, f 0 , and the sampling frequency, f s (the number of signal samples).
Further details related to this fact can be found in [16].Squared errors for the method of [7] (red dotted line) and the proposed method (blue dashed line).Example 2: This example considers the reconstruction of an audio signal recording that consists of five spoken words 'hello' at certain intervals.This signal is 6.9 seconds long and sampled at 44.1 kHz.In line with the reconstruction purpose, three noisy LR versions of this audio signal with SNR = 0 dB are synthetically produced by the model given in Fig. 1.The original audio signal and its three LR (down-sampled) waveforms are shown in Fig. 6.For this example, on the point of using in the EMD-IT based denoising procedure, the multiplication constant C required in (3) is determined as 0.4 by grid searching and the values of M 1 and M 2 in (5) are taken as 2 and K i -2, respectively.
Using the LR signals shown in Fig. 6b along with the proposed method and the method of [7] lead to the recon-  struction results and the squared errors depicted in Fig. 7.It appears from this figure that, all in all, the proposed method introduces better performance than that of method [7].The respective NMSE values calculated for the proposed method and the method of [7] are 0.156 and 0.223, confirm this claim as well.In the NMSE sense, the reason of providing the proposed method better result than that of the method of [7] can be explained in the following manner.
When we compare Fig. 7a to Fig. 7c, it can clearly be seen that, at the out of the time intervals corresponding to   'hello' sayings (idle channel case), the proposed method reconstructs the audio signal with less noise.The NMSE values presented in Tab. 1 reflect this case.Note that the numerical values tabulated in Tab. 1 are obtained by detecting the time intervals related to in and out of 'hello' sayings (speeches) in the original and reconstructed audio signals.In order to fulfill this detection, a rectangular window having a size of 100 samples is moved on the original signal.After each window movement, variance of 100 samples is calculated.The samples whose variance values are greater than 10 4 are considered as speeches.
in out Method of [7] 0.140 34.607 Proposed method 0.144 5.021 Tab. 1. Computed NMSE values that correspond to the time intervals in and out of "hello" sayings.

Conclusions
A new method for the reconstruction of an HR signal from a set of its available noisy LR observations is proposed.Unlike popular methods that are based on the principle of ME or Wiener filter theory; the proposed method does not need of knowing any second-order statistics related to the desired HR signal and its LR versions.To know the value of down-sampling rate alone is enough for fulfilling the signal reconstruction by using the proposed method.Furthermore, it is demonstrated by several simulations that satisfactory performance has been provided by the proposed method.Thus, the method introduced in this paper can be an option to available methods for HR signal reconstruction.Applying this method to image reconstruction will be a subject of our further direction.
On the other hand, throughout this paper, LR signals {x i [m i ]} are assumed as maximally decimated versions of d[n], i.e.L i = L, m i = m, and M = L. i
HR signal reconstruction based on the proposed method.

Fig. 3 .
Fig. 3. Block diagram representation of the proposed multirate HR signal reconstruction method.

Example 1 :
In this example, we try to reconstruct the HR triangular waveform signal d[n] with N = 4000 samples from its three noisy LR versions obtained by the model shown in Fig. 1, where zero-mean WGN is added on d[n] to make SNR = 5 dB.The following Matlab code with signal frequency f 0 = 4 Hz and sampling frequency f s = 4000 Hz is used to generate signal d[n]: n = 0:1:fs; d = sawtooth(2*pi*f0*(n/fs+0.625),0.5); Three noisy LR versions of d[n].
by the method of[7] (red dotted line) and the proposed method (blue dashed line).

Fig. 5 .
Fig. 5. Example 1: Reconstruction of the triangular waveform signal shown in Fig. 4a from its three noisy LR versions shown in Fig. 4b.
Three noisy LR versions of d[n] in (a).

Fig. 6 . 2 :
Fig. 6.Example 2:The waveform of the audio signal and its three noisy LR forms.
Reconstruction result obtained by the method of[7].Squared error resulted from using the method of[7].Reconstruction result obtained by the proposed method.
Squared error resulted from using the proposed method.

Fig. 7 .
Fig. 7. Example 2: Reconstruction of the HR audio signal shown in Fig. 6a from its three noisy LR versions shown in Fig. 6b.