阅读说明：本技术 一种改进的近似复指数基扩展的联合反馈k-means信道估计方法 (Improved joint feedback k-means channel estimation method for approximate complex exponential basis expansion ) 是由 姜斌 周志杰 包建荣 刘超 唐向宏 于 2020-12-25 设计创作，主要内容包括：本发明公开一种改进的近似复指数基扩展的联合反馈k-means信道估计方法,包括步骤：S1.对发射机发送的数据采用FTCA-CE-BEM方法,得到该信道模型下的接收信号；S2、利用迫零均衡方法,得到消除ICI后的接收信号；S3.对消除ICI的接收信号,采用最小二乘法计算CIR的估计h-(LS)(n),并采用密度参数,删除孤立CIR,得非孤立CIR：h′-(LS)(n),S4.将h′-(LS)(n)划分为噪声类和信号类,并计算各自类的初始聚类中心；S5.设置判别函数,以对所有CIR进行判决并重新分类以及计算聚类中心；S6.判断聚类结果是否改变,若是则根据改变后的结果回到步骤S5,若否则将判为噪声类的h’-(LS)(n)置零,得到时域信道函数h′-(k-means)(n)；S7.对h′-(k-means)(n)执行DFT计算,得到反馈后的频域信道函数H′-(k-means)(k)。本发明适合高速移动通信中的高精度信道估计且估计准确度高。(The invention discloses an improved joint feedback k-means channel estimation method of approximate complex exponential base expansion, which comprises the following steps: s1, obtaining a receiving signal under the channel model by adopting an FTCA-CE-BEM method for data sent by a transmitter; s2, obtaining a received signal with ICI eliminated by using a zero-forcing equalization method; s3, calculating the estimation h of the CIR by adopting a least square method for the received signal for eliminating the ICI LS (n), and deleting the isolated CIR by adopting a density parameter to obtain a non-isolated CIR: h' LS (n), S4. mixing h' LS (n) dividing the signal into a noise class and a signal class, and calculating initial clustering centers of the respective classes; s5, setting a discrimination function to judge all CIRs, reclassifying and calculating a clustering center; s6, judging whether the clustering result is changed or not, if so, returning to the original clustering result according to the changed resultStep S5, if not, judging the h 'as noise' LS (n) setting zero to obtain a time domain channel function h' k‑means (n); s7, h 'is paired' k‑means (n) performing DFT calculation to obtain a frequency domain channel function H 'after feedback' k‑means (k) In that respect The method is suitable for high-precision channel estimation in high-speed mobile communication and has high estimation accuracy.)

1. An improved k-means channel estimation method with approximate complex exponential base spreading, which is characterized by comprising the following steps:

s1, aiming at data sent by a transmitter, obtaining a receiving signal under a channel model by adopting a fractional tap channel approximate complex exponential basis expansion model method, wherein the fractional tap channel approximate complex exponential basis expansion model is FTCA-CE-BEM;

s2, obtaining a received signal with ICI eliminated by using a zero-forcing equalization method, wherein ICI is sub-carrier interference;

s3, aiming at the ICI-eliminated received signal obtained in the step S2, adopting the minimum twoMultiplication of the estimate h of the CIR_LS(n), simultaneously, deleting the isolated CIR by adopting a density parameter to obtain a non-isolated CIR: h'_LS(n), wherein the density parameter is the number of CIRs within a distance r from a certain CIR, the CIR is the channel impulse response, and r is the radius of the selected sphere;

s4, h 'obtained in step S3'_LS(n) dividing the signal into a noise class and a signal class, and calculating initial clustering centers of the respective classes;

s5, setting a discrimination function of the distance from the CIR to the clustering center to judge and reclassify all the CIRs and calculate the clustering center;

s6, judging whether the clustering result is changed or not, if so, returning to the step S5 according to the changed result, and if not, judging the clustering result to be h 'of the noise class'_LS(n) set to zero, and h 'determined to be signal-based by bonding'_LS(n) obtaining a time domain channel function h'_k-means(n)；

S7, h 'is paired'_k-means(n) performing a DFT calculation to obtain a fed back frequency domain channel function H'_k-means(k) I.e., the final channel estimation result, DFT is a discrete fourier transform.

2. The improved k-means channel estimation method for joint feedback of approximate complex exponential base extension according to claim 1, wherein in step S1, the FTCA-CE-BEM method is adopted to obtain the received signal under the channel model, which comprises the following steps:

s1.1, transmit data X ═ X (0), X (1), …, X (N-1) to the transmitter]^TPerforming an N-point inverse discrete fourier transform to obtain a time domain signal;

s1.2, performing discrete Fourier transform on time domain data y (n) received by a receiving end to obtain a frequency domain signal, and calculating to obtain k subcarrier receiving signals:

h (n, L) is a real number and is expressed as a sampling value of the fast time-varying channel impulse response at the nth time and the first path, L is the total path number of information transmission, and x (n-L) is the input of the n-L time; w (n) is the mean value of zero and the variance of σ²Is a desired channel function, w (k) is frequency domain channel additive noise, i (k) represents ICI caused by a fast time varying channel, and is expressed as:

wherein, the expression of the channel function H (k, m) is:

s1.3, adopting a fractional tap channel approximation method, and introducing a fractional weighting factor K_αTo simulate the non-sampling interval channel to obtain the frequency domain function H of the actual channel_FTCA(k)：

H_FTCA(k)＝H₁(k)+H_e(k) (16)

Wherein H₁(k)＝DFT[h₁(τ)]；H_e(k) Channel estimation error for using FTCA method; g (l) represents the tap coefficients of the filter;indicating the number of taps of the FTCA filter,denotes rounding up, τ_maxRepresenting the maximum time delay, f_sRepresenting the sampling frequency, T_sRepresents a sampling interval;

s1.4, according to equations (15) and (17), the channel function H (k, m) is expressed as:

s1.5, representing the tap coefficient g (n, l) in equation (18) in the BEM form, then:

wherein the base coefficient g_q(l) A weighting coefficient representing the Q-th complex exponential base on the l-th path, Q being a natural number, andrepresenting the order of the basis extension model, f_maxRepresents the maximum doppler shift;

s1.6, substituting formula (20) into formula (18) to obtain:

when m-k + Q-Q/2 is an integer, H (k, m) ≠ 0, and 0. ltoreq. m, k. ltoreq. N-1, so equation (21) is simplified as:

given k (0. ltoreq. k < N) by equation (22), each Q (0. ltoreq. Q. ltoreq. Q) corresponds to m, so that H (k, m) is not zero, and the non-zero element position of the k-th row of the channel matrix corresponding to H (k, m) is defined as:

s1.7, obtaining a signal expression received by a channel model in the FTCA-CE-BEM method according to the formula (22), the formula (23) and the formula (13), wherein the signal expression is as follows:

3. the improved k-means channel estimation method with approximate complex exponential-based spread joint feedback according to claim 2, wherein CIR of multipath channel is expressed as:wherein h (l) is channel gain, T_sRepresents the sampling interval, i_lT_sThe channel delay of the l-th path is represented, and δ (n) represents an impulse function.

4. The improved k-means channel estimation method with joint feedback approximating complex exponential base spreading as claimed in claim 2, wherein the base coefficient g in formula (24) in step S1.7_q(l) The calculation steps are as follows:

A. the complex exponential basis coefficients are represented by a vector g of dimension (Q +1) mx 1, and are expressed as:

wherein, g_q＝[g_q(0),g_q(1),…,g_q(M-1)]^TQ is more than or equal to 0 and less than or equal to Q, and represents a coefficient vector with dimension of M multiplied by 1 corresponding to the qth basis function;

B. matrix with dimension N × NF_MRepresenting a matrix consisting of the first M columns of F, then F_MThe k-th row of expressions of (a) is:

C. the equation (24) is derived as:

D. using P pilot signals to estimate g, and their positions are k (1), k (2), …, k (P), respectively, according to equation (27), P linear equations can be expressed as:

E. defining a pilot sequence:

andis a P x 1 dimensional real number vector, Q ∈ {0,1, …, Q }, and equation (27) is changed to:

Y＝[diag(X₀)F₀,…,diag(X_Q)F_Q]g+W＝Ag+W (29)

wherein the content of the first and second substances,

A＝[diag(X₀)F₀,…,diag(X_Q)F_Q],Y＝[Y(k(1)),…,Y(k(P))]^T,W＝[W(k(1)),…,W(k(P))]^T，

diag[·]representing a diagonal matrix, i.e. the elements not on the main diagonal are all zero, the elements on the main diagonal andq belongs to the element one-to-one correspondence in {0,1, …, Q };

F. from equation (29), a least squares estimate of g is obtained as:

g_LS＝(A^HA)^-1A^HY (30)。

5. the improved k-means channel estimation method by approximate complex exponential base expansion joint feedback according to claim 4, wherein in step S2, the zero-forcing equalization method is used to obtain the received signal after ICI cancellation, which includes the following steps:

s2.1, obtaining a base coefficient least square estimation quantity g according to the formula (30) to substitute the formula (22) and obtaining a channel function H (k, m);

s2.2 according to zero-forcing equalisation methods, i.e. X ═ H^-1(k, m) Y, estimating an input signal X from the received signal Y, and substituting X into formula (14) to obtain an inter-subcarrier interference value I' (k);

s2.3, using FTCA-CE-BEM method, the received signal Y (k) is obtained from equation (13), and then the received signal Y' (k) with ICI removed is calculated as:

Y′(k)＝Y(k)-I′(k)＝H(k)X(k)+O(k) (31)

where o (k) ═ w (k) + I (k) -I' (k) denotes the remaining ICI and channel noise.

6. The improved k-means channel estimation method for approximation of complex exponential-based spread joint feedback according to claim 5, wherein in step S3, the least square method is used to calculate the estimated h of CIR_LS(n), comprising the steps of:

s3.1, channel estimation according to the least-squares method, i.e.And according to formula (31) in step S2.3, the channel function expression is:

s3.2. for H_LS(k) Performing an inverse discrete fourier transform, namely:

n is more than or equal to 0 and less than or equal to N-1, so as to obtain:

h_LS(n)＝h(n)+o(n)，0≤n≤N-1 (33)

wherein o (n) ═ IDFT [ o (k)/x (k) ].

7. The improved k-means channel estimation method for joint feedback of approximate complex exponential base spread according to claim 6, wherein the step S4 includes: h 'are'_LS(n) is divided into three parts, n is more than or equal to 0 and less than or equal to L_cp-1 and N-L_cpN is not more than N and not more than N'_LS(n) as signal class training samples, apply L_cp≤n≤N-L_cp-1 part of h'_LS(n) as noise-like training samples, L_cpIndicating the cyclic prefix length.

8. The improved k-means channel estimation method with approximate complex exponential based spread combination as claimed in claim 7, wherein in step S4:

computing signal class initial clustering centersThe expression of (a) is:

computing noise-like initial cluster centersThe expression of (a) is:

9. the improved k-means channel estimation method with approximate complex exponential base spreading combination as claimed in claim 8, wherein the step S5 comprises the following steps:

s5.1, respectively calculating signal class h'_LS(n) and noise-like h'_LS(n) distances to the various cluster centers, respectivelyAndn is more than or equal to 0 and less than or equal to N-1. The expression is as follows:

s5.2, setting a discriminant function r (n), wherein the expression is as follows:

s5.3, if r (n) is less than or equal to 0, corresponding h'_LS(n) is judged to be a signal type, and if r (n) > 0, corresponding h'_LSAnd (n) judging the cluster centers as noise classes, and recalculating the updated cluster centers respectively.

10. The improved k-means channel estimation method for joint feedback of approximate complex exponential base spreading as claimed in claim 9, wherein the step S7 specifically includes:

for h 'obtained in step S5'_k-means(n) performing a discrete fourier transform, namely:

0≤k≤N-1，

obtaining a frequency domain channel function H 'after feedback'_k-means(k) I.e. the final channel estimation result, is:

H′_k-means(k)＝DFT[h′_k-means(n)] (39)。

Technical Field

The invention belongs to the field of digital communication, and particularly relates to an improved joint feedback k-means channel estimation method based on approximate complex exponential expansion.

Background

Channel estimation is a process of estimating model parameters of a certain assumed channel model from received data. If the channel is linear, then the channel estimate is an estimate of the system impulse response. It is emphasized that channel estimation is a mathematical representation of the effect of the channel on the input signal, while "good" channel estimation is an estimation algorithm that minimizes some estimation error.

The performance of a wireless communication system is greatly affected by wireless channels, such as shadow fading and frequency selective fading, so that the propagation path between a transmitter and a receiver is very complicated. Wireless channels are not fixed and predictable as wired channels, but rather have a large degree of randomness, which presents a significant challenge to the design of a receiver. In coherent detection of an OFDM system, a channel needs to be estimated, and the accuracy of channel estimation directly affects the performance of the whole system. In order to accurately recover a transmission signal at a receiving end, people adopt various measures to resist the influence of multipath effect on a transmission signal, and the realization of a channel estimation technology needs to know information of a wireless channel, such as the order of the channel, the doppler shift, the multipath delay or the impulse response of the channel. Therefore, channel parameter estimation is a key technology for implementing a wireless communication system. Whether detailed channel information can be obtained or not is an important index for measuring the performance of a wireless communication system, so that a transmitting signal can be correctly demodulated at a receiving end. Therefore, the research on the channel parameter estimation algorithm is a significant work.

Most of the traditional channel estimation methods are directed at time-invariant channels, but with the development of high-speed railways in China, particularly when the relative moving speed reaches more than 300km/h, wireless communication cannot be completed in high quality, and in a fast time-variant OFDM system, the orthogonality of subcarriers is destroyed, so that the interference among the subcarriers is caused, and therefore a new method needs to be researched to adapt to the fast time-variant channel environment.

Disclosure of Invention

The invention discloses a joint feedback k-MEANS (FTCA-CE-BEM-k-MEANS) channel estimation method based on improved approximate complex exponential basis expansion in a fast time-varying environment and an Orthogonal Frequency Division Multiplexing (OFDM) system, aiming at the problem of performance deterioration of a wireless communication system caused by subcarrier interference (ICI) due to Doppler frequency shift increase in a high-speed mobile environment.

The invention adopts the following technical scheme: an improved joint feedback k-means channel estimation method of approximate complex exponential base expansion comprises the following steps:

s1, aiming at data sent by a transmitter, obtaining a receiving signal under a channel model by adopting a fractional tap channel approximate complex exponential basis expansion model method, wherein the fractional tap channel approximate complex exponential basis expansion model is FTCA-CE-BEM;

s2, obtaining a received signal with ICI eliminated by using a zero-forcing equalization method, wherein ICI is sub-carrier interference;

s3, aiming at the received signal obtained in the step S2 after ICI elimination, calculating the estimation h of the CIR by adopting a least square method_LS(n), simultaneously, deleting the isolated CIR by adopting a density parameter to obtain a non-isolated CIR: h'_LS(n), wherein the density parameter is the number of CIRs within a distance r from a certain CIR, the CIR is the channel impulse response, and r is the radius of the selected sphere;

s4, h 'obtained in step S3'_LS(n) dividing the signal into a noise class and a signal class, and calculating initial clustering centers of the respective classes;

s5, setting a discrimination function of the distance from the CIR to the clustering center to judge and reclassify all the CIRs and calculate the clustering center;

s6, judging whether the clustering result is changed, if so, returning to the step S5 according to the changed result, otherwise,h 'judged to be noise'_LS(n) set to zero, and h 'determined to be signal-based by bonding'_LS(n) obtaining a time domain channel function h'_k-means(n)；

S7, h 'is paired'_k-means(n) performing a DFT calculation to obtain a fed back frequency domain channel function H'_k-means(k) I.e., the final channel estimation result, DFT is a discrete fourier transform.

As a preferred scheme, in step S1, obtaining the received signal under the channel model by using the FTCA-CE-BEM method specifically includes the following steps:

s1.1, transmit data X ═ X (0), X (1), …, X (N-1) to the transmitter]^TPerforming an N-point inverse discrete fourier transform to obtain a time domain signal;

s1.2, performing discrete Fourier transform on time domain data y (n) received by a receiving end to obtain a frequency domain signal, and calculating to obtain k subcarrier receiving signals:

h (n, L) is a real number and is expressed as a sampling value of the fast time-varying channel impulse response at the nth time and the first path, L is the total path number of information transmission, and x (n-L) is the input of the n-L time; w (n) is the mean value of zero and the variance of σ²Is a desired channel function, w (k) is frequency domain channel additive noise, i (k) represents ICI caused by a fast time varying channel, and is expressed as:

wherein, the expression of the channel function H (k, m) is:

s1.3, adopting a fractional tap channel approximation method, and introducing a fractional weighting factor K_αTo simulate the non-sampling interval channel to obtain the frequency domain function H of the actual channel_FTCA(k)：

H_FTCA(k)＝H₁(k)+H_e(k) (16)

Wherein H₁(k)＝DFT[h₁(τ)]；H_e(k) Channel estimation error for using FTCA method; g (l) represents the tap coefficients of the filter;indicating the number of taps of the FTCA filter,denotes rounding up, τ_maxRepresenting the maximum time delay, f_sRepresenting the sampling frequency, T_sRepresents a sampling interval;

s1.4, according to equations (15) and (17), the channel function H (k, m) is expressed as:

s1.5, representing the tap coefficient g (n, l) in equation (18) in the BEM form, then:

wherein the base coefficient g_q(l) A weighting coefficient representing the Q-th complex exponential base on the l-th path, Q being a natural number, andrepresenting the order of the basis extension model, f_maxRepresents the maximum doppler shift;

s1.6, substituting formula (20) into formula (18) to obtain:

when m-k + Q-Q/2 is an integer, H (k, m) ≠ 0, and 0. ltoreq. m, k. ltoreq. N-1, so equation (21) is simplified as:

given k (0. ltoreq. k < N) by equation (22), each Q (0. ltoreq. Q. ltoreq. Q) corresponds to m, so that H (k, m) is not zero, and the non-zero element position of the k-th row of the channel matrix corresponding to H (k, m) is defined as:

s1.7, obtaining a signal expression received by a channel model in the FTCA-CE-BEM method according to the formula (22), the formula (23) and the formula (13), wherein the signal expression is as follows:

as a preferred approach, the CIR of the multipath channel is expressed as:wherein h (l) is channel gain, T_sRepresents the sampling interval, i_lT_sThe channel delay of the l-th path is represented, and δ (n) represents an impulse function.

Preferably, in step S1.7, the base number g in formula (24)_q(l) The calculation steps are as follows:

A. the complex exponential basis coefficients are represented by a vector g of dimension (Q +1) mx 1, and are expressed as:

wherein, g_q＝[g_q(0),g_q(1),…,g_q(M-1)]^TQ is more than or equal to 0 and less than or equal to Q, and represents a coefficient vector with dimension of M multiplied by 1 corresponding to the qth basis function;

B. matrix with dimension N × NF_MRepresenting a matrix consisting of the first M columns of F, then F_MThe k-th row of expressions of (a) is:

C. the equation (24) is derived as:

D. using P pilot signals to estimate g, and their positions are k (1), k (2), …, k (P), respectively, according to equation (27), P linear equations can be expressed as:

E. defining a pilot sequence:

andis a P x 1 dimensional real number vector, Q ∈ {0,1, …, Q }, and equation (27) is changed to:

Y＝[diag(X₀)F₀,…,diag(X_Q)F_Q]g+W＝Ag+W (29)

wherein the content of the first and second substances,

A＝[diag(X₀)F₀,…,diag(X_Q)F_Q],Y＝[Y(k(1)),…,Y(k(P))]^T,W＝[W(k(1)),…,W(k(P))]^T，diag[·]representing a diagonal matrix, i.e. the elements not on the main diagonal are all zero, the elements on the main diagonal andq belongs to the element one-to-one correspondence in {0,1, …, Q };

F. from equation (29), a least squares estimate of g is obtained as:

g_LS＝(A^HA)^-1A^HY (30)。

preferably, in step S2, obtaining the received signal after ICI cancellation by using a zero-forcing equalization method, specifically including the following steps:

s2.1, obtaining a base coefficient least square estimation quantity g according to the formula (30) to substitute the formula (22) and obtaining a channel function H (k, m);

s2.2 according to zero-forcing equalisation methods, i.e. X ═ H^-1(k, m) Y, estimating an input signal X from the received signal Y, and substituting X into formula (14) to obtain an inter-subcarrier interference value I' (k);

s2.3, using FTCA-CE-BEM method, the received signal Y (k) is obtained from equation (13), and then the received signal Y' (k) with ICI removed is calculated as:

Y′(k)＝Y(k)-I′(k)＝H(k)X(k)+O(k) (31)

where o (k) ═ w (k) + I (k) -I' (k) denotes the remaining ICI and channel noise.

Preferably, in step S3, the estimate h of the CIR is calculated using a least squares method_LS(n), comprising the steps of:

s3.1, channel estimation according to the least-squares method, i.e.And according to formula (31) in step S2.3, the channel function expression is:

s3.2. for H_LS(k) Performing an inverse discrete fourier transform, namely:

to obtain:

h_LS(n)＝h(n)+o(n)，0≤n≤N-1 (33)

wherein o (n) ═ IDFT [ o (k)/x (k) ].

Preferably, step S4 includes: h 'are'_LS(n) is divided into three parts, n is more than or equal to 0 and less than or equal to L_cp-1 and N-L_cpN is not more than N and not more than N'_LS(n) as signal class training samples, apply L_cp≤n≤N-L_cp-1 part of h'_LS(n) as noise-like training samples, L_cpIndicating the cyclic prefix length.

Preferably, in step S4:

computing signal class initial clustering centersThe expression of (a) is:

computing noise-like initial cluster centersThe expression of (a) is:

preferably, step S5 includes the following steps:

s5.1, respectively calculating signal class h'_LS(n) and noise-like h'_LS(n) distances to the cluster centers of the respective classes,are respectively marked asAndn is more than or equal to 0 and less than or equal to N-1. The expression is as follows:

s5.2, setting a discriminant function r (n), wherein the expression is as follows:

s5.3, if r (n) is less than or equal to 0, corresponding h'_LS(n) is judged to be a signal type, and if r (n) > 0, corresponding h'_LSAnd (n) judging the cluster centers as noise classes, and recalculating the updated cluster centers respectively.

Preferably, step S7 specifically includes:

for h 'obtained in step S5'_k-means(n) performing a discrete fourier transform, namely:

obtaining a frequency domain channel function H 'after feedback'_k-means(k) I.e. the final channel estimation result, is:

H′_k-means(k)＝DFT[h′_k-means(n)] (39)。

the invention has the beneficial effects that: on the basis of an FTCA-CE-BEM method, ICI is eliminated through feedback, then a k-means method is introduced to reduce the influence of ICI and noise on communication, and density parameters are introduced to exclude isolated points in the k-means method, so that the estimation accuracy of the method is further improved.

Drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to the drawings without creative efforts.

FIG. 1 is a flow chart of an improved method of joint feedback k-means channel estimation approximating complex exponential basis expansion;

FIG. 2 is a flow chart of the FTCA-CE-BEM method for calculating the received signal after ICI cancellation;

FIG. 3 is a flow chart of the method of least squares to obtain CIR after removing isolated points;

FIG. 4 is a flow chart of the initial clustering center obtained by the modified k-means method;

FIG. 5 is a flow chart of recalculating the cluster centers using the discriminant function;

FIG. 6 is a flow chart of the iterative method to obtain the time-domain CIR;

FIG. 7 is a flow chart of obtaining a final frequency domain CIR;

FIG. 8 is an OFDM system model;

FIG. 9 is a non-sampling interval CIR energy distribution;

fig. 10 is a graph of the relationship between the signal-to-noise ratio and the bit error rate of the four algorithms.

Detailed Description

The following description of the embodiments of the present invention is provided by way of specific examples, and other advantages and effects of the present invention will be readily apparent to those skilled in the art from the disclosure herein. The invention is capable of other and different embodiments and of being practiced or of being carried out in various ways, and its several details are capable of modification in various respects, all without departing from the spirit and scope of the present invention. It is to be noted that the features in the following embodiments and examples may be combined with each other without conflict.

An improved joint feedback k-means channel estimation method of approximate complex exponential basis expansion is taught below, and is realized based on a Basis Expansion Model (BEM) method, a Fractional Tap Channel Approximation (FTCA) method, a k-means (k-means) method, a least square method (LS) method, a Discrete Fourier Transform (DFT) method, Gaussian distribution, a zero forcing equalization method, a Discrete Fourier Transform (DFT) and an Inverse Discrete Fourier Transform (IDFT). Wherein:

basis Extension Model (BEM) method:

the base-extended model (BEM) is a fast time-varying channel estimation method in wireless communications. The essence is to approximate a fast time-varying channel with a weighted superposition of a small number of basis functions, i.e. tap coefficients g (n, l) with basis function b_q(n) and base coefficient g_q(l) The process is represented as:wherein the base coefficient g_q(l) A weighting coefficient indicating the q-th complex exponential base on the l-th path; q is a natural number, andrepresenting the order of the basis expansion model; f. of_maxIndicating the maximum doppler shift.

If the basis functions are expressed as Fourier functions, i.e. b_q(n)＝e^{j2π(q-Q/2)n/N}N is more than or equal to 0 and less than or equal to N-1, Q is more than or equal to 0 and less than or equal to Q, the method is called a complex exponential basis expansion model (CE-BEM), wherein e is a natural constant 2.718128.

Fractional Tap Channel Approximation (FTCA) method:

the Fractional Tap Channel Approximation (FTCA) method is a fast time varying channel estimation method in wireless communications. The essence is that a fractional weighting factor K is introduced_α(0＜K_α1) to simulate non-sampling interval channels, i.e. filters h₁The tap interval of (τ) is a fraction of the sampling interval. The actual channel frequency domain function is then: h_FTCA(k)＝H₁(k)+H_e(k) In which H is₁(k)＝DFT[h₁(τ)]And is andH_e(k) channel estimation error for using FTCA method; g (l) represents the tap coefficients of the filter;indicating the number of taps of the FTCA filter,denotes rounding up, τ_maxRepresenting the maximum time delay, f_sRepresenting the sampling frequency.

k-means (k-means) method:

the k-means (k-means) method is an iterative solution clustering method in data analysis. The essence is that the distance of the data sample to the cluster center is minimal. The method comprises the following implementation steps: 1. selecting k data samples as initial clustering centers; 2. assigning each data sample to the nearest cluster center, forming k clusters; 3. recalculating the cluster center of each cluster; 4. and stopping iteration until the cluster does not send changes or the maximum iteration times are reached, and obtaining a final clustering result. Otherwise, the step 2 is executed again.

Least Squares (LS) method:

in a fast time-varying Orthogonal Frequency Division Multiplexing (OFDM) system, the received signal y (k) has the expression: y (k) ═ H (k, k) x (k) + i (k) + w (k), where H (k, k) is the desired channel function, w (k) is the frequency domain channel additive noise, and i (k) represents the sub-carrier interference (ICI) caused by the fast time-varying channel. The LS channel estimation method estimates parameters H (k, k) to minimize an objective function J, which is expressed as:

J＝[Y(k)-H(k,k)X(k)]^T[Y(k)-H(k,k)X(k)] (1)

let J solve the partial derivative of H (k, k), the expression is:

let J pair H (k, k)) And (5) solving a second-order partial derivative:since the second order partial derivative is greater than zero, there is a minimum value for J. Order toThe expression of the channel function under the available LS channel estimation method is:

discrete Fourier Transform (DFT) method:

the DFT channel estimation method is realized by firstly obtaining a channel function at a pilot frequency position by using the least square method (LS) channel estimation method, wherein the expression is as follows:

where H (k, k) is the desired channel function, W₁(k) I (k) + w (k). To pairPerforming a discrete fourier transform (IDFT), the expression:

wherein H (n) ═ IDFT [ H (k, k)]，w₁(n)＝IDFT[W₁(k)/X(k)]And n is a natural number and represents a subscript.

Distributed mainly on both sides according to the non-sampling interval CIR and with a cyclic prefix length L_cpGenerally, the length of the CIR is not less than the length of the CIR, and the time domain channel response after noise elimination and zero padding can be written, and the expression is:

to pairPerforming Discrete Fourier Transform (DFT) to obtain the channel frequency domain response under the DFT channel estimation method, wherein the expression is as follows:wherein N and k represent the following table, N represents the number of OFDM subcarriers, and N, k and N are all natural numbers.

Gaussian distribution:

the gaussian distribution is also called normal distribution, if a random variable X obeys a location parameter of μ, a scale parameter of σ, and its probability density function is:x is called a normal random variable, and the distribution to which X obeys is called a normal distribution and can be recorded as X-N (mu, sigma)²). When μ is 0 and σ is 1, the distribution to which X obeys is referred to as a standard normal distribution, denoted as X to N (0,1), where —, denotes obedience, and N denotes a normal distribution.

The zero forcing equalization method comprises the following steps:

the zero-forcing equalization method can obtain an equalization estimation value of a transmission data vector X by using a peak distortion criterion as follows: x ═ H (H)^{- 1}Y, where Y represents a received data vector, H represents a channel response matrix, and the peak distortion criterion is defined as:where y (0) denotes the sample value at the time when t is 0, y_kThe tap coefficient determined by the formula can obtain the best equalization effect by representing the intersymbol interference at the sampling time of k, wherein k is a natural number.

Discrete Fourier Transform (DFT):

DFT can transform signal from time domain to frequency domain, and then study signal spectrum structure and change rule. For a sequence of N points { x (N) },0 ≦ N ≦ N-1, its DFT may be expressed as:wherein n and k are both natural numbers and both represent subscripts.

Inverse Discrete Fourier Transform (IDFT):

the IDFT may transform the signal from the frequency domain to the time domain. For a sequence of N points { x (k) },0 ≦ k ≦ N-1, its IDFT may be expressed as:wherein n and k are both natural numbers and both represent subscripts.

Referring to fig. 1, the present embodiment provides an improved method for estimating a k-means channel by joint feedback of approximate complex exponential base spreading, including the steps of:

s1, aiming at data sent by a transmitter, obtaining a receiving signal under a channel model by adopting a fractional tap channel approximate complex exponential basis expansion model method, wherein the fractional tap channel approximate complex exponential basis expansion model is FTCA-CE-BEM;

s2, obtaining a received signal with ICI eliminated by using a zero-forcing equalization method, wherein ICI is sub-carrier interference;

s3, aiming at the received signal obtained in the step S2 after ICI elimination, calculating the estimation h of the CIR by adopting a least square method_LS(n), simultaneously, deleting the isolated CIR by adopting a density parameter to obtain a non-isolated CIR: h'_LS(n), wherein the density parameter is the number of CIRs within a distance r from a certain CIR, the CIR is the channel impulse response, and r is the radius of the selected sphere;

s4, h 'obtained in step S3'_LS(n) dividing the signal into a noise class and a signal class, and calculating initial clustering centers of the respective classes;

s5, setting a discrimination function of the distance from the CIR to the clustering center to judge and reclassify all the CIRs and calculate the clustering center;

s6, judging whether the clustering result is changed or not, if so, returning to the step S5 according to the changed result, and if not, judging the clustering result to be h 'of the noise class'_LS(n) set to zero, and h 'determined to be signal-based by bonding'_LS(n) obtaining a time domain channel function h'_k-means(n)；

S7. pairh′_k-means(n) performing a DFT calculation to obtain a fed back frequency domain channel function H'_k-means(k) I.e., the final channel estimation result, DFT is a discrete fourier transform.

Specifically, the method comprises the following steps:

referring to fig. 2, the specific steps of calculating the received signal after ICI cancellation are as follows:

in step S1, obtaining the received signal under the channel model by using the FTCA-CE-BEM method, specifically including the following steps:

s1.1, transmit data X ═ X (0), X (1), …, X (N-1) to the transmitter]^TPerforming an N-point inverse discrete fourier transform to obtain a time domain signal;

namely, it isThe resulting time domain signal is represented as:

wherein N and N are natural numbers and respectively represent the serial number of elements in the x vector and the number of OFDM subcarriers; x ═ x (0), x (1), …, x (N-1)]^TThe real number vector of Nx 1 dimension represents the signal data vector of N points of each frame; f₁∈R^N×NA real number matrix representing dimensions N × N; w ═ e^-j2π/NAnd e is a natural constant 2.718128. Moreover, the vector x can also be expressed as:

wherein, F₂For normalizing Discrete Fourier Transform (DFT) matrices, i.e. F₂The sum of squares of each column (row) element is 1, and its (u, v) -th element can be represented as:

s1.2, y (n) is time domain data received by a receiving end. The product corresponding to the time domain is convolved by the convolution theorem, i.e. the frequency domain. y (n) can be expressed as:

wherein h (n, l) is a real number and is expressed as a sampling value of a fast time-varying Channel Impulse Response (CIR) at the nth time and the l path; l is the total path number of information transmission; x (n-l) is the input at time n-l; w (n) is the mean value of zero and the variance of σ²White gaussian noise. I.e. it satisfies N (0, σ)²) A gaussian distribution.

The CIR of a multipath channel may be expressed as:wherein h (l) is channel gain, T_sRepresents the sampling interval, i_lT_sThe channel delay of the l-th path is represented, and δ (n) represents an impulse function. If i_lA positive integer, the multipath channel is a sampling interval channel. At this time, the time domain channel matrix H₁From a vector of real numbers of dimension Nx 1Cyclic shift composition, where N represents the number of OFDM subcarriers, H₁The expression is as follows:

if i_lAnd is not an integer greater than 0, the multipath channel is a non-sampling interval channel. At this time, the time domain channel matrix H₂It will no longer be a column cyclic property matrix. At this time, H₂Can be expressed as:

DFT processing is performed on y (n), i.e.A frequency domain signal Y is obtained and can be expressed as:

wherein the content of the first and second substances,Y＝[Y(0),Y(1),…,Y(N-1)]^Tis a real number vector with dimension Nx 1; [. the]^HRepresenting a matrix transposition; w ═ W (0), W (1), …, W (N-1)]^TIs a real number vector of dimension N × 1.

Taking Y (n) obtained from equation (9) into Y obtained from equation (12), the k-th subcarrier reception signal is:

wherein H (k, k) is the desired channel function; w (k) is frequency domain channel additive noise; i (k) denotes ICI due to fast time varying channels, expressed as:

wherein, the expression of the channel function H (k, m) is:

s1.3, adopting a fractional tap channel approximation method, and introducing a fractional weighting factor K_α(0＜K_α1) to simulate a non-sampling interval channel, and K_αAre real numbers. I.e. the filter h₁Tap spacing K of (tau)_αT_sFor a sampling interval T_sMultiple of the fraction of. The actual channel frequency domain function H_FTCA(k) Expressed as:

H_FTCA(k)＝H₁(k)+H_e(k) (16)

wherein H₁(k)＝DFT[h₁(τ)]；H_e(k) Channel estimation error for using FTCA method; g (l) represents the tap coefficients of the filter;indicating the number of taps of the FTCA filter,denotes rounding up, τ_maxRepresenting the maximum time delay, f_sRepresenting the sampling frequency.

S1.4, according to equations (15) and (17), the channel function H (k, m) is expressed as:

the Basis Extension Model (BEM) can approximate a fast time-varying channel with a weighted superposition of a small number of basis functions. That is, the tap coefficient g (n, l) can be used as the basis function b_q(n) and base coefficient g_q(l) Expressed as:

wherein the base coefficient g_q(l) A weighting coefficient indicating the q-th complex exponential base on the l-th path; q is a natural number, andrepresenting the order of the basis expansion model; f. of_maxIndicating the maximum doppler shift.

S1.5, representing the tap coefficient g (n, l) in equation (18) in the BEM form, then:

s1.6, substituting formula (20) into formula (18) to obtain:

when m-k + Q-Q/2 is an integer, H (k, m) ≠ 0, and 0. ltoreq. m, k. ltoreq. N-1, so equation (21) is simplified as:

given k (0. ltoreq. k < N) by equation (22), each Q (0. ltoreq. Q. ltoreq. Q) corresponds to m, so that H (k, m) is not zero, and the non-zero element position of the k-th row of the channel matrix corresponding to H (k, m) is defined as:

s1.7, obtaining a signal expression received by a channel model in the FTCA-CE-BEM method according to the formula (22), the formula (23) and the formula (13), wherein the signal expression is as follows:

the expression of the signal received by the FTCA-CE-BEM method is obtained from the expression (24) in step 1.7. Wherein, the base coefficient g in the OFDM symbol period_q(l) The pilot frequency information needs to be inserted for estimation, and the specific calculation steps are as follows:

A. to simplify the calculation, the complex exponential basis coefficients are represented by a vector g with dimension (Q +1) mx 1, and are expressed as:

wherein, g_q＝[g_q(0),g_q(1),…,g_q(M-1)]^TQ is more than or equal to 0 and less than or equal to Q, and represents a coefficient vector with dimension of M multiplied by 1 corresponding to the qth basis function;

B. matrix with dimension N × NF_MRepresenting a matrix consisting of the first M columns of F, then F_MThe k-th row of expressions of (a) is:

C. the equation (24) is derived as:

D. as can be seen from equation (27), in the fast time varying channel, the received signal y (k) is affected not only by the k-th subcarrier but also by Q subcarriers adjacent to the k-th subcarrier. Therefore, if x (k) is a pilot signal, then the transmit signals affecting the Q subcarriers of the received signal y (k) should also be pilot signals in order to estimate g.

In a fast time varying system, P pilot signals are used to estimate g, and their positions are k (1), k (2), …, k (P), respectively. From equation (27), the P linear equations can be expressed as:

as can be seen from equation (28), calculating g requires knowing the associated transmitted signals for the P pilot locations that affect signal reception:q ∈ {0,1, … Q }. Therefore, it is estimated that g needs to transmit at least P × (Q +1) pilots.

E. Defining a pilot sequence:

andis a P x 1 dimensional real number vector, Q ∈ {0,1, …, Q }, and equation (27) is changed to:

Y＝[diag(X₀)F₀,…,diag(X_Q)F_Q]g+W＝Ag+W (29)

wherein the content of the first and second substances,

A＝[diag(X₀)F₀,…,diag(X_Q)F_Q],Y＝[Y(k(1)),…,Y(k(P))]^T,W＝[W(k(1)),…,W(k(P))]^T，diag[·]representing a diagonal matrix, i.e. the elements not on the main diagonal are all zero, the elements on the main diagonal andq belongs to the element one-to-one correspondence in {0,1, …, Q };

F. from equation (29), a least squares estimate of g is obtained as:

g_LS＝(A^HA)^-1A^HY (30)。

in step S2, obtaining the received signal after ICI cancellation by using a zero-forcing equalization method, specifically including the following steps:

s2.1, obtaining a base coefficient least square estimation quantity g according to the formula (30) to substitute the formula (22) and obtaining a channel function H (k, m);

s2.2 according to zero-forcing equalisation methods, i.e. X ═ H^-1(k, m) Y, estimating an input signal X from the received signal Y, and substituting X into formula (14) to obtain an inter-subcarrier interference value I' (k);

s2.3, using FTCA-CE-BEM method, the received signal Y (k) is obtained from equation (13), and then the received signal Y' (k) with ICI removed is calculated as:

Y′(k)＝Y(k)-I′(k)＝H(k)X(k)+O(k) (31)

where o (k) ═ w (k) + I (k) -I' (k) denotes the remaining ICI and channel noise.

Referring to fig. 3, the CIR after the outliers are deleted is obtained using the least square method. The method is mainly completed by the following steps:

in step S3, the least square method is used to calculate the estimated CIR h_LS(n), comprising the steps of:

s3.1, channel estimation according to the least-squares method, i.e.And according to formula (31) in step S2.3, the channel function expression is:

s3.2. for H_LS(k) Performing an inverse discrete fourier transform, namely:

to obtain:

h_LS(n)＝h(n)+o(n)，0≤n≤N-1 (33)

wherein o (n) ═ IDFT [ o (k)/x (k) ].

Isolated point definition: that is, the distance from a certain CIR to the corresponding cluster center is greater than the average distance from other CIRs in the cluster to the cluster center, and the CIRs with sparse distribution are called isolated points. Therefore, the density parameter is introduced to describe the sparsity of CIR distribution. The density parameter is defined as the number of CIRs within a distance r from the CIR. Where r is the radius of a given sphere. Then, the density parameter S for each CIR can be calculated_nAnd compared to a given minimum density Q. If S_nIf Q is less than or equal to Q, the corresponding CIR distribution is considered to be sparse, and the CIR meeting the definition of the isolated point is set to zero. The result after the treatment was recorded as h'_LS(n)。

Referring to fig. 4, an initial cluster center is obtained using the modified k-means method. The method is mainly completed by the following steps:

in step S4, the k-means method may classify the data according to the similarity between the data objects. The groups with larger similarity are divided into one group, and the number of the classification groups is the k value. If k is set to 2, i.e.Received signals are classified into noise and signal classes. Distributed mainly on both sides according to the non-sampling interval CIR and with a cyclic prefix length L_cpGenerally not less than CIR length, h'_LS(n) is divided into three parts. N is more than or equal to 0 and less than or equal to L_cp-1 and N-L_cpN is not more than N and not more than N'_LS(n) as signal class training samples.

In step S4:

computing signal class initial clustering centersThe expression of (a) is:

computing noise-like initial cluster centersThe expression of (a) is:

referring to fig. 5, in order to use the discriminant function, a decision result is obtained and the clustering center is recalculated. The method is mainly completed by the following steps:

in step S5, the method includes the steps of:

s5.1, respectively calculating signal class h'_LS(n) and noise-like h'_LS(n) distances to the various cluster centers, respectivelyAndn is more than or equal to 0 and less than or equal to N-1. The expression is as follows:

s5.2, setting a discriminant function r (n), wherein the expression is as follows:

s5.3, if r (n) is less than or equal to 0, corresponding h'_LS(n) is judged to be a signal type, and if r (n) > 0, corresponding h'_LSAnd (n) judging the cluster centers as noise classes, and recalculating the updated cluster centers respectively.

Referring to fig. 6, the proposed method time-domain CIR is obtained by using an iterative method. The method is mainly completed by the following steps:

s6, judging whether the clustering result is changed or not, if so, returning to the step S5 according to the changed result, and if not, judging the clustering result to be h 'of the noise class'_LS(n) set to zero, and h 'determined to be signal-based by bonding'_LS(n) obtaining a time domain channel function h'_k-means(n)；

Referring to FIG. 7, is pair h'_k-means(n) performing DFT to obtain a final frequency domain CIR. The method is mainly completed by the following steps:

in step S7, the method specifically includes:

for h 'obtained in step S5'_k-means(n) performing a discrete fourier transform, namely:

obtaining a frequency domain channel function H 'after feedback'_k-means(k) I.e. the final channel estimation result, is:

H′_k-means(k)＝DFT[h′_k-means(n)] (39)。

referring to fig. 8, which is a schematic diagram of a system model in this embodiment, it can be seen that an analog signal sent by a transmitter is QPSK modulated, then is serial-to-parallel converted, and pilot is inserted to form M parallel data paths. After the IDFT is performed on these data, N time-domain discrete signals are obtained, and then a Cyclic Prefix (CP) is inserted, for example, assuming that the CP length is 2, the last three bits of N ═ 6 time-domain discrete signals {1,2,3,4,5,6} are taken as the CP, so that the CP-inserted signals become {5,6,1,2,3,4,5,6 }. The purpose of inserting the CP is to make the digital signal cyclic, so as to remove the inter-symbol interference. Then the digital signal after inserting CP is changed by parallel-serial, and the converted digital signal is converted into analog signal and sent out. The transmitted signal passes through a fast time varying channel and is also affected by white gaussian noise. After receiving the transmitted signal, the receiver converts the analog signal into a digital signal, then carries out serial-to-parallel conversion, and then removes the CP according to the method of inserting the CP. The CP time is required to be larger than the maximum delay of the channel, otherwise the inter-symbol interference cannot be completely eliminated. The DFT is performed on the signal without the CP to obtain a frequency domain signal, where the received signal may be represented as Y ═ HX + W, where Y is the received signal, X is the transmitted signal, H is the frequency domain response function of the channel, and W is gaussian white noise. The received signal is subjected to frequency domain equalization, then to parallel-to-serial conversion, and finally to QPSK demodulation, so as to obtain the original transmitted signal.

Referring to fig. 9, which is a schematic diagram of the energy distribution of the non-sampling interval CIR in this embodiment, a theoretical expression selected by simulation is as follows: h (T) δ (T-2.5T)_s)+0.8δ(t-6.5T_s) Sampling interval f_s1 μ s, the number of multipaths L is 5, as shown in fig. 9: the non-sampling interval CIR energy is mainly concentrated on both sides. Using this property, the proposed method will be h 'after removal of ICI and isolated points'_LS(n) is divided into three parts, and n is more than or equal to 0 and less than or equal to L_cp-1 and N-L_cpN is less than or equal to N is used as a signal, L is_cp≤n≤N-L_cpPart-1 is referred to as noise.

Fig. 10 is a graph showing a comparison of the relationship between the signal-to-noise ratio and the bit error rate at a moving speed v of 300Km/h according to the channel estimation method using FTCA-CE-BEM-k-MEANS, the channel estimation method using FTCA-CE-BEM, the channel estimation method using CE-BEM, and the channel estimation method using BEM provided in this embodiment. The data used for the simulation are shown in table 1 below:

TABLE 1 simulation data sheet

Parameter(s)	Numerical value
		Carrier frequency f0	5GHz
Sampling interval Ts	5.75μs
		Velocity v of mobile station	300Km/h
Number of subcarriers N	128
		Modulation system	QPSK
Number of multipaths L	5
		Normalized time delay ND	4
Doppler maximum frequency offset fmax	1388Hz
		Cyclic prefix length Lcp	32

From the table, the base extension orderAccording to a weighting factor K_αAnd normalized maximum time delay N_DThe value range of the weighting factor is K_α∈[0.36,0.57]. K is obtained in simulation experiment of the invention_α0.49. The number of taps of the filter approximated by the algorithmic modelWherein f is_s＝1/T_sThe total number of pilot frequencies required for simulation is P_totalThe pilot number of each subsequence is 32, and the corresponding set of positions of each Q-Q/2 pilot subsequence is K_Q/22,10,18,26,34, …,250, which is used for channel estimation, other pilot positions can be determined.

All basis coefficients g when generating a channel using a basis extension model_q(l) Obey a complex gaussian random distribution with variance:wherein the content of the first and second substances,distribution function, tau, representing the multipath intensity of the channel_rms＝4T_sThe mean-squared delay of the path is denoted 23 μ s, and the doppler power spectral density is:

referring to fig. 10, it can be seen that the error rate decreases greatly with increasing signal-to-noise ratio in the FTCA-CE-BEM and FTCA-CE-BEM-k-MEANS methods. The reason is that: the BEM method is based on estimating a small number of basis coefficients g_q(l) By calculated tap coefficientsThe fast time-varying channel is simulated, but energy leakage is generated, so the estimation effect of the method is the worst, and the system can hardly work normally. The CE-BEM method adopts a Fourier function b on the basis of retaining the advantages of the BEM method_q(n)＝e^{j2π(q-Q/2)n/N}As a basis function. Thus, the method is simple in structure and easy to implement, but b_qThe introduction of (n) can create spectral leakage problems. Therefore, the process performance is only slightly better than the BEM process. The FTCA-CE-BEM method introduces a weighting factor K after inheriting the advantages of the CE-BEM method_α(0＜K_αLess than or equal to 1) to estimate not only the sampling interval channel but also the non-sampling interval channel, so as to obtain a better estimation fast time-varying channel. Therefore, the method performance is greatly improved compared with CE-BEM. Based on the original algorithm, the proposed FTCA-CE-BEM-k-MEANS method utilizes the formula (99) to eliminate ICI and introduce the density parameter S through feedback_nThe improved k-means method further eliminates the influence of ICI and noise on the received signal by using an iterative idea, effectively reduces CIR energy leakage and further improves the system performance.

The above-mentioned embodiments are merely illustrative of the preferred embodiments of the present invention, and do not limit the scope of the present invention, and various modifications and improvements of the technical solutions of the present invention by those skilled in the art should fall within the protection scope of the present invention without departing from the design spirit of the present invention.