Channel estimation method of NOMA (Non-Orthogonal Multiple Access) system based on vector approximate message passing

Abstract

The invention discloses a channel estimation method of an NOMA (Non-Orthogonal Multiple Access) system based on vector approximate message passing, which comprises the following steps of: firstly, based on the sparsity of channel pulse response in a time domain, adopting Bernoulli-Gaussian distribution as prior distribution of the channel pulse response; then obtaining the joint probability of the received signal and the channel pulse response according to the Bayesian theory, and designing a factor graph model of channel estimation according to the joint probability; and finally, iteratively solving the posterior probability of the channel impulse response according to the factor graph and a sum-product algorithm, and updating unknown hyper-parameters in prior distribution through an expectation maximization algorithm in the iteration process, thereby obtaining an optimal channel estimation value. The method has the advantages that accurate channel prior information does not need to be provided, the mean square error and the bit error rate of channel estimation are reduced, and the calculation complexity is low.

Description

Technical Field

[0001] The present invention relates to a channel estimation technology, and in particular to a channel estimation method of a NOMA (Non-Orthogonal Multiple Access) system based on vector approximate message passing (VAMP). Background Art

[0002] Non-Orthogonal Multiple Access (NOMA) technology is an emerging 5G technology that can meet the growing bandwidth requirements of next-generation networks. The core advantage of the NOMA system is that it allows multiple users to communicate with different powers on the same time and frequency resources. This resource sharing advantage improves spectrum efficiency and user capacity, which is of great significance for high-density user scenarios and large-scale connected IoT devices, and can meet the growing communication needs. In the NOMA system, the signal will be subject to various interferences during the propagation process. When the signal reaches the receiving end from the transmitter, the amplitude, phase and frequency of the signal will change greatly, and the role of channel estimation and channel equalization is to restore the signal as accurately as possible.

[0003] Existing channel estimation methods can be divided into two categories. The first category is blind / semi-blind channel estimation methods, which obtain "channel state information" from the structure and statistical information of the received signal, without or with very little training sequence. Its advantage is that it reduces resource overhead, but the computational complexity is high, and problems such as phase ambiguity may occur, which to a certain extent limits its practicality; the second category is reference signal-based channel estimation methods, which need to send pilot signals or training sequences in addition to data symbols at the transmitter. At present, least squares (LS), linear minimum mean square error (LMMSE) and other training sequence-based channel estimation methods are widely used. Generally speaking, the overall performance of reference signal-based channel estimation methods is better than that of blind / semi-blind channel estimation methods.

[0004] Since the wireless channel is sparse, that is, most of the tap energies of the channel impulse response are approximately zero, but the traditional LS and LMMSE do not utilize the sparsity of the wireless channel, which will cause a waste of resources. The most representative sparse channel estimation method is the channel estimation method based on compressed sensing. In the context of compressed sensing, existing work mainly uses sparse Bayesian learning (SBL) theory to estimate sparse channels. The SBL algorithm is relatively stable, but each iteration process involves a matrix inversion process, and the computational complexity is relatively high. Donoho et al. proposed an approximate message passing (AMP) algorithm based on factor graphs that can effectively solve compressed sensing problems. Subsequently, AMP was extended to generalized approximate message passing (GAMP). Some scholars applied GAMP to joint channel estimation and decoding, but GAMP can only converge when the elements in the observation matrix obey independent and identically distributed Gaussian distributions, which is unstable and has poor channel estimation performance. In addition, AMP algorithms require accurate prior information of the quantity to be estimated. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a channel estimation method for a NOMA system based on vector approximate message passing, which does not require accurate channel prior information and can obtain a better channel estimation value in the NOMA system.

[0006] The technical solution adopted by the present invention to solve the above technical problems is: a channel estimation method for a NOMA system based on vector approximate message transmission, characterized by comprising the following steps:

[0007] Step 1: At the transmitting end of the NOMA downlink communication system, the transmitting end sends a time domain signal t to the receiving end. The time domain signal t passes through the channel, where the dimension of t is N×1, N represents the total number of subcarriers in the NOMA system, N=N d +N p , N d Indicates the number of data subcarriers, N p Indicates the number of pilot subcarriers;

[0008] Step 2: At the receiving end of the NOMA downlink communication system, the receiving end receives the time domain signal r, r = Gt + ε; then r is converted from the time domain to the frequency domain through discrete Fourier transform to obtain the frequency domain received signal R, Then, according to R, we get the frequency domain received signal R corresponding to the pilot subcarrier. p ,

[0009] Where G represents an N×N circulant matrix, the first L elements in the first column vector of G correspond to the L elements in the L×1 channel impulse response vector h, and the last NL elements are all 0, L represents the delay spread length of the channel, h=[h0,h1,…,h L-1 ] T , the element with index l in h is h l , l=0,1,...,L-1, the superscript “T” represents the transpose of a matrix or vector, ε represents the additive white Gaussian noise in the time domain, F represents the normalized discrete Fourier transform matrix, the dimension of F is N×N, F * is the Hermitian transform of F. The two binary transmission sequences at the transmitter are encoded by the encoder and then mapped into frequency domain symbols through QAM. The two frequency domain symbols are superimposed into signal u after power allocation. X is a diagonal matrix, X = diag(u), diag() means to find the diagonal matrix, and the diagonal elements of X are u1,u2,…,u N , u1 represents the first element in u, u2 represents the second element in u, u N represents the Nth element in u, F L represents the matrix composed of the first L column vectors of F, w represents the background noise vector in the frequency domain, w=Fε, satisfying CN(w;0,z w -1 I N ), CN is the symbol for complex Gaussian distribution, I N represents the N-order identity matrix, z w -1 is the variance of w, X p represents a matrix consisting of row vectors corresponding to the pilot subcarrier index positions in X, w p represents a vector consisting of row elements in w corresponding to the pilot subcarrier index positions,

