Method for recovering original signal in reduced complexity DS-CDMA system

Abstract

Disclosed is a method for recovering an original signal in a DS-CDMA system based on complexity reduction. In such a method, first, a partial derivative for rMLL is calculated by using a partial derivative generated by applying fast Fourier transform (FFT) to a reduced negative marginal log likelihood (rMLL) obtained by applying a law of log determinant to a Gaussian process regression (GPR) scheme used for the multi-user detection and thereafter, integrating stochastic gradient descent (SGD). Thereafter, the rMLL is calculated by using the partial derivative for the rMLL and a hyper-parameter is updated to a convergence point until an error gap is converged by repeated calculation of the rMLL. Next, a kernel function used for estimating a matched filter is calculated by using the hyper-parameter estimated through the convergence and the original signal for each of multi-users is recovered by using the kernel function.

Description

TECHNICAL FIELD

The present invention relates to the field of communication, machine learning, and optimization, and more particularly, to a method for recovering an original signal in a direct sequence-code division multiple access (DS-CDMA) system based on complexity reduction using Gaussian process regression.

BACKGROUND ART

A direct sequence-code division multiple access (DS-CDMA) system distinguishes a user by a signal of a frequency band. Unfortunately, inter-signal interference occurs even among a small number of users and is recognized as multiple access inference (MAI). This noise problem causes a critical increase in bit error rate (BER) under a near/far effect, which is a sensitive issue in the DS-CDMA system.

To mitigate such a problem, a multiuser detection (MUD) technique has been applied to eliminate the interference. In such a technique, a known optimization solving means for the MUD may be acquired through minimizing the means square error (MMSE). Nonetheless, a lot of computational resources and training effort are required for performing the calculation. Obviously, such a method is not suitable for implementation in most communication devices. In order to solve such a problem, many approaches including multilayer perceptron, a support vector machine, a wavelet neuron network, and Gaussian process regression (GPR) are proposed. Of the machine learning approaches, the GPR is considered the most promising tool in terms of flexibility and accuracy.

In fact, a Gaussian process is widely used for prediction and classification in many research areas such as data communications, networking, and signal processing. Rather than determining parameters of a model from a scratch, the Gaussian process can help to adopt the parameters to represent an actual underlying function. As such, the Gaussian process is a suitable choice for noise, corrupted or erroneous data. However, such a method has a disadvantage of high complexity. In a standard implementation, the GPR requires complexity of O(n³) for computation and O(n²) for storage when computing n training points in a dataset. Even in the application of a sparse spectrum Gaussian process, if m is the number of basic functions, the complexity is still O(mn²) for computation and O(nm) for storage.

Accordingly, there is a need for another scheme capable of reducing the complexity of the GPR.

DISCLOSURE

Technical Problem

The present invention has been made in an effort to provide a method for recovering an original signal in a DS-CDMA system which can reduce complexity of GPR.

Technical Solution

An exemplary embodiment of the present invention provides a method for recovering an original signal through multi-user detection in a synchronous mobile communication system, including: calculating a partial derivative for rMLL by using a partial derivative generated by applying fast Fourier transform (FFT) to a reduced negative marginal log likelihood (rMLL) obtained by applying a law of log determinant to a Gaussian process regression (GPR) scheme used for the multi-user detection and thereafter, integrating stochastic gradient descent (SGD); calculating the rMLL by using the partial derivative for the rMLL; updating a hyper-parameter to a convergence point until an error gap is to converged by repeated calculation of the rMLL; calculating a kernel function used for estimating a matched filter by using the hyper-parameter estimated through the convergence; and recovering the original signal for each of multi-users by using the kernel function. Where the partial derivative may follow an equation below.

$\frac{\partial }{\partial {\theta }_1}{ℱ}_{\mathrm{rMLL}}=-\frac{\sqrt{\frac{2}{\pi }}{\tilde{x}}_k^2\left(i\right)\exp \left(2{\pi }^2{\theta }_2^2{\omega }^2\right)}{{\theta }_2{\theta }_1^3}$

Where

${ℱ}_{\mathrm{rMLL}}\left(\theta \right)=\frac{1}{2n}\sum _k\frac{{\tilde{x}}_k^2\left(i\right)}{{ℱ}_{\mathrm{SE}}\left({\omega }_k\right)}$ and represents Fourier transform for the rMLL, θ₁ represents an output-scale size, θ₂ represents a time-scale of y from one moment up to a next moment, ω is a frequency representation of a received signal y in a frequency domain, a tilde symbol represents Fourier transform of an original signal x(i) in {tilde over (x)}(i), and _SE(ω)=θ₂θ₁²√{square root over (2π)}exp(−2π²ω²θ₂²) and represents Fourier transform for the kernel function. Further, the error gap (RMSE) by repeated computation of the rMLL may be evaluated through a relationship equation below.

$\mathrm{RMSE}=\sqrt{\frac{\sum _{j=1}^n{ℱ}_j^{\left(k\right)}-{ℱ}_j^{\left(k-1\right)}}{n}}$

Where _j^(k) and _j^(k-1) represent a value of the rMLL at a target location j after k^th and (k−1)^th repetitions, respectively, and n represents the number of repetitions. In addition, the update may be performed by using the relationship equation below.

${\theta }_2^{\left(k\right)}\leftarrow {\theta }_2^{\left(k-1\right)}+\alpha \left(k\right)\frac{\partial }{\partial {\theta }_2^{\left(k-1\right)}}{ℱ}_{\mathrm{rMLL}}$ ${\theta }_1^{\left(k\right)}\leftarrow {\theta }_1^{\left(k-1\right)}+\alpha \left(k\right)\frac{\partial }{\partial {\theta }_1^{\left(k-1\right)}}{ℱ}_{\mathrm{rMLL}}$

Where α(k)=1/(k+1) represents a Robbins-Monroe decay function for k^th repetitions. The kernel function follows the relationship equation below.

$k\left(y_i,y_j\right)={\theta }_1^2\exp \left(-\frac{{y_i-y_j}^2}{2{\theta }_2^2}\right)$

Further, applying the law of log determinant to the Gaussian process regression (GPR) scheme may be convergence to a constant using a law of log determinant log |P| while minimizing the following relationship equation which is a negative log likelihood for estimating a hyper-parameter set constituted by

${\theta }_1\mathrm{and}{\theta }_2-\log p\left(x\left(i\right)|y,\theta \right)=\frac{1}{2}{x}^T\left(i\right)P^{-1}x\left(i\right)+\frac{1}{2}\log P+\frac{n}{2}\log \left(2\pi \right)$

Where p(x(i)|y,w)=(x(i)|Φ^Tw,σ_noise²I) and represents a joint Gaussian distribution, the matched filter w is a random variable and follows the Gaussian distribution like p(w)=(w|0,σ_w²I) which is a probability of w, σ_noise and σ_w represent standard deviations of noise and the matched filter, respectively, I_n represents a unit matrix having a size of n, Φ=[ϕ(y₁), ϕ(y₂)], . . . , ϕ(y_n) represents a vector of non-linear mapping to a higher dimensional space for a received signal x(i), and P represents a covariance matrix.

Advantageous Effects

According to the present invention, the complexity of the GPR can be reduced to become O(n log n) for computation and O(n) for storage when n is the number of training points. Consequently, through the improvement, by comparing the Gaussian process in the related art and MMSE estimation, regression can be greatly accelerated so as to maintain preferable BER.

DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram illustrating a synchronous DS-CDMA system model to which the present invention is applied.

FIG. 2 is a flowchart of a method for recovering an original signal in a DS-CDMA system according to an exemplary embodiment of the present invention.

MODE FOR INVENTION

In the following detailed description, only certain exemplary embodiments of the present invention have been shown and described, simply by way of illustration. As those skilled in the art would realize, the described embodiments may be modified in various different ways, all without departing from the spirit or scope of the present invention. Accordingly, the drawings and description are to be regarded as illustrative in nature and not restrictive. Like reference numerals designate like elements throughout the specification.

Throughout the specification, unless explicitly described to the contrary, the word “comprise” and variations such as “comprises” or “comprising”, will be understood to imply the inclusion of stated elements but not the exclusion of any other elements. In addition, the terms “-er”, “-or” and “module” described in the specification mean units for processing at least one function and operation and can be implemented by hardware components or software components and combinations thereof.

First, an exemplary embodiment of the present invention may be applied to a synchronous DS-CDMA system illustrated in FIG. 1. In such a system, all users transmit symbols at the same symbol rate. In such a system, symbol k is transmitted at moment j. Each user's symbol s_j(l) is multiplied by a spreading code c_l. In this case, the spreading code c_l is a sequence of n_S pseudorandom binary values considered as a chip. Resulting signals are amplified by different gains a_l through a multiplier 11 of a transmitter 10 and then summed and transmitted in various forms via a summer 12.

The transmitted signal is received via a channel response h(z) 20 in a chip period. In this case, additive white Gaussian noise (AWGN) is added during transmission and is indicated through a summer 12 in a system model in FIG. 1.

Meanwhile, a multi-user detector (MUD) 31 of a receiver 30 restores a transmission bit for a specific user by using the received signal, that is, a chip signal output through the summer 12, and outputs the restored transmission bit as

A. Assumption

An input data set having n training points is referred to as D={x_t,y_t|_{t=1, . . . , n}}. Herein, {x_t∈^d} represents a vector of an original signal and y_t represents a column vector of the received signal at a time step t. A relationship of the vectors is shown in [Equation 1].y_t=SAx_t+n_t  [Equation 1]

Where S represents a U×V matrix (each column of such a matrix containing a spreading code for each user). A represents a V×V diagonal matrix and contains a size for a user. The size represents a fading degree of the signal transmitted over a channel (the fading degree indicates how far the user is from the receiver). Last, n_t represents the additive white Gaussian noise (AWGN) added to the channel with time progress.

In the receiver, an original signal x_t(i) of an i-th i^th user needs to be restored as shown in [Equation 2].{circumflex over (x)}_t(i)=sgn(w_i^Ty_t)  [Equation 2]

Where w_i represents is a matched filter for the i^th user . Event though w_i is an optimal non-linear vector, such a vector may be estimated as shown in [Equation 3] by using an MMSE method.

$\begin{array}{cc}{w}_i^{*}=\underset{w_i}{\mathrm{argmin}}𝔼\left[{\left(x_t\left(i\right)-w_i^T{y}_t\right)}^2\right]=C_{\mathrm{yy}}^{-1}{C}_{\mathrm{yx}}& \left[\mathrm{Equation}3\right]\end{array}$

Where C_yy=[y_ty_t^T] represents an autocorrelation of the received vector and C_yx=[y_tx_t(i)] represents a cross-correlation between the received vector and the original. [Equation 3] is known as distributed MMSE and may be solved without finding spreading sequences of other users. However, the problem of the solving means is that a huge training data set for an inverse matrix is required and high computational complexity is also required.

B. Deriving MUD of Gaussian Process

Φ=[ϕ(y₁), ϕ(y₂)], . . . , ϕ(y_n) represents the vector of non-linear mapping to a higher dimensional space of the received signal and ϕ(⋅) represents a corresponding mapping function. When an original signal vector x(i)=[x₁(i), x₂(i), . . . , x_n(i)]^T is applied to the received signal vector y=[y₁, y₂, . . . , y_n], the original signal vector to the received signal vector has a joint Gaussian distribution as shown in [Equation 4].p(x(i)|y,w)=(x(i)|Φ^Tw,σ_noise²I)  [Equation 4]

Where the matched filter w is actually a random variable. A probability of w follows a Gaussian distribution like p(w)=(w|0,σ_w²I) and σ_noise and σ_w represent standard deviations of noise and the matched filter, respectively and I_n represents a unit matrix having a size of n. By applying a Bayes' law to [Equation 4], a posterior distribution of w is computed as shown in [Equation 5].

$\begin{array}{cc}p\left(w|x\left(i\right),y\right)=\frac{p\left(w\right)p\left(x\left(i\right)|w,y\right)}{p\left(x\left(i\right)|y\right)}& \left[\mathrm{Equation}5\right]\end{array}$

Theoretically, [Equation 3] may be transformed into a nonlinear form using maximum a posterior (MAP) estimation of the random variable w for [Equation 5]. The transformation is expressed as shown in [Equation 6].

$\begin{array}{cc}{w}^{*}=\underset{w}{\mathrm{argmin}}\left\{{x\left(i\right)-{\Phi }^{T}w}^2+\lambda {w}^2\right\}& \left[\mathrm{Equation}6\right]\end{array}$

Where λ=(σ_noise/σ_w)². A equation term λ∥w∥² is integrated into a MAP as a regularizer to skip an over-fitting problem. By discovering w*, the estimation of an original signal {circumflex over (x)}(i) may be obtained as shown in [Equation 7].{circumflex over (x)}(i)=k^TP⁻¹x(i)  [Equation 7]

Where in the case of k(y_i,y_j)=(ϕ(y_i)^Tϕ(y_j)), k=[k(y,y₁), k(y,y₂), . . . , k(y,y_n)] represents a kernel function of the non-linear transformation and P=K+σ_noise²I (where K represents a covariance matrix which is K_ij=k(y_i,y_j))). Due to a priority of a processing speed, a square exponential kernel function is adopted to compute the estimate of the original signal. Such a kernel function is shown in [Equation 8].

$\begin{array}{cc}k\left(y_i,y_j\right)={\theta }_1^2\exp \left(-\frac{{y_i-y_j}^2}{2{\theta }_2^2}\right)& \left[\mathrm{Equation}8\right]\end{array}$

Where θ₁ represents an output-scale size and θ₂ represents a time-scale of y from one moment to a next moment. A set θ={θ₁,θ₂} is known as a set of hyper-parameters. Next, the estimation of θ may be computed as shown in [Equation 9].

$\begin{array}{cc}{\theta }^{*}=\underset{\theta }{\mathrm{argmin}}p\left(x\left(i\right)|y,\theta \right)& \left[\mathrm{Equation}9\right]\end{array}$

In order to solve [Equation 9], O(n³) for the computational complexity and O(n²) for storage when the size of the storage is n are assumed. Most of the complexity is caused by computing the inverse matrix and a log matrix. Obviously, such a disadvantage is a burden on the DS-CDMA system. In order to solve such a problem, a complexity reduction method is proposed to significantly accelerate a computation process. By using such a method, the computation and storage complexity corresponds to O(n log n) and O(n), respectively.

C. Reduction of Complexity

The proposed complexity reduction method is a combination of three techniques and three techniques are fast Fourier transform (FFT), a law of log determinant and application of stochastic gradient descent (SGD). By definition, the hyper-parameter set θ* may be estimated as shown in [Equation 10] by minimizing negative log likelihood −log p(x(i)|y,θ).

$\begin{array}{cc}-\log p\left(x\left(i\right)|y,\theta \right)=\frac{1}{2}{x}^T\left(i\right)P^{-1}x\left(i\right)+\frac{1}{2}\log P+\frac{n}{2}\log \left(2\pi \right)& \left[\mathrm{Equation}10\right]\end{array}$

Due to the high complexity for computing the inverse matrix P⁻¹, it is necessary to develop an approximation method. It may be a better solving means to approximately minimize an upper bound of the equation term in terms of the approximation, instead of trying to minimize the negative log likelihood. Analytically, in [Equation 10], a better computation focuses on two terms, that is, a data-fit term expressed as x^T(i)P⁻¹x(i) and the log determinant log |P|. First, a simplification derivation needs to be applied to the terms to reduce the equation. To this end, the law of the log determinant is used to calculate the log determinant {circumflex over (P)} of the empirical covariance matrix {circumflex over (P)}, which simplifies [Equation 10] to [Equation 11].

$\begin{array}{cc}-\log p\left(x\left(i\right)|y,\theta \right)=\frac{1}{2}{x}^T\left(i\right)P^{-1}x\left(i\right)+\frac{1}{2}\log \hat{P}+\frac{n}{2}\log \left(2\pi \right)& \left[\mathrm{Equation}11\right]\end{array}$

Where {circumflex over (P)} is calculated according to [Equation 12] based on the empirical covariance matrix and a constant τ.

$\begin{array}{cc}\hat{P}=\log \hat{P}-\tau w.r.t\hat{P}=\frac{1}{n}\sum _{k=1}^{n+1}{\left[x_k\left(i\right)-\stackrel{\_}{x}\left(i\right)\right]\left[x_k\left(i\right)-\stackrel{\_}{x}\left(i\right)\right]}^T\tau =\gamma \left(\frac{n}{2}\right)-\log \left(\frac{n}{2}\right)& \left[\mathrm{Equation}12\right]\end{array}$

Where γ(⋅) represents a Digamma function. After multiple recalculation, the term {circumflex over (P)} a converges to the constant according to a central limit theorem. Minimizing the negative log marginal likelihood after the passage of time due to such convergence may involve only minimizing the reduced negative marginal log likelihood (rMLL) as shown in [Equation 13].

$\begin{array}{cc}-\log p\left(x\left(i\right)|y,\theta \right)=\frac{1}{2}{x}^T\left(i\right)P^{-1}x\left(i\right)& \left[\mathrm{Equation}13\right]\end{array}$

A gap of the approximation step with respect to the actual solving means may be very small and may be measured using a mean square error of [Equation 14].

$\begin{array}{cc}{𝔼\left(\hat{P}-\log P\right)}^2=-2\log \left(1-\frac{1}{n}\right)& \left[\mathrm{Equation}14\right]\end{array}$

Nevertheless, solving the inverse matrix P in [Equation 13] still costs computationally more. Therefore, another method is required in order to achieve the object. Since the covariance matrix P has a positive-definite value, it is possible to perform the transformation using FFT. Such a technique is used for taking the calculation from a spatial-temporal domain to a frequency domain. It is worth mentioning that the cost of the FFT is only O(n log n). Obviously, the cost is much better than the method in the related art. Hereinafter, such transformation will be described in detail.

First, the squared exponent kernel k(y_i,y_j) in Equation (8) needs to be rewritten into a Fourier transform representation as shown in [Equation 15]._SE(ω)=θ₂θ₁²√{square root over (2π)}exp(−2π²ω²θ₂²)  [Equation 15]

Where ω is a frequency representation of the received signal y in the frequency domain. It is assumed that Ψ generates a function {tilde over (P)}=P⁻¹. Under the frequency domain, the Parseval theorem is applied to derive the Fourier transform for [Equation 13].

$\begin{array}{cc}{ℱ}_{\mathrm{rMLL}}\left(\theta \right)=ℱ\left(-\log p\left(x\left(i\right)|y,\theta \right)\right)=\frac{1}{2n}{x}^T\left(i\right)\Psi *x_o\left(i\right)\hspace{1em}& \left[\mathrm{Equation}16\right]\end{array}$

Where a tilde symbol represents the Fourier transform of x(i) in {tilde over (x)}(i) and x_o(i) represents a data vector in a periodic domain. In a next step, a convolution theorem is continuously applied to a limit Ψ_SE≡1 and last Fourier transform of rMLL may be expressed as shown in [Equation 17].

$\begin{array}{cc}{ℱ}_{\mathrm{rMLL}}\left(\theta \right)=\frac{1}{2n}\sum _n{\Psi }_k*{\tilde{x}}_k^2\left(i\right)=\frac{1}{2n}\sum _k\frac{{\hat{x}}_k^2\left(i\right)}{{ℱ}_{\mathrm{SE}}\left({\omega }_k\right)}& \left[\mathrm{Equation}17\right]\end{array}$

In the form of [Equation 17], θ which is the set of hyper-parameters may be estimated using a slope-based technique. In this case, the stochastic gradient descent (SGD) is selected because of a characteristic that the SGD is less sensitive to fast convergence and local minima. In order to integrate the SGD, a partial derivative of [Equation 17] is required for each hyper-parameter. Such an equation is given by [Equation 18].

$\begin{array}{cc}\frac{\partial }{\partial {\theta }_2}{ℱ}_{\mathrm{rMLL}}={\tilde{x}}_k^2\left(i\right)\exp \left(2{\pi }^2{\theta }_2^2{\omega }^2\right)\left(\frac{2\sqrt{2}{\pi }^{3/2}{\omega }^2}{{\theta }_1^2}-\frac{1}{\sqrt{2\pi }{\theta }_2^2{\theta }^{2^1}}\right)\frac{\partial }{\partial {\theta }_1}{ℱ}_{\mathrm{rMLL}}=-\frac{\sqrt{\frac{2}{\pi }}{\tilde{x}}_k^2\left(i\right)\exp \left(2{\pi }^2{\theta }_2^2{\omega }^2\right)}{{\theta }_2{\theta }_1^3}& \left[\mathrm{Equation}18\right]\end{array}$

Subsequently, an update procedure is applied to update the hyper-parameters to corresponding convergence points.

Such a process is expressed as shown in [Equation 19].

$\begin{array}{cc}{\theta }_2^{\left(k\right)}\leftarrow {\theta }_2^{\left(k-1\right)}+\alpha \left(k\right)\frac{\partial }{\partial {\theta }_2^{\left(k-1\right)}}{ℱ}_{\mathrm{rMLL}}{\theta }_1^{\left(k\right)}\leftarrow {\theta }_1^{\left(k-1\right)}+\alpha \left(k\right)\frac{\partial }{\partial {\theta }_1^{\left(k-1\right)}}{ℱ}_{\mathrm{rMLL}}& \left[\mathrm{Equation}19\right]\end{array}$

Where α(k)=1/(k+1) represents a Robbins-Monroe decay function for k^th repetition. Such a function is chosen instead of exact line search or backtracking line search or backtracking line search primarily due to a performance problem. Further, in order to control the number of update repetitions, an error function is defined based on a root mean square error (RMSE) method to measure the convergence. It should be noted that the RMSE method is stricter than the means square error (MSE) method which is generally used. By using such an error function, an error gap between a current repetition value and a previous value may be evaluated as shown in [Equation 20].

$\begin{array}{cc}\mathrm{RMSE}=\sqrt{\frac{{\sum }_{j=1}^n{ℱ}_j^{\left(k\right)}-{ℱ}_j^{\left(k-1\right)}}{n}}& \left[\mathrm{Equation}20\right]\end{array}$

Where _j^(k) and _j^(k-1) represent a value of rMLL at a target location j after k^th and (k−1)^th repetitions, respectively. In the proposed method, an RMSE threshold is limited to 10⁻¹¹ which yields a solving means that is close to an actual value. Obviously, the computation may be performed without the inverse matrix. By the end of the method, a required set θ* of hyper-parameters is obtained. In addition, due to the fact that the method does not need to maintain a historical covariance matrix, which is primarily used to inverse the covariance matrix, the computational complexity is only O(n log n) and the storage complexity is O(n).

Hereinafter, a method for recovering the original signal in the DS-CDMA system using various equations derived through the above process will be described.

FIG. 2 is a flowchart of a method for recovering an original signal in a DS-CDMA system according to an exemplary embodiment of the present invention.

Referring to FIG. 2, first, a receiver 30 collects a reception chip received through a channel 20 from a transmitter 10 (S100).

Thereafter, the hyper-parameter to be estimated is initialized in the exemplary embodiment of the present invention (S110).

Next, a partial derivative of rMLL is calculated using a finally derived [Equation 18] by applying fast Fourier transform, a law of log determinant, and stochastic gradient descent (SGD) to GPR according to the exemplary embodiment of the present invention (S120).

Subsequently, the rMLL is calculated using the calculated partial derivative result (S130).

Then, the RMSE is evaluated using [Equation 20] that calculates the error gap between the current repetition value and the previous value (S140).

After determining whether an RMSE evaluation result converges to a preset threshold value (S150), if the RMSE evaluation result does not converge to the preset threshold value, the update is performed using [Equation 19] to update the RMSE evaluation result to the convergence point corresponding to the hyper-parameter (S160).

However, when it is determined in step S150 that the RMSE evaluation result converges to the preset threshold value, the kernel function is calculated using [Equation 8] (S170) and the original signal for a specific user is finally restored using the calculated kernel function (S180).

As described above, in the exemplary embodiment of the present invention, by applying the fast Fourier transform, the law of the log determinant, and the stochastic gradient descent (SGD) to the existing GPR, the regression may be significantly accelerated so that the computational complexity and storage complexity in the case of using the existing GPR are reduced and better BER is maintained.

While this invention has been described in connection with what is presently considered to be practical example embodiments, it is to be understood that the invention is not limited to the disclosed embodiments, but, on the contrary, is intended to cover various modifications and equivalent arrangements included within the spirit and scope of the appended claims.

Claims

1. A method for recovering an original signal through multi-user detection in a synchronous mobile communication system, comprising: calculating a partial derivative of a reduced negative marginal log likelihood (rMLL), thereafter, integrating stochastic gradient descent (SGD) by applying a partial derivative to a fast Fourier transform (FFT) of the rMLL, the rMLL is obtained by applying a law of log determinant to a Gaussian process regression (GPR) scheme used for the multi-user detection;calculating the rMLL by using the partial derivative for the rMLL;updating a hyper-parameter to a convergence point until an error gap is converged by repeated calculation of the rMLL;calculating a kernel function used for estimating a matched filter by using the hyper-parameter estimated through the convergence point; andrecovering the original signal for each of multi-users by using the kernel function.

2. The method of claim 1, wherein: the partial derivative follows an relationship equation below,

3. The method of claim 2, wherein: the error gap (RMSE) by repeated computation of the rMLL is evaluated through a relationship equation below,

4. The method of claim 3, wherein: the update is performed by using the relationship equation below,

5. The method of claim 4, wherein: the kernel function follows the relationship equation below,

6. The method of claim 4, wherein: applying the law of log determinant to the Gaussian process regression (GPR) scheme is convergence to a constant using a law of log determinant log |P| while minimizing the following relationship equation which is a negative log likelihood for estimating a hyper-parameter set (θ) constituted by θ1 and θ2,