[0010] Step 3: According to the sparsity of h, h is assumed to obey an independent and identically distributed Bernoulli-Gaussian distribution, and the prior probability p(h|θ1) about h is obtained. Then according to Bayesian theory, we get p ,h joint probability p(R p ,h|θ), Reintroducing h 1 and h 2 , and let h 1 =h 2 =h, about R p ,h joint probability p(R p ,h|θ) can be rewritten as p ,h 1 ,h 2 The joint probability p(Rp ,h 1 ,h 2 |θ),

[0011] Where θ1=[λ,z h ] T ,λ represents the sparsity of h,λ∈(0,1),z h represents the inverse variance of the non-zero elements in h, p(h l |θ1) represents the relationship between h l The prior probability of δ() is the Dirac function, CN(h l ; 0,z h -1 ) means h l Subject to mean 0 and variance z h -1 The complex Gaussian distribution, θ = [θ1 T ,z w ] T , p(R p |h,θ) means R under the condition that h,θ are known p The likelihood function of Represents R p Subject to mean Φh and variance The complex Gaussian distribution of N p The identity matrix of order, p(h 1 |θ1) represents the relationship between h 1 The prior probability of Represents R p The mean is Φh 2 And the variance is The complex Gaussian distribution of

[0012] Step 4: According to Draw the factor graph of channel estimation, which has two variable nodes h 1 and h 2 、Three factor nodes p(h 1 |θ1) and δ(h 1 -h 2 )and

[0013] Step 5: Based on the factor graph, use the vector approximate message passing algorithm to approximate the estimated value of h. The specific process is as follows:

[0014] Step 5.a1: Let k represent the number of iterations, the initial value of k is 1; let K represent the maximum number of iterations;

[0015] Step 5.a2: During the kth iteration, calculate the factor node δ(h1 -h 2 ) is passed to the variable node h 1 News Then according to p(h 1 |θ1) and Using the Bayesian method, we can get 1 The posterior probability p(h 1 |ω 1,k ,z 1,k ,θ1), Then Substitution And simplify it to get Further get h 1 The posterior expectation β 1,k and the mean variance Finally, we get h 1 The approximate belief b app (h 1 ),

[0016] Among them, CN(h 1 ;ω 1,k ,z 1,k -1 I L ) means h 1 The mean is ω 1,k And the variance is z 1,k -1 I L The complex Gaussian distribution, ω 1,k represents the kth iteration h 1 Expectations, ω 1,1 =0, z 1,k -1 represents the kth iteration h 1 The variance of z 1,1 =1,I L represents the identity matrix of order L, Indicates h 1 The element with index l in Indicates about The prior probability of express Subject to mean 0 and variance z h -1 The complex Gaussian distribution, ρ k (l) represents ρ k The element with index l in ρ k is the introduced intermediate quantity, ω 1,k (l) represents ω 1,k The element with index l in CN(ω 1,k (l); 0, z1,k -1 ) represents ω 1,k (l) has a mean of 0 and a variance of z 1,k -1 The complex Gaussian distribution, CN(ω 1,k (l); 0, (z 1,k -1 +z h -1 )) represents ω 1,k (l) has a mean of 0 and a variance of z 1,k -1 +z h -1 The complex Gaussian distribution of k (l) indicates a k The element with index l in a k is the introduced intermediate quantity, v k is the intermediate quantity introduced, express The mean is a k (l) and the variance is v k The complex Gaussian distribution of 1,k is the intermediate quantity introduced, Indicates h 1 The mean is β 1,k And the variance is The complex Gaussian distribution of

[0017] Step 5.a3: During the kth iteration, calculate the variable node h 1 Passed to the factor node δ(h 1 -h 2 ) According to the sum-product algorithm, we can get

[0018] Among them, ω 2,k represents the kth iteration h 2 expectations, z 2,k -1 represents the kth iteration h 2 The variance of CN(h 1 ;ω 2,k ,z 2,k -1 I L ) means h 1 The mean is ω 2,k And the variance is z 2,k -1 I L The complex Gaussian distribution of

[0019] Step 5.a4: During the kth iteration, the factor node δ(h 1 -h 2 ) is passed to the variable node h 2 News No change, obtained Then according to and Using the Bayesian method, we can get 2 The posterior probability And simplify to get h 2 The posterior expectation β 2,k and the mean variance β 2,k =(z w Φ H Φ+z 2,k I L ) -1 (z w Φ H R p +z 2,k ω 2,k ), Finally, we get h 2 The approximate belief b app (h 2 ),

[0020] Among them, CN(h 2 ;ω 2,k ,z 2,k -1 I L ) means h 2 The mean is ω 2,k And the variance is z 2,k -1 I L The superscript "H" indicates the conjugate transpose of the matrix. 2,k is the intermediate quantity introduced, Tr[] represents matrix trace, Indicates h 2 The mean is β 2,k And the variance is The complex Gaussian distribution of

[0021] Step 5.a5: Determine whether k is less than K. If so, execute step 5.a6; otherwise, according to the minimum mean square error estimation criterion, obtain the minimum mean square error estimate of h. in, represents the estimated value of h, represents the true value of h, || ||2 represents the l2 norm, Indicates that in R p Under known conditions The posterior probability of Indicates that in R p Under known conditions expectations, Expressing expectation; based on β 1,K As an estimate of h

[0022] Step 5.a6: Calculate variable node h 2 Passed to the factor node δ(h 1 -h 2 ) According to the sum-product algorithm, we can get Among them, CN(h 2 ;ω 1,k+1 ,z 1,k+1 -1 I L ) means h 2 The mean is ω 1,k+1 And the variance is z 1,k+1 -1 I L The complex Gaussian distribution of Then for β 2,k With φ 2,k Perform singular value decomposition on Φ in the expression of , and get Φ = USV H , further obtained Among them, S is a diagonal matrix, U and V are both unitary matrices, Y represents the rank of Φ, and d k is the intermediate quantity introduced, d k =z w diag(z w S H S+z 2,k 1 L ) -1 S H S, 1 L represents an L×L dimensional all-one matrix, d k (l) indicates d k The element with index l in ; then let k = k + 1, return to step 5.a2 to continue iteration, where the "=" in k = k + 1 is an assignment symbol;

[0023] The expectation maximization algorithm is used simultaneously to calculate the estimated value of θ. The specific process is as follows:

[0024] Step 5.b1: Let k represent the number of iterations, the initial value of k is 1; let K represent the maximum number of iterations;

[0025] Step 5.b2: During the kth iteration, perform E-step: 1,k-1 Under known conditions, calculate the joint probability p(R p ,h|θ) in the posterior probability p(h|ω 1,k ,z 1,k ,θ 1,k-1 ) In z w,k-1 Under known conditions, calculate the joint probability p(R p ,h|θ) in the posterior probability p(h|ω 2,k ,z 2,k ,z w,k-1 ) When k = 1, θ 1,k-1 represents the initial value of θ1. When k>1, θ 1,k-1 represents the estimated value of θ1 at the k-1th iteration, and p(h 1 |ω 1,k ,z 1,k ,θ1) in h 1 Substituting h for p(h|ω 1,k ,z 1,k ,θ1), and p(h|ω 1,k ,z 1,k ,θ1) is replaced by θ 1,k-1 We get p(h|ω 1,k ,z 1,k ,θ 1,k-1 ), ln() represents the natural logarithm function, p(R p |h,z w ) means in h,z w Under known conditions, R p The likelihood function of exp() represents an exponential function with the natural constant e as the base. When k = 1, z w,k-1 Indicates z w The initial value of z is w,k-1 represents the k-1th iteration z w The estimated value of p(h 2 |ω 2,k ,z 2,k ,z w ) 2 Substituting h for p(h|ω 2,k ,z 2,k ,z w ), p(h|ω 2,k ,z 2,k ,z w ) w Replace with zw,k-1 We get p(h|ω 2,k ,z 2,k ,z w,k-1 );

[0026] Execute M-step: by maximizing Get the estimated value θ of θ1 at the kth iteration 1,k , Then we get the estimated value λ of the kth iteration k and z h The estimated value of z h,k , By maximizing Qz w (z w |z w,k-1 ), get the kth iteration z w The estimated value of z w,k , And simplify to get Among them, || represents the modulus value, represents the set of non-zero singular values ​​in Φ, s i represents the i-th non-zero singular value in Φ;

[0027] Step 5.b3: Determine whether k is less than K. If so, set k=k+1 and return to step 5.b2 to continue iterating; otherwise, end the iteration process; where the "=" in k=k+1 is an assignment symbol.

[0028] In step 1, the process of obtaining the time domain signal t is as follows: two binary transmission sequences are respectively recorded as b1 and b2; b1 is encoded by the encoder and then mapped to the frequency domain symbol x1 through QAM, and b2 is encoded by the encoder and then mapped to the frequency domain symbol x2 through QAM, wherein the encoder used for encoding b1 and b2 is the same, the QAM mapping method after b1 is encoded by the encoder and b2 is encoded by the encoder is the same, and the dimensions of x1 and x2 are N×1; x1 and x2 are superimposed as a signal u after power allocation, Among them, α1 and α2 are power allocation coefficients, α1>α2, α1, α2∈(0,1), α1+α2=1, and the dimension of u is N×1; the signal u is converted into a time domain signal t by inverse discrete Fourier transform, t=F * u, where F represents the normalized discrete Fourier transform matrix, the dimension of F is N×N, and F * is the Hermitian transform of F.

[0029] Compared with the prior art, the advantages of the present invention are:

[0030] 1) The method of the present invention uses the vector approximate message passing (VAMP) algorithm to approximate the estimated value of the channel impulse response vector based on the factor graph of the channel estimation. The singular value decomposition of the matrix is ​​used in the iterative process so that the matrix is ​​directly inverted only once at the beginning of the iteration, and the subsequent iterations do not need to perform the matrix direct inversion operation, which significantly reduces the computational complexity of the matrix direct inversion. At the same time, experiments show that compared with the traditional channel estimation method based on compressed sensing, such as the channel estimation method based on sparse Bayesian learning (SBL), the method of the present invention has the best channel estimation value and lower computational complexity.

[0031] 2) The method of the present invention uses the expectation maximization algorithm to calculate unknown hyperparameters. Compared with the traditional channel estimation method based on the approximate message passing (AMP) algorithm, it is no longer necessary to provide accurate channel prior information, which improves the applicability of the method of the present invention in practical applications.

[0032] 3) Compared with the traditional channel estimation method based on the generalized approximate message passing (GAMP) algorithm, the method of the present invention avoids the algorithm divergence problem caused by the elements in the observation matrix not satisfying the zero-mean Gaussian distribution, and improves the stability of the NOMA system. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 A system implementation block diagram of the method of the present invention;

[0034] Figure 2 A schematic diagram of a factor graph drawn according to joint probability in the method of the present invention;

[0035] Figure 3 It is a block diagram for realizing the vector approximate message passing algorithm in the method of the present invention;

[0036] Figure 4 The message passing flowchart in factor graph form for the vector approximate message passing algorithm;

[0037] Figure 5 It is a schematic diagram of the change of mean square error (MSE) of channel estimation under different signal-to-noise ratio (SNR) conditions when the total number of subcarriers is 256, the number of pilot subcarriers is 76, and the number of data subcarriers is 180 for the method of the present invention (EM-VAMP), the channel estimation method based on EM-SBL, and the channel estimation method based on GAMPSBL;

[0038] Figure 6 It is a schematic diagram of the change of the bit error rate (BER) of channel estimation under different signal-to-noise ratio (SNR) conditions when the total number of subcarriers is 256, the number of pilot subcarriers is 76, and the number of data subcarriers is 180 for the method of the present invention (EM-VAMP), the channel estimation method based on EM-SBL, and the channel estimation method based on GAMPSBL;

[0039] Figure 7 It is a schematic diagram of the change curve of the bit error rate (BER) of the downlink two-user NOMA system with the signal-to-noise ratio (SNR) when the total number of subcarriers is 256, the number of pilot subcarriers is 76, and the number of data subcarriers is 180 for the method of the present invention (EM-VAMP) and the channel estimation method based on GAMPSBL;

[0040] Figure 8 It is a schematic diagram of the change of mean square error (MSE) of channel estimation under different numbers of pilot subcarriers for the method of the present invention (EM-VAMP), the channel estimation method based on EM-SBL and the channel estimation method based on GAMPSBL when the signal-to-noise ratio (SNR) is 10 dB, the total number of subcarriers is 256, and the number of data subcarriers is 180;

[0041] Fig. 9 The schematic diagram is a diagram showing the variation of the algorithm running time with the number of iterations for the method of the present invention (EM-VAMP), the channel estimation method based on EM-SBL and the channel estimation method based on GAMPSBL when the signal-to-noise ratio (SNR) is 10 dB, the total number of subcarriers is 256, the number of pilot subcarriers is 76, and the number of data subcarriers is 180. DETAILED DESCRIPTION

[0042] The present invention is further described in detail below with reference to the accompanying drawings.

[0043] The present invention proposes a channel estimation method for a NOMA system based on vector approximate message transmission, and its overall implementation block diagram is as follows: Figure 1 As shown, it includes the following steps:

[0044] Step 1: At the transmitting end of the NOMA downlink communication system, the transmitting end sends a time domain signal t to the receiving end. The time domain signal t passes through the channel, where the dimension of t is N×1, N represents the total number of subcarriers in the NOMA system, N=N d +N p , N d Indicates the number of data subcarriers, N p Indicates the number of pilot subcarriers.

[0045] In this specific embodiment, the process of obtaining the time domain signal t is as follows: two binary transmission sequences are respectively recorded as b1 and b2; b1 is encoded by the encoder and then mapped to the frequency domain symbol x1 through QAM (quadrature amplitude modulation), and b2 is encoded by the encoder and then mapped to the frequency domain symbol x2 through QAM, wherein the encoder used for encoding b1 and b2 is the same, the QAM mapping method after b1 is encoded by the encoder and b2 is encoded by the encoder is the same, and the dimensions of x1 and x2 are N×1; x1 and x2 are superimposed as a signal u after power allocation, Among them, α1 and α2 are power allocation coefficients, α1>α2, α1, α2∈(0,1), α1+α2=1, and the dimension of u is N×1; the signal u is converted into a time domain signal t by an inverse discrete Fourier transform (IDFT), t=F * u, where F represents the normalized discrete Fourier transform matrix, the dimension of F is N×N, and F * is the Hermitian transform of F.

[0046] Step 2: At the receiving end of the NOMA downlink communication system, the receiving end receives the time domain signal r, r = Gt + ε; then r is converted from the time domain to the frequency domain through Discrete Fourier Transform (DFT) to obtain the frequency domain received signal R, In order to achieve the purpose of channel estimation, the frequency domain received signal R corresponding to the pilot subcarrier is obtained according to R p ,

[0047] Where G represents an N×N circulant matrix, the first L elements in the first column vector of G correspond to the L elements in the L×1 channel impulse response vector h, and the last NL elements are all 0, L represents the delay spread length of the channel, h=[h0,h1,...,h L-1 ] T , h0 represents the element with index 0 in h, h1 represents the element with index 1 in h, and h L-1 Indicates that the element with index L-1 in h is the element with index l in h. l , l=0,1,...,L-1, the superscript "T" indicates the transpose of a matrix or vector, ε indicates additive white Gaussian noise in the time domain, X is a diagonal matrix, X=diag(u), diag() indicates the diagonal matrix, and the diagonal elements of X are u1,u2,...,u N , u1 represents the first element in u, u2 represents the second element in u, u N represents the Nth element in u, F L represents the matrix composed of the first L column vectors of F, w represents the background noise vector in the frequency domain, w=Fε, satisfying CN(w;0,z w -1 I N ), indicating that w has a mean of 0 and a variance of z w -1 I N The complex Gaussian distribution, CN is the symbol for complex Gaussian distribution, I N represents the N-order identity matrix, z w-1 is the variance of w, X p represents a matrix consisting of row vectors corresponding to the pilot subcarrier index positions in X, w p represents a vector consisting of row elements in w corresponding to the pilot subcarrier index positions,

[0048] Step 3: According to the sparsity of h, h is assumed to obey an independent and identically distributed Bernoulli-Gaussian distribution, and the prior probability p(h|θ1) about h is obtained. Then according to Bayesian theory, we get p ,h joint probability p(R p ,h|θ), Reintroducing h 1 and h 2 , and let h 1 =h 2 =h, about R p ,h joint probability p(R p ,h|θ) can be rewritten as p ,h 1 ,h 2 The joint probability p(R p ,h 1 ,h 2 |θ),

[0049] Where θ1=[λ,z h ] T ,λ represents the sparsity of h,λ∈(0,1),z h represents the inverse variance of the non-zero elements in h, p(h l |θ1) represents the relationship between h l The prior probability of δ() is the Dirac function, CN(h l ; 0,z h -1 ) means h l Subject to mean 0 and variance z h -1 complex Gaussian distribution, θ is a hyperparameter set, p(R p |h,θ) means R under the condition that h,θ are known p The likelihood function of Represents R p Subject to mean Φh and variance The complex Gaussian distribution of N p The identity matrix of order, p(h 1 |θ1) represents the relationship between h 1 The prior probability of Represents Rp The mean is Φh 2 And the variance is The complex Gaussian distribution of .

[0050] Step 4: According to Plot the factor graph of the channel estimate, such as Figure 2 As shown, there are two variable nodes h in the factor graph 1 and h 2 、Three factor nodes p(h 1 |θ1) and δ(h 1 -h 2 )and

[0051] Step 5: Based on the factor graph, use the vector approximate message passing algorithm to approximate the estimated value of h. The specific process is as follows:

[0052] Step 5.a1: Let k represent the number of iterations, the initial value of k is 1; let K represent the maximum number of iterations.

[0053] Step 5.a2: During the kth iteration, calculate the factor node δ(h 1 -h 2 ) is passed to the variable node h 1 News Then according to p(h 1 |θ1) and Using the Bayesian method, we can get 1 The posterior probability p(h 1 |ω 1,k ,z 1,k ,θ1), Then Substitution And simplify it to get Further get h 1 The posterior expectation β 1,k and the mean variance Finally, we get h 1 The approximate belief b app (h 1 ),

[0054] Among them, CN(h 1 ;ω 1,k ,z 1,k -1 I L ) means h 1 The mean is ω 1,k And the variance is z 1,k -1 I L The complex Gaussian distribution, ω 1,krepresents the kth iteration h 1 Expectations, ω 1,1 =0, z 1,k -1 represents the kth iteration h 1 The variance of z 1,1 =1,I L represents the identity matrix of order L, Indicates h 1 The element with index l in Indicates about The prior probability of express Subject to mean 0 and variance z h -1 The complex Gaussian distribution, ρ k (l) represents ρ k The element with index l in ρ k is the introduced intermediate quantity, ω 1,k (l) represents ω 1,k The element with index l in CN(ω 1,k (l); 0, z 1,k -1 ) represents ω 1,k (l) has a mean of 0 and a variance of z 1,k -1 The complex Gaussian distribution, CN(ω 1,k (l); 0, (z 1,k -1 +z h -1 )) represents ω 1,k (l) has a mean of 0 and a variance of z 1,k -1 +z h -1 The complex Gaussian distribution of k (l) indicates a k The element with index l in a k is the introduced intermediate quantity, v k is the intermediate quantity introduced, express The mean is a k (l) and the variance is v k The complex Gaussian distribution of 1,k is the intermediate quantity introduced, Indicates h 1 The mean is β 1,k And the variance is The complex Gaussian distribution of .

[0055] Step 5.a3: During the kth iteration, calculate the variable node h 1 Passed to the factor node δ(h 1 -h 2 ) According to the sum-product algorithm, we can get

[0056] Among them, ω 2,k represents the kth iteration h 2 expectations, z 2,k -1 represents the kth iteration h 2 The variance of CN(h 1 ;ω 2,k ,z 2,k -1 I L ) means h 1 The mean is ω 2，k And the variance is z 2,k -1 I L The complex Gaussian distribution of .

[0057] Step 5.a4: During the kth iteration, the factor node δ(h 1 -h 2 ) is passed to the variable node h 2 News No change, obtained Then according to and Using the Bayesian method, we can get 2 The posterior probability And simplify to get h 2 The posterior expectation β 2,k and the mean variance β 2,k =(z w Φ H Φ+z 2,k I L ) -1 (z w Φ H R p +z 2,k ω 2,k ), Finally, we get h 2 The approximate belief b app (h 2 ),

[0058] Among them, CN(h 2 ;ω 2,k ,z2,k -1 I L ) means h 2 The mean is ω 2,k And the variance is z 2,k -1 I L The superscript "H" indicates the conjugate transpose of the matrix. 2,k is the intermediate quantity introduced, Tr[] represents matrix trace, Indicates h 2 The mean is β 2,k And the variance is The complex Gaussian distribution of .

[0059] Step 5.a5: Determine whether k is less than K. If so, execute step 5.a6; otherwise, obtain the minimum mean square error estimate of h according to the minimum mean square error (MMSE) estimation criterion. in, represents the estimated value of h, represents the true value of h, || ||2 represents the l2 norm, Indicates that in R p Under known conditions The posterior probability of Indicates that in R p Under known conditions expectations, Expressing expectation; based on β 1,K As an estimate of h

[0060] Step 5.a6: Calculate variable node h 2 Passed to the factor node δ(h 1 -h 2 ) According to the sum-product algorithm, we can get Among them, CN(h 2 ;ω 1,k+1 ,z 1,k+1 -1 I L ) means h 2 The mean is ω 1,k+1 And the variance is z 1,k+1 -1 I L The complex Gaussian distribution of In order to reduce β 2,k With φ 2,k The computational complexity caused by matrix inversion in the expression of , and then β 2,k With φ 2,kThe singular value decomposition (SVD) of Φ in the expression of Φ is obtained to obtain Φ = USV H , further obtained Where S is a diagonal matrix, U and V are both unitary matrices, Y represents the rank of Φ, Y = rank (Φ), rank () means to find the rank of the matrix, d k is the intermediate quantity introduced, d k =z w diag(z w S H S+z 2,k 1 L ) -1 S H S, 1 L represents an L×L dimensional all-one matrix, d k (l) indicates d k The element with index l in ; then let k=k+1, return to step 5.a2 to continue iteration, where the "=" in k=k+1 is the assignment symbol.

[0061] Figure 3 The implementation block diagram of the vector approximate message passing algorithm in the method of the present invention is given. Figure 4 The message passing flowchart in factor graph form of the vector approximate message passing algorithm is given.

[0062] The expectation maximization (EM) algorithm is used to calculate the estimated value of θ. The specific process is as follows:

[0063] Step 5.b1: Let k represent the number of iterations, with the initial value of k being 1; let K represent the maximum number of iterations.

[0064] Step 5.b2: During the kth iteration, perform E-step: 1,k-1 Under known conditions, calculate the joint probability p(R p ,h|θ) in the posterior probability p(h|ω 1,k ,z 1,k ,θ 1,k-1 ) θ1 (θ1|θ 1,k-1 ), In z w,k-1 Under known conditions, calculate the joint probability p(R p ,h|θ) in the posterior probability p(h|ω 2,k ,z 2,k ,z w,k-1 ) When k = 1, θ 1,k-1represents the initial value of θ1, which is a random number. When k>1, θ 1,k-1 represents the estimated value of θ1 at the k-1th iteration, and p(h 1 |ω 1,k ,z 1,k ,θ1) in h 1 Substituting h for p(h|ω 1,k ,z 1,k ,θ1), and p(h|ω 1,k ,z 1,k ,θ1) is replaced by θ 1,k-1 We get p(h|ω 1,k ,z 1,k ,θ 1,k-1 ), ln() represents the natural logarithm function, p(R p |h,z w ) means in h,z w Under known conditions, R p The likelihood function of π=3.14…, exp() represents the exponential function with the natural constant e as the base, e=2.71…, when k=1, z w,k-1 Indicates z w The initial value of z w The initial value of z is preset, when k>1 w,k-1 represents the k-1th iteration z w The estimated value of p(h 2 |ω 2,k ,z 2,k ,z w ) 2 Substituting h for p(h|ω 2,k ,z 2,k ,z w ), p(h|ω 2,k ,z 2,k ,z w ) w Replace with z w,k-1 We get p(h|ω 2,k ,z 2,k ,z w,k-1 ).

[0065] Execute M-step: by maximizing Get the estimated value θ of θ1 at the kth iteration 1,k ,

[0066] Then we get the estimated value λ of the kth iteration k and z h The estimated value of z h,k ,

[0067] By maximizing Get the kth iteration z w The estimated value of z w,k , And simplify to get Among them, || represents the modulus value, represents the set of non-zero singular values ​​in Φ, s i represents the i-th non-zero singular value in Φ.

[0068] Step 5.b3: Determine whether k is less than K. If so, set k=k+1 and return to step 5.b2 to continue iterating; otherwise, end the iteration process; where the "=" in k=k+1 is an assignment symbol.

[0069] In order to further illustrate the effectiveness and feasibility of the method of the present invention, a computer simulation was performed on the method of the present invention, and the detailed simulation parameters are listed in Table 1.

[0070] Table 1 Simulation parameters list

[0071]

[0072] In order to compare the performance of different channel estimation methods, the bit error rate (BER) of channel estimation is defined. b e Indicates the number of error bits, b all represents the total number of bits; defines the mean square error (MSE) of channel estimation, represents the true value of h, represents the estimated value of h; the signal-to-noise ratio (SNR) of the NOMA system is defined as the ratio of the normalized transmitted signal power to the total noise power of the system, that is, Under the same simulation environment, the performance of the method of the present invention, the channel estimation method based on expectation maximization-sparse Bayesian learning (EM-SBL) and the channel estimation method based on generalized approximate message passing sparse Bayesian learning (GAMPSBL) are compared in terms of the mean square error (MSE) of channel estimation and the bit error rate (BER) of channel estimation.

[0073] Figure 5The variation of mean square error (MSE) of channel estimation under different signal-to-noise ratio (SNR) conditions is given for the method of the present invention (EM-VAMP), the channel estimation method based on EM-SBL and the channel estimation method based on GAMPSBL when the total number of subcarriers is 256, the number of pilot subcarriers is 76 and the number of data subcarriers is 180. Figure 5 It can be seen that the MSE of the three methods all show a downward trend with the increase of the signal-to-noise ratio (SNR). The method of the present invention shows a lower MSE under all SNRs, demonstrating its superior performance in channel estimation. For example, when MSE=10 -4 When the signal-to-noise ratio gain of the method of the present invention is about 5 dB compared with the channel estimation method based on EM-SBL, and is about 7 dB compared with the channel estimation method based on GAMPSBL.

[0074] Figure 6 The changes of the bit error rate (BER) of the channel estimation under different signal-to-noise ratio (SNR) conditions are given for the method of the present invention (EM-VAMP), the channel estimation method based on EM-SBL and the channel estimation method based on GAMPSBL when the total number of subcarriers is 256, the number of pilot subcarriers is 76 and the number of data subcarriers is 180. Figure 6 It can be seen that the BER of the three methods decreases with the increase of the signal-to-noise ratio (SNR), and the method of the present invention shows a lower BER under all SNRs. -2 When the signal-to-noise ratio gain of the method of the present invention is about 1 dB compared with the channel estimation method based on EM-SBL, and is about 1.5 dB compared with the channel estimation method based on GAMPSBL.

[0075] Figure 7 The variation curves of the bit error rate (BER) of the downlink two-user NOMA system with the signal-to-noise ratio (SNR) for the method of the present invention (EM-VAMP) and the channel estimation method based on GAMPSBL are given when the total number of subcarriers is 256, the number of pilot subcarriers is 76, and the number of data subcarriers is 180. Figure 7 It can be seen that the BER of long-distance user 1 is lower than that of short-distance user 2. This is because long-distance user 1 is far from the base station and has poor channel conditions, so it obtains a larger power allocation coefficient and has no accumulated error in the SIC (Successive Interference Cancellation) process. In the two-user scenario, the BER of the method of the present invention is the lowest in long-distance user 1 or short-distance user 2. At BER=10 -2When , the signal-to-noise ratio gain of the method of the present invention in the long-distance user 1 is about 1.5 dB compared with the channel estimation method based on GAMPSBL, and the signal-to-noise ratio gain of the method in the short-distance user 2 is about 1.5 dB compared with the channel estimation method based on GAMPSBL.

[0076] Figure 8 The variation of mean square error (MSE) of channel estimation under different numbers of pilot subcarriers for the method of the present invention (EM-VAMP), the channel estimation method based on EM-SBL and the channel estimation method based on GAMPSBL is given when the signal-to-noise ratio (SNR) is 10dB, the total number of subcarriers is 256 and the number of data subcarriers is 180. Figure 8 It can be seen that the channel estimation performance of the method of the present invention under different pilot subcarriers is better than that of other methods. When the number of pilot subcarriers is small, the observation values ​​that can be effectively used by the observation matrix of compressed sensing are small, and the method of the present invention has better performance with a small number of observation values. As the number of pilot subcarriers gradually increases, the MSE of all methods is improved, but the method of the present invention still has the best channel estimation performance.

[0077] Fig. 9 The following is a schematic diagram showing the variation of the running time of the method (EM-VAMP), the channel estimation method based on EM-SBL and the channel estimation method based on GAMPSBL with the number of iterations when the signal-to-noise ratio (SNR) is 10 dB, the total number of subcarriers is 256, the number of pilot subcarriers is 76, and the number of data subcarriers is 180. Fig. 9 It can be seen that the running time of all channel estimation methods gradually increases with the increase of the number of iterations. Compared with the channel estimation method based on EM-SBL, the running time of the method of the present invention is greatly reduced. In all iteration number regions, the average running time of the method of the present invention is about 9 seconds shorter than the channel estimation method based on EM-SBL and about 0.8 seconds shorter than the channel estimation method based on GAMPSBL. Among all channel estimation methods, the method of the present invention shows the lowest computational complexity and the shortest running time, which verifies the low time complexity and obvious real-time advantages of the method of the present invention in channel estimation.

Claims

1. A channel estimation method for NOMA system based on vector approximate message passing, characterized in that The following steps are involved: Step 1: At the transmitting end of the NOMA downlink communication system, the transmitting end sends a time domain signal t to the receiving end. The time domain signal t passes through the channel, where the dimension of t is N×1, N represents the total number of subcarriers in the NOMA system, N=N d +N p , N d Indicates the number of data subcarriers, N p Indicates the number of pilot subcarriers; Step 2: At the receiving end of the NOMA downlink communication system, the receiving end receives the time domain signal r, r = Gt + ε; then r is converted from the time domain to the frequency domain through discrete Fourier transform to obtain the frequency domain received signal R, Then, according to R, we get the frequency domain received signal R corresponding to the pilot subcarrier. p , Where G represents an N×N circulant matrix, the first L elements in the first column vector of G correspond to the L elements in the L×1 channel impulse response vector h, and the last NL elements are all 0, L represents the delay spread length of the channel, h=[h0,h1,…,h L-1 ] T , the index in h is The elements of The superscript "T" indicates the transpose of a matrix or vector, ε indicates additive white Gaussian noise in the time domain, and F indicates the normalized discrete Fourier transform matrix. The dimension of F is N×N. * is the Hermitian transform of F. The two binary transmission sequences at the transmitter are encoded by the encoder and then mapped into frequency domain symbols through QAM. The two frequency domain symbols are superimposed into signal u after power allocation. X is a diagonal matrix, X = diag(u), diag() means to find the diagonal matrix, and the diagonal elements of X are u1,u2,…,u N , u1 represents the first element in u, u2 represents the second element in u, u N represents the Nth element in u, F L represents the matrix composed of the first L column vectors of F, w represents the background noise vector in the frequency domain, w=Fε, satisfying CN(w;0,z w -1 I N ), CN is the symbol for complex Gaussian distribution, I N represents the N-order identity matrix, z w -1 is the variance of w, X p represents a matrix consisting of row vectors corresponding to the pilot subcarrier index positions in X, w p represents a vector consisting of row elements in w corresponding to the pilot subcarrier index positions, Step 3: According to the sparsity of h, h is assumed to obey an independent and identically distributed Bernoulli-Gaussian distribution, and the prior probability p(h|θ1) about h is obtained. Then according to Bayesian theory, we get p ,h joint probability p(R p ,h|θ), Reintroducing h 1 and h 2 , and let h 1 =h 2 =h, about R p ,h joint probability p(R p ,h|θ) can be rewritten as p ,h 1 ,h 2 The joint probability p(R p ,h 1 ,h 2 |θ), Where θ1=[λ,z h ] T ,λ represents the sparsity of h,λ∈(0,1),z h represents the inverse variance of the non-zero elements in h, Indicates about The prior probability of , δ() is the Dirac function, Indicates h l It has a mean of 0 and a variance of z. h -1 The complex Gaussian distribution, θ = [θ1 T ,z w ] T , p(R p |h,θ) means R under the condition that h,θ are known p The likelihood function of It means that Rp has a mean of Φh and a variance of The complex Gaussian distribution of N p The identity matrix of order, p(h 1 |θ1) represents the relationship between h 1 The prior probability of Represents R p The mean is Φh 2 And the variance is The complex Gaussian distribution of Step 4: According to Draw the factor graph of channel estimation, which has two variable nodes h 1 and h 2 、Three factor nodes p(h 1 |θ1) and δ(h 1 -h 2 )and Step 5: Based on the factor graph, use the vector approximate message passing algorithm to approximate the estimated value of h. The specific process is as follows: Step 5.a1: Let k represent the number of iterations, the initial value of k is 1; let K represent the maximum number of iterations; Step 5.a2: During the kth iteration, calculate the factor node δ(h 1 -h 2 ) is passed to the variable node h 1 News Then according to p(h 1 |θ1) and Using the Bayesian method, we can get 1 The posterior probability p(h 1 |ω 1,k ,z 1,k ,θ1), Then Substitution And simplify it to get Further get h 1 The posterior expectation β 1,k and the mean variance Finally, we get h 1 The approximate belief b app (h 1 ), Among them, CN(h 1 ;ω 1,k ,z 1,k -1 I L ) means h 1 The mean is ω 1,k And the variance is z 1,k -1 I L The complex Gaussian distribution, ω 1,k represents the kth iteration h 1 Expectations, ω 1,1 =0, z 1,k -1 represents the kth iteration h 1 The variance of z 1,1 =1,I L represents the identity matrix of order L, Indicates h 1 The element with index l in Indicates about The prior probability of express It has a mean of 0 and a variance of z. h -1 The complex Gaussian distribution of Represents ρ k The element with index l in ρ k is the intermediate quantity introduced, Represents ω 1,k The element with index l in express It has a mean of 0 and a variance of z. 1,k -1 The complex Gaussian distribution of 0,(z 1,k -1 +z h -1 ))express It has a mean of 0 and a variance of z. 1,k -1 +z h -1 The complex Gaussian distribution of Indicates a k The index is Elements of a k is the introduced intermediate quantity, v k is the intermediate quantity introduced, express The mean is And the variance is v k The complex Gaussian distribution of 1,k is the intermediate quantity introduced, Indicates h 1 The mean is β 1,k And the variance is The complex Gaussian distribution of Step 5.a3: During the kth iteration, calculate the variable node h 1 Passed to the factor node δ(h 1 -h 2 ) According to the sum-product algorithm, we can get Among them, ω 2,k represents the kth iteration h 2 expectations, z 2,k -1 represents the kth iteration h 2 The variance of CN(h 1 ;ω 2,k ,z 2,k -1 I L ) means h 1 The mean is ω 2,k And the variance is z 2,k -1 I L The complex Gaussian distribution of Step 5.a4: During the kth iteration, the factor node δ(h 1 -h 2 ) is passed to the variable node h 2 News No change, obtained Then according to and Using the Bayesian method, we can get 2 The posterior probability p(h 2 |ω 2,k ,z 2,k ,z w ), And simplify to get h 2 The posterior expectation β 2,k and the mean variance β 2,k =(z w Φ H Φ+z 2,k I L ) -1 (z w Φ H R p +z 2,k ω 2,k ), Finally, we get h 2 The approximate belief b app (h 2 ), Among them, CN(h 2 ;ω 2,k ,z 2,k -1 I L ) means h 2 The mean is ω 2,k And the variance is z 2,k -1 I L The superscript "H" indicates the conjugate transpose of the matrix. 2,k is the intermediate quantity introduced, Tr[] represents matrix trace, Indicates h 2 The mean is β 2,k And the variance is The complex Gaussian distribution of Step 5.a5: Determine whether k is less than K. If so, execute step 5.a6; otherwise, according to the minimum mean square error estimation criterion, obtain the minimum mean square error estimate of h. in, represents the estimated value of h, represents the true value of h, || ||2 represents the l2 norm, Indicates that in R p Under known conditions The posterior probability of Indicates that in R p Under known conditions expectations, Expressing expectation; based on β 1,K As an estimate of h Step 5.a6: Calculate variable node h 2 Passed to the factor node δ(h 1 -h 2 ) According to the sum-product algorithm, we can get Among them, CN(h 2 ;ω 1,k+1 ,z 1,k+1 -1 I L ) means h 2 The mean is ω 1,k+1 And the variance is z 1,k+1 -1 I L The complex Gaussian distribution of Then for β 2,k With φ 2,k Perform singular value decomposition on Φ in the expression of , and get Φ = USV H , further obtained Among them, S is a diagonal matrix, U and V are both unitary matrices, Y represents the rank of Φ, and d k is the intermediate quantity introduced, 1 L represents an L×L dimensional all-one matrix, Indicates d k The element with index l in ; then let k = k + 1, return to step 5.a2 to continue iterating, where "=" in k = k + 1 is the assignment symbol; The expectation maximization algorithm is used simultaneously to calculate the estimated value of θ. The specific process is as follows: Step 5.b1: Let k represent the number of iterations, the initial value of k is 1; let K represent the maximum number of iterations; Step 5.b2: During the kth iteration, perform E-step: 1,k-1 Under known conditions, calculate the joint probability p(R p ,h|θ) in the posterior probability p(h|ω 1,k ,z 1,k ,θ 1,k-1 ) In z w,k-1 Under known conditions, calculate the joint probability p(R p ,h|θ) in the posterior probability p(h|ω 2,k ,z 2,k ,z w,k-1 ) When k = 1, θ 1,k-1 represents the initial value of θ1. When k>1, θ 1,k-1 represents the estimated value of θ1 at the k-1th iteration, and p(h 1 |ω 1,k ,z 1,k ,θ1) in h 1 Substituting h for p(h|ω 1,k ,z 1,k ,θ1), and p(h|ω 1,k ,z 1,k ,θ1) is replaced by θ 1,k-1 We get p(h|ω 1,k ,z 1,k ,θ 1,k-1 ), ln() represents the natural logarithm function, p(R p |h,z w ) means in h,z w The likelihood function of Rp under known conditions is: exp() represents an exponential function with the natural constant e as the base. When k = 1, z w,k-1 Indicates z w The initial value of z is w,k-1 represents the k-1th iteration z w The estimated value of p(h 2 |ω 2,k ,z 2,k ,z w ) 2 Substituting h for p(h|ω 2,k ,z 2,k ,z w ), p(h|ω 2,k ,z 2,k ,z w ) w Replace with z w,k-1 We get p(h|ω 2,k ,z 2,k ,z w,k-1 ); Execute M-step: by maximizing Get the estimated value θ of θ1 at the kth iteration 1,k , Then we get the estimated value λ of the kth iteration k and z h The estimated value z h,k , By maximizing Get the kth iteration z w The estimated value z w,k , And simplify to get Among them, || represents the modulus value, represents the set of non-zero singular values ​​in Φ, s i represents the i-th non-zero singular value in Φ; Step 5.b3: Determine whether k is less than K. If so, set k=k+1 and return to step 5.b2 to continue iterating; otherwise, end the iteration process; where "=" in k=k+1 is an assignment symbol.

2. The channel estimation method of the NOMA system based on vector approximate message passing according to claim 1 is characterized in that In step 1, the process of obtaining the time domain signal t is as follows: two binary transmission sequences are respectively recorded as b1 and b2; b1 is encoded by the encoder and then mapped to the frequency domain symbol x1 through QAM, and b2 is encoded by the encoder and then mapped to the frequency domain symbol x2 through QAM, wherein the encoder used for encoding b1 and b2 is the same, the QAM mapping method after b1 is encoded by the encoder and b2 is encoded by the encoder is the same, and the dimensions of x1 and x2 are N×1; x1 and x2 are superimposed as a signal u after power allocation, Among them, α1 and α2 are power allocation coefficients, α1>α2, α1, α2∈(0,1), α1+α2=1, and the dimension of u is N×1; the signal u is converted into a time domain signal t by inverse discrete Fourier transform, t=F * u, where F represents the normalized discrete Fourier transform matrix, the dimension of F is N×N, F * is the Hermitian transform of F.

Images