Non-stationary channel simulator design method based on subspace reconstruction

Abstract

The invention discloses a method for designing a non-stationary channel simulator based on subspace reconstruction. The method comprises the following steps: 1) configuring basic parameters such as an environment and an antenna in a channel model; 2) generating a channel transmission function by using a three-dimensional space-time frequency non-stationary geometric random channel model; 3) projecting a channel transmission function to a subspace expanded by a group of Chirp signals, and sending and storing a projection coefficient and a Chirp signal parameter to the FPGA; 4) realizing low-complexity channel transmission function reconstruction on the FPGA according to the parameters; and 5) performing frequency domain processing on the sending signal by the reconstructed channel transmission function to obtain a receiving signal after channel fading is applied. According to the non-stationary channel simulator design method established by the invention, simulation of a space-time frequency non-stationary MIMO channel is supported, the blank in the field of channel simulators is filled up, and the proposed hardware architecture has important reference value for the design of the non-stationary channel simulator based on a geometric random channel model.

Description

Technical Field

[0001] The invention belongs to the field of channel modeling and channel simulators, and in particular relates to a non-stationary channel simulator design method based on subspace reconstruction. Background Art

[0002] Before a new communication system is commercialized, its performance needs to be verified by actual measurements in the corresponding environment. Although it is possible to conduct field tests on the system at different locations, this method is costly and easily affected by the surrounding environment. In addition, the comparative analysis of simulation results makes it difficult to reproduce the propagation conditions of the real field. A more practical approach is to create a suitable and stable simulation environment for the channel and then evaluate the performance of the communication system in these environments. Channel simulators can simulate channel environments in a controllable and repeatable manner and can be used for consistency testing, performance testing, and interoperability testing of communication systems. This means that field testing is not required and time and cost efficiency can be greatly improved. At present, due to the high resource consumption of multiple-input multiple-output (MIMO) channel simulators based on tapped-delay-line (TDL) architecture, frequency domain processing channel simulators suitable for high-order MIMO do not yet support non-stationary geometric random channel models with high computational complexity, and lack pre-processing algorithms for non-stationary geometric random channel models to reduce their computational complexity, so as to apply them to large-scale MIMO channel simulators based on frequency domain processing.

[0003] In summary, it is very necessary to establish a low-complexity MIMO channel simulator using a non-stationary geometric random channel model. Summary of the invention

[0004] The purpose of the present invention is to provide a non-stationary channel simulator design method based on subspace reconstruction, which uses a geometric random channel model to perform high-precision MIMO channel simulation to accurately and stably test large-scale MIMO communication systems in mobile scenarios, thereby solving the problem that the channel simulator based on the geometric random channel model needs to calculate a large number of time-varying subpath parameters.

[0005] In order to solve the above technical problems, the specific technical solutions of the present invention are as follows:

[0006] A method for designing a non-stationary channel simulator based on subspace reconstruction, characterized in that it comprises the following steps:

[0007] Step S1: Based on the universal geometric random channel model, the environment and basic parameters of the transmitting and receiving antennas in the channel model are configured on the host computer;

[0008] Step S2, generating a three-dimensional space-time-frequency non-stationary channel transfer function based on a universal geometric random channel model;

[0009] Step S3: Project the channel transfer function onto the subspace spanned by the linear frequency modulation Chirp signal to obtain a set of Chirp signal basis parameters;

[0010] Step S4, sending and storing the basic parameters obtained in step S3 to a field programmable gate array FPGA, and obtaining a reconstructed channel transfer function by superimposing the Chirp signal on the FPGA;

[0011] Step S5, performing frequency domain processing on the MIMO channel simulator input signal according to the channel transfer function reconstructed in step S4;

[0012] Step S6: verify the consistency of the channel transfer function, Doppler power spectral density and delay power spectral density of the reconstructed channel and the original channel.

[0013] 2. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 1, wherein step S1 specifically comprises:

[0014] The environment and basic parameters of the transceiver antenna configured on the host computer include simulation time, time sampling interval, carrier frequency, channel bandwidth, frequency sampling interval, initial position coordinates of the transceiver, moving speed of the transceiver, center position coordinates of the dual clusters of the transceiver, moving speed of the dual clusters of the transceiver, birth and death rate of the dual clusters of the transceiver, number of antennas and antenna interval of the transceiver, azimuth angle and elevation angle of the antennas of the transceiver.

[0015] 3. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 1, wherein step S2 specifically comprises:

[0016] Step S201, calculate the initial position coordinates of the pth antenna unit at the transmitting end:

[0017]

[0018] in, is the initial position coordinate of the first antenna unit at the transmitting end, δ T is the transmitting antenna spacing, and are the azimuth and elevation angles of the transmitting antenna; calculate the coordinates of the qth antenna unit at the receiving end;

[0019]

[0020] in, is the initial position coordinate of the first antenna unit at the receiving end, δ R is the receiving antenna spacing, and is the azimuth and elevation angle of the receiving antenna;

[0021] Step S202: Generate scatterer position coordinates. The scatterer position coordinates (x′, y′, z′) are distributed according to the ellipsoidal Gaussian distribution with the center of the dual cluster of the transmitting and receiving ends as the Cartesian coordinate origin. The joint probability density of the scatterer position coordinates (x′, y′, z′) is:

[0022]

[0023] Among them, σ AS represents the cluster azimuthal spread, σ ES represents the cluster pitch angle spread, σ DS represents the cluster delay spread; the Cartesian coordinates of the mth scatterer in the nth transmitting cluster relative to the first transmitting antenna element The coordinate transformation is obtained by the following formula:

[0024]

[0025] in, are the spherical coordinates of the centers of the n originating clusters;

[0026] Cartesian coordinates of the mth scatterer in the nth receiving cluster relative to the first receiving antenna element The coordinate transformation is obtained by the following formula:

[0027]

[0028] in, are the spherical coordinates of the centers of the n clusters at the receiving end;

[0029] Step S203, calculate the multipath propagation distance and delay between the antenna units at the transmitting and receiving ends, and obtain the distance from the p-th transmitting antenna unit to the m-th scatterer in the n-th cluster at time t based on the initial coordinates and moving speed of the antenna units at the transmitting and receiving ends, the initial coordinates of the scatterers in the dual clusters, and the moving speed of the clusters. The calculation method is expressed as:

[0030]

[0031] in, and are the moving speed of the transmitting cluster and the moving speed of the transmitting antenna respectively; the distance from the qth receiving antenna unit to the mth scatterer in the nth cluster at time t is obtained, and the calculation method is expressed as:

[0032]

[0033] in, and are the moving speed of the receiving cluster and the moving speed of the receiving antenna respectively; the sub-path delay caused by multipath propagation between the antenna units at the transmitting and receiving ends is calculated as follows:

[0034]

[0035] Where c is the speed of light, is the distance between the two clusters at the transmitting and receiving ends, τ link is the virtual link delay that obeys exponential distribution;

[0036] Step S204: Generate the power of each sub-path according to the sub-path delay obtained in step S203. The calculation method is expressed as:

[0037]

[0038] Among them, f c is the carrier frequency, DS is the delay spread, r τ and γ are factors that depend on the time delay and frequency domain, respectively;

[0039] Step S205: Generate a channel transfer function according to the sub-path delay and sub-path power obtained in step S203 and step S204, and the calculation method is:

[0040]

[0041] Among them, f c is the carrier frequency, N qp (t) and M n (t) are the number of clusters and the number of scatterers in the cluster respectively; according to the simulation time, channel bandwidth, antenna array size, and based on the birth and death rate of the dual clusters at the transmitting and receiving ends in the space-time-frequency domain, the index of the surviving cluster is determined, and according to step S202, the position information of the newly generated cluster that changes with the antenna unit interval, time interval, and frequency interval is calculated, and then according to steps S203 and S204, the sub-path delay and sub-path power of the corresponding cluster are calculated, and then according to step S205, the non-stationary channel transfer function updated with the space-time-frequency evolution is calculated.

[0042] 4. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 3, characterized in that the step S3 specifically comprises:

[0043] Step S301: When the transceiver and the dual clusters move at a fixed speed and direction, the distance calculation in step S203 is approximately:

[0044]

[0045] in, is the fixed speed of the originating cluster, is the fixed speed of the transmitting antenna, is the initial distance from the first antenna unit at the transmitting end to the mth scatterer in the nth cluster at the transmitting end, and is a parameter that depends on the antenna array, the moving speed and the angle of the scatterer;

[0046]

[0047] in, The fixed speed of the receiving cluster is is the fixed speed of the receiving antenna, is the initial distance from the first antenna unit at the receiving end to the mth scatterer in the nth cluster at the receiving end, and is a parameter that depends on the antenna array, the moving speed and the angle of the scatterer;

[0048] Step S302: Substituting the approximate calculation relationship in step S301 into step S204, the channel transfer function of the MIMO channel is approximately calculated as:

[0049]

[0050] Where L(t) = max{N qp (t)·M n (t)},η qp，l ,ξ qp，l , t qp，l , a qp，l , b qp，l are all constant factors; j is the imaginary unit; l is the summation subscript index; the approximate calculation method of the Doppler frequency of the channel model is expressed as:

[0051]

[0052] dt is the differential symbol, T and R are the superscripts of the transmitting and receiving ends respectively.

[0053] Step S303: According to the approximate relationship in step S302, the Doppler frequency of the channel model is approximated by a linear time-varying function, and the channel transfer function is obtained by a linear combination of multiple Chirp signals; thus, a K-dimensional subspace is constructed using the Chirp signal as a basis function, K<<L(t); the k-th basis function of the K-dimensional subspace is expressed as:

[0054]

[0055] Among them, α k represents the initial frequency, β krepresents the Chirp rate; according to the approximate relationship in step S302, considering the case where the transmitting end is stationary and the receiving end is moving, α k and β k The value range of is expressed as:

[0056]

[0057]

[0058] Thus α k and β k The value of is obtained from the statistical distribution of the coefficient corresponding to the approximately linear time-varying Doppler frequency in the channel model;

[0059] Step S304, determine the projection coefficient of the channel transfer function on the basis function in step S303; consider any subchannel in the MIMO channel, and express the channel fading at different frequencies as:

[0060] h=[H(t,f1),H(t,f2),…,H(t,f I )]

[0061] Among them, H(t, f i ), i = 1, 2, ..., I is the H(t, f) of the MIMO subchannel at frequency f i The sample value on ;

[0062] The basis function of the subspace in step S303 is expressed as:

[0063] Ψ=[b1(t),b2(t),…,b K (t)]

[0064] Among them, b k (t), k = 1, 2, ..., K is the basis function in step S303;

[0065] Assume that Ψ is orthogonalized from matrix G to matrix Φ, that is, Φ = Ψ·G, where G is expressed as:

[0066] in,

[0067] <,> represents the inner product operation, b k (t) is the kth column vector of Ψ, g j (t) is the jth column vector of Φ; the projection X of h on Φ is expressed as:

[0068]

[0069] x ji is the jth row and ith column element of X; thus, the projection A of the channel fading h on the Chirp basis function Ψ is expressed as:

[0070] A=G·X.

[0071] 4. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 3, characterized in that the step S4 specifically comprises:

[0072] The channel matrix preprocessing results based on subspace projection in step S303 and step S304 are Send and store to FPGA, so that this set of basic parameters in the time window T w The channel transfer function is reconstructed in , and the calculation method is expressed as:

[0073]

[0074] in, is the reconstructed channel transfer function, that is

[0075]

[0076] Among them, a k (f i ) is the kth row and ith column element of A, t0 is the initial time of channel simulation, T w is the duration of the channel simulation; therefore, a low-complexity channel transfer function calculation is completed on the FPGA by superposition of time-division Chirp signals.

[0077] 5. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 4, characterized in that step S5 specifically comprises:

[0078] Step S501: The channel transfer function calculated in step S4 Get the channel time sampling interval T ch The channel frequency response within H qp (f i );According to the predefined channel bandwidth and frequency sampling interval in the channel model, the actual channel frequency response is a discrete H qp (m);

[0079] Step S502: From the input signal stream s p Take N from (n) s point signals p，i (n) put into cache, N s The calculation method is expressed as:

[0080] N s =N H -N a

[0081] Among them, NH is the number of channel frequency response points, means round down, τ m =max{τ|h qp (t, τ) ≥ 0}, h(t, τ) represents the channel impulse response, and is the channel transfer function H qp The inverse Fourier transform of (t, f) with respect to f;

[0082] Step S503: In the cache N s The input signal of point s p，i (n)N a Zeros and perform fast Fourier transform FFT to get N points H The input signal spectrum S p，i (m); then the q-th output signal spectrum of the P×QMIMO channel simulator is expressed as:

[0083]

[0084] Step S504: q，i (m)N H The inverse fast Fourier transform IFFT of the point is used to obtain the Q-channel output signal r q，i (n); Using the overlap-add method, the qth continuous output signal stream of the channel simulator is expressed as

[0085]

[0086] Step S505: According to step S4, the channel frequency response H is refreshed. qp (m), and repeat steps S501-S504; channel update time T ch , that is, the channel time sampling interval, the calculation method is expressed as:

[0087] T ch =N s ·T s

[0088] Among them, T s is the input signal sampling period.

[0089] 6. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 5, characterized in that step S6 specifically comprises:

[0090] Step S601: The average error between the channel transfer function based on the subspace reconstruction and the channel transfer function based on the geometric random channel model is calculated as follows:

[0091]

[0092] Where B represents the channel bandwidth, H(t, f) and They represent the channel transfer function based on the geometric random channel model and the channel transfer function based on subspace reconstruction respectively; df is the calculus operator.

[0093] Step S602: Doppler power spectrum density is a short-time Fourier transform of the time autocorrelation function, and the calculation method is expressed as:

[0094] S qp (t, v) = ∫R qp (t, Δt)w(t-Δt)e -j2πvΔt dΔt

[0095] Among them, R qp (t, Δt) represents the time autocorrelation function of the channel transfer function at the zero frequency point, and w(t-Δt) represents the time window function;

[0096] Step S603: The time delay power spectrum density is the inverse Fourier transform of the frequency correlation function, which is also the square of the channel impulse response modulus. The calculation method is expressed as:

[0097] S qp (t, τ) = ∫R qp (t,Δf)e j2πτΔf dΔf=||h qp (t, τ)|| 2

[0098] Among them, R qp (t, Δf) represents the frequency correlation function of the channel transfer function, h qp (t, τ) represents the time-varying impulse response of the channel.

[0099] The present invention provides a non-stationary channel simulator design method based on subspace reconstruction, which has the following advantages: the present invention proposes a non-stationary geometric random channel model preprocessing method that projects a time-domain non-stationary fading sequence onto a subspace spanned by a Chirp signal, thereby performing channel reconstruction on an FPGA using a set of basis parameters, solving the problem that a channel simulator based on a geometric random channel model needs to calculate a large number of time-varying subpath parameters; further, the present invention proposes a non-stationary MIMO channel simulator based on subspace reconstruction and frequency domain processing, first performing subspace projection preprocessing on the channel transfer function calculated by the non-stationary geometric random channel model, and then reconstructing the channel transfer function on the FPGA using basis parameters obtained by subspace projection, thereby performing frequency-domain processing-based channel application on the channel simulator input signal using the reconstructed channel transfer function, solving the problem of high resource consumption of large-scale MIMO channel simulators under traditional TDL architectures. The channel simulator proposed by the present invention can be used to accurately and stably test large-scale MIMO communication systems in mobile scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 A schematic flow chart of the method for designing a non-stationary MIMO channel simulator based on subspace reconstruction provided in Example 1;

[0101] Figure 2 A schematic diagram of a non-stationary MIMO channel model based on geometric randomness provided in Example 1;

[0102] Figure 3 A schematic diagram of the hardware architecture of the non-stationary MIMO channel simulator based on subspace reconstruction provided in Example 1;

[0103] Figure 4 A schematic diagram of a frequency domain processing module of a non-stationary MIMO channel simulator hardware architecture based on subspace reconstruction provided in Example 1;

[0104] Figure 5 is a schematic diagram of the channel transfer function verification result provided in Example 1, wherein f c =2.6GHz, B=60MH z , N=23, M=20, K=30;

[0105] Figure 6 Schematic diagram of the verification result of the average error of the channel transfer function provided in Example 1, where f c =2.6GHz, B=60MHz, N=23, M=20, K=30;

[0106] Figure 7Schematic diagram of the simulation results of the Doppler power spectrum density based on the geometric random channel model and the simulation results based on the subspace reconstruction provided in Example 1, wherein f c =2.6GHz, v R =(10,0,0)m / s;

[0107] Figure 8 Schematic diagram of the simulation results of the delay power spectrum density based on the geometric random channel model and the simulation results based on the subspace reconstruction provided in Example 1, wherein f c =2.6GHz, v R =(10,0,0)m / s. DETAILED DESCRIPTION

[0108] In order to better understand the purpose, structure and function of the present invention, the following is a further detailed description of a non-stationary channel simulator design method based on subspace reconstruction of the present invention in conjunction with the accompanying drawings.

[0109] Example 1

[0110] See also Figure 1-Figure 8 This embodiment provides a design method for a non-stationary channel simulator based on subspace reconstruction, and the channel simulator design method includes the following steps:

[0111] Step 1: Configure the channel model environment and basic parameters such as the transmitting and receiving antennas on the host computer;

[0112] Specifically, in this embodiment, the carrier frequency is set to 2.6 GHz, and the channel bandwidth is set to 60 MHz. The initial position coordinates of the transmitting end are (0, 0, 35) m, the transmitting end remains stationary, and the receiving end moves at a speed of (10, 0, 0) m / s. The spherical coordinates of the center positions of the dual clusters of the transmitting and receiving ends are (d n ,φ A，n ,φ E，n ), where φ A，n and φ A，n According to the angle parameter configuration of CDL-A in 3GPP standard document 38901, d n Generated by exponential distribution. The transmitting antenna is a ULA with 128 antenna elements and a half wavelength spacing. The receiving antenna is a ULA with 4 antenna elements and a half wavelength spacing. The transmitting antenna azimuth is The elevation angle is 0, the receiving antenna azimuth is 0, and the elevation angle is Both the transmitting antenna and the receiving antenna are omnidirectional antennas.

[0113] Step 2: Generate scatterer position, subpath delay, and subpath power using a universal geometric random channel model to calculate the channel transfer function;

[0114] Specifically, in this embodiment, the model adopts a non-stationary MIMO channel model. The specific channel model schematic diagram is shown in Figure 2 In this embodiment, generating the channel transfer function specifically includes the following steps:

[0115] Step S201, calculate the initial position coordinates of the pth antenna unit at the transmitting end:

[0116]

[0117] in, is the initial position coordinate of the first antenna unit at the transmitting end, δ T is the transmitting antenna spacing, and are the azimuth and elevation angles of the transmitting antenna. Similarly, the coordinates of the qth antenna unit at the receiving end can be calculated.

[0118] Step S202: Generate scatterer position coordinates. The scatterer coordinates (x′, y′, z′) are distributed according to the ellipsoidal Gaussian distribution with the center of the dual cluster of the transmitting and receiving ends as the Cartesian coordinate origin. The joint probability density of the scatterer coordinates (x′, y′, z′) is:

[0119]

[0120] Among them, σ AS represents the cluster azimuthal spread, σ ES represents the cluster pitch angle spread, σ DS = represents the cluster delay spread. Thus, the Cartesian coordinates of the mth scatterer in the nth transmitting (receiving) cluster relative to the first antenna element of the transmitting (receiving) end are It can be obtained by the following coordinate transformation:

[0121]

[0122] Among them, (d n ,φ A，n ,φ E，n ) is the spherical coordinate of the center of the nth cluster at the transmitting (receiving) end.

[0123] Step S203, calculate the multipath propagation distance and delay between the antenna units at the transmitting and receiving ends. Based on the initial coordinates and moving speed of the antenna units at the transmitting and receiving ends, the initial coordinates of the scatterers in the dual clusters, and the moving speed of the clusters, the distance from the p-th transmitting antenna unit to the m-th scatterer in the n-th cluster at time t can be obtained. The calculation method is expressed as:

[0124]

[0125] in, and are the moving speed of the transmitting cluster and the moving speed of the transmitting antenna respectively. Similarly, the distance from the qth receiving antenna unit to the mth scatterer in the nth cluster at time t can be obtained. Therefore, the calculation method of the multipath propagation delay between the antenna units at the transmitting and receiving ends is:

[0126]

[0127] Where c is the speed of light, is the distance between the two clusters at the transmitting and receiving ends, τ link is the virtual link delay that follows exponential distribution.

[0128] Step S204: Generate the power of each sub-path according to the sub-path delay obtained in step S203. The calculation method is expressed as:

[0129]

[0130] Among them, f c is the carrier frequency, DS is the delay spread, r τ and γ are factors that depend on the time delay and frequency domain, respectively.

[0131] Step S205: Generate a channel transfer function according to the sub-path delay and sub-path power obtained in step S203 and step S204, and the calculation method is:

[0132]

[0133] Among them, f c is the carrier frequency, N qp (t) and M n (t) are the number of clusters and the number of scatterers in the cluster, respectively. According to the simulation time, channel bandwidth, antenna array size, and based on the birth and death rate of the dual clusters at the transmitting and receiving ends in space and time, the index of the surviving cluster is determined, and the position information of the newly generated clusters that change with the antenna unit interval, time interval, and frequency interval is calculated according to step S202, and then the sub-path delay and sub-path power of the corresponding cluster are calculated according to steps S203 and S204, and then the non-stationary channel transfer function that evolves with space and time and frequency is calculated.

[0134] Step 3: Project the channel transfer function generated in the previous step onto the subspace spanned by the Chirp signal to obtain a set of Chirp signal basis parameters;

[0135] In this embodiment, step 3 specifically includes:

[0136] Step S301: When the transceiver and the dual clusters move at a fixed speed and direction, the distance calculation in step S203 can be approximated as:

[0137]

[0138] in, is the initial distance from the first antenna unit at the transmitting end to the mth scatterer in the nth cluster at the transmitting end, and It is a parameter that depends on the antenna array, the moving speed and the angle of the scatterer.

[0139] Step S302: Substituting the approximate calculation relationship in step S301 into step S204, the subchannel transfer function of the MIMO channel can be approximately calculated as:

[0140]

[0141] Among them, η l ,ξ l , t l , a l , b l are all constant factors. Similarly, the Doppler frequency can be approximated and the calculation method is expressed as:

[0142]

[0143] Step S303: According to the approximate relationship in step S302, the Doppler frequency of the channel model can be approximated by a linear time-varying function, and the channel transfer function can be obtained by a linear combination of multiple Chirp signals. Therefore, a K-dimensional subspace is constructed using the Chirp signal as a basis function, K<<L(t). Specifically, in this embodiment, K=30. The kth basis function of the K-dimensional subspace is expressed as:

[0144]

[0145] Among them, α k represents the initial frequency, β k represents the Chirp rate. According to the approximate relationship in step S302, considering the case where the transmitting end is stationary and the receiving end is moving, α k and β k The value range of can be expressed as:

[0146]

[0147]

[0148] Thus α k and β kThe value of can be obtained by the orthogonal matching pursuit algorithm according to its value range or by the statistical distribution of the coefficient corresponding to the approximately linear time-varying Doppler frequency in the channel model.

[0149] Step S304: Determine the projection coefficient of the channel transfer function on the basis function in step S303. Consider any subchannel in the MIMO channel and express the channel fading at different frequencies as:

[0150] h=[H(t,f1),H(t,f2),…,H(t,f I )]

[0151] The basis function of the subspace in step S303 is expressed as:

[0152] Ψ=[b1(t),b2(t),…,b k (t)]

[0153] Assume that Ψ can be orthogonalized from matrix G to matrix Φ, that is, Φ = Ψ·G, where G can be expressed as:

[0154] in,

[0155] <·,·> represents the inner product operation, b k (t) is the kth column vector of Ψ, g j (t) is the jth column vector of Φ. The projection x of h on Φ is expressed as:

[0156]

[0157] x ji is the jth row and ith column element of X. Therefore, the projection A of the channel fading h on the Chirp basis function Ψ is expressed as:

[0158] A=G·X

[0159] Step 4: Send and store the basic parameters obtained in the previous step to the FPGA, and obtain the reconstructed channel transfer function by superimposing the Chirp signal on the FPGA;

[0160] In this embodiment, step 4 specifically includes:

[0161] The channel matrix preprocessing results based on subspace projection in step S303 and step S304 are Send and store to FPGA, so that this set of basic parameters in the time window T w The channel transfer function is reconstructed in , and the calculation method is expressed as:

[0162]

[0163] in, is the reconstructed channel transfer function.

[0164]

[0165] Among them, a k (f i ) is the kth row and ith column element of A. Therefore, low-complexity channel transfer function calculation can be completed on FPGA by superposition of time-division Chirp signals.

[0166] Step 5: Perform frequency domain processing on the MIMO channel simulator input signal according to the channel transfer function reconstructed in the previous step;

[0167] The hardware block diagram of the MIMO channel simulator architecture and frequency domain processing based on subspace reconstruction and frequency domain processing is shown in Figure 3 and Figure 4 First, the channel transfer function calculated by the non-stationary geometric random channel model is preprocessed by subspace projection, and then the channel transfer function is reconstructed by the basis parameters obtained by the subspace projection on the FPGA, so that the channel simulator input signal is subjected to channel application based on frequency domain processing by the reconstructed channel transfer function. In this embodiment, step 5 specifically includes:

[0168] Step S501: The channel transfer function calculated in step S4 Get the channel time sampling interval T ch The channel frequency response within H qp (f) According to the predefined channel bandwidth and frequency sampling interval in the channel model, the actual channel frequency response is a discrete H qp (m).

[0169] Step S502: From the input signal stream s p Take N from (n) s point signals p，i (n) put into cache, N s The calculation method is expressed as:

[0170] N s =N H -N α

[0171] Among them, N H is the number of channel frequency response points, means round down, τ m =max{τ|h qp (t, τ) ≥ 0}, h(t, τ) represents the channel impulse response, and is the channel transfer function H qp The inverse Fourier transform of (t,f) with respect to f.

[0172] Step S503: In the cache N s The input signal of point s p，i (n)N a Zeros and perform FFT, the number of points is N H The input signal spectrum S p，i (m). Then the spectrum of the qth output signal of the P×QMIMO channel simulator is expressed as:

[0173]

[0174] Step S504: q，i (m)N H IFFT of the points, and get the Q-channel output signal r q，i (n). Using the overlap-add method, the qth continuous output signal stream of the channel simulator is expressed as

[0175]

[0176] Step S505: According to step S4, the channel frequency response H is refreshed. qp (m), and repeat steps S501-S504. Channel update time T ch , that is, the channel time sampling interval, the calculation method is expressed as:

[0177] T ch =N s ·T s

[0178] Among them, T s is the input signal sampling period.

[0179] Step 6: Verify the consistency of the channel transfer function, Doppler power spectral density and delay power spectral density of the reconstructed channel and the original channel;

[0180] In this embodiment, step 6 specifically includes:

[0181] Step S601: The average error between the channel transfer function based on the subspace reconstruction and the channel transfer function based on the geometric random channel model is calculated as follows:

[0182]

[0183] Where B represents the channel bandwidth, H(t, f) and They represent the channel transfer function based on the geometric random channel model and the channel transfer function based on subspace reconstruction respectively.

[0184] Step S602: Doppler power spectrum density is a short-time Fourier transform of the time autocorrelation function, and the calculation method is expressed as:

[0185] S qp (t, v) = ∫R qp (t, Δt)w(t-Δt)e -j2πvΔt dΔt

[0186] Among them, R qp (t, Δt) represents the time autocorrelation function of the channel transfer function at the zero frequency point, and w(t-Δt) represents the time window function.

[0187] Step S603: The time delay power spectrum density is the inverse Fourier transform of the frequency correlation function, which is also the square of the channel impulse response modulus. The calculation method is expressed as:

[0188] S qp (t, τ) = ∫R qp (t,Δf)e j2πτΔf dΔf=||h qp (t, τ)|| 2

[0189] Among them, R qp (t, Δf) represents the frequency correlation function of the channel transfer function, h qp (t, τ) represents the time-varying impulse response of the channel.

[0190] In order to verify the correctness of the method provided in this embodiment, an experiment was conducted, specifically:

[0191] Verify the channel transfer function of the sub-channel between the transmitting and receiving antennas. The simulation results of the sub-channel transfer function between the fourth transmitting antenna and the first receiving antenna are as follows: Figure 5 (a) shows the simulation results. Figure 5 (b) shows that the amplitude of the channel transfer function varies with time and frequency, indicating the time selectivity and frequency selectivity of the wireless channel. Figure 6 is the average error between the channel transfer function based on subspace reconstruction and the channel transfer function based on the geometric random channel model, and its error is less than -90dB. It can be seen that the channel transfer function obtained from subspace reconstruction can well approximate the channel transfer function generated by the geometric random channel model. Due to the consistency between the subchannel transfer functions, the MIMO channel can be accurately simulated with low computational complexity by the subchannel transfer function based on subspace reconstruction between each antenna pair.

[0192] Verify the Doppler power spectrum density at different times. The results are as follows Figure 7As shown. The Doppler power spectrum density output by the channel simulator at t=0s, t=5s, and t=10s is compared with the results of the channel model simulation. The difference in Doppler power spectrum density at different times indicates that the time autocorrelation function is different in different time intervals, thus reflecting the time domain non-stationary characteristics of the channel. In addition, the channel simulation results are consistent with the channel model simulation results, indicating that the designed channel simulator can reproduce the time domain characteristics of the channel.

[0193] Verify the delay power spectrum density at different times. The channel model simulation results are as follows: Figure 8 (a) shows the channel simulation results. Figure 8 (b) As shown in the figure, the channel simulation results are in good agreement with the channel model simulation results. The delay spread of the delay axis in the figure shows the frequency selectivity of the multipath channel. The time-varying characteristics of the delay power spectrum density in the figure show the birth and death process of the effective subpath in the channel.

[0194] In summary, the present invention provides a method for designing a non-stationary MIMO channel simulator based on subspace reconstruction. Reconstructing the channel transfer function based on the subspace spanned by the Chirp signal can avoid the need for the non-stationary geometric random channel model to generate a large number of subpath sums to calculate the channel transfer function, thereby reducing the computational complexity of obtaining the channel transfer function. The designed non-stationary MIMO channel simulator directly uses the channel transfer function based on subspace reconstruction to apply channel fading to the input signal in the frequency domain, abandoning the traditional TDL-based architecture, and can achieve higher resource utilization and scalability for high-order MIMO channel simulation. In addition, the statistical characteristics of the simulation are of great reference value for the design of non-stationary channel simulators.

[0195] It is to be understood that the present invention is described by some embodiments, and it is known to those skilled in the art that various changes or equivalent substitutions may be made to these features and embodiments without departing from the spirit and scope of the present invention. In addition, under the teachings of the present invention, these features and embodiments may be modified to adapt to specific circumstances and materials without departing from the spirit and scope of the present invention. Therefore, the present invention is not limited by the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are within the scope of protection of the present invention.

Claims

1. A non-stationary channel simulator design method based on subspace reconstruction, characterized in that: It includes the following steps: Step S1: Based on the general geometric random channel model, configure the environment and basic parameters of the transceiver antennas in the upper computer for the channel model; Step S2: Generate a three-dimensional spatio-temporal-frequency non-stationary channel transfer function based on the general geometric random channel model; Step S3: Project the channel transfer function onto the subspace spanned by the chirp signals to obtain a set of chirp signal basis parameters; Step S4: Transmit and store the basis parameters obtained in Step S3 onto the field programmable gate array (FPGA), and obtain the reconstructed channel transfer function by superimposing chirp signals on the FPGA; Step S5: Perform frequency domain processing on the input signal of the MIMO channel simulator according to the reconstructed channel transfer function in Step S4; Step S6: Verify the consistency of the channel transfer function, Doppler power spectral density, and delay power spectral density between the reconstructed channel and the original channel.

2. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 1, characterized in that: The specific content of Step S1 includes: The environment and basic parameters of the transceiver antennas configured in the upper computer include simulation time, time sampling interval, carrier frequency, channel bandwidth, frequency sampling interval, initial position coordinates of the transceiver, moving speed of the transceiver, position coordinates of the double cluster center of the transceiver, moving speed of the double cluster of the transceiver, birth and death rates of the double cluster of the transceiver, number of transceiver antennas and antenna spacing, azimuth angle and elevation angle of the transceiver antennas.

3. The non-stationary channel simulator design method based on subspace reconstruction according to claim 1 is characterized in that: The specific content of Step S2 includes: Step S201: Calculate the initial position coordinates of the p-th antenna element at the transmitting end; in, is the initial position coordinate of the first antenna unit at the transmitting end, δ T is the transmitting antenna spacing, and are the azimuth and elevation angles of the transmitting antenna; calculate the coordinates of the qth antenna unit at the receiving end; in, is the initial position coordinate of the first antenna unit at the receiving end, δ R is the receiving antenna spacing, and is the azimuth and elevation angle of the receiving antenna; Step S202: Generate the position coordinates of the scatterers. The position coordinates of the scatterers (x', y', z') follow an ellipsoidal Gaussian distribution with the center of the double cluster of the transceiver as the origin of the Cartesian coordinates. The joint probability density of the position coordinates of the scatterers (x', y', z') is: Among them, σ AS represents the cluster azimuthal spread, σ ES represents the cluster pitch angle spread, σ DS represents the cluster delay spread; the Cartesian coordinates of the mth scatterer in the nth transmitting cluster relative to the first transmitting antenna element The coordinate transformation is obtained by the following formula: in, are the spherical coordinates of the centers of the n originating clusters; Cartesian coordinates of the mth scatterer in the nth receiving cluster relative to the first receiving antenna element The coordinate transformation is obtained by the following formula: in, are the spherical coordinates of the centers of the n clusters at the receiving end; Step S203: Calculate the multipath propagation distance and delay between each antenna element at the transceiver. Based on the initial coordinates and moving speed of the transceiver antenna elements, and the initial coordinates and cluster moving speed of the scatterers in the double cluster, obtain the distance from the p-th transmitting antenna element to the m-th scatterer in the n-th cluster at time t. The calculation method is expressed as: in, and are the moving speed of the transmitting cluster and the moving speed of the transmitting antenna respectively; the distance from the qth receiving antenna unit to the mth scatterer in the nth cluster at time t is obtained, and the calculation method is expressed as: in, and are the moving speed of the receiving cluster and the moving speed of the receiving antenna respectively; the sub-path delay caused by multipath propagation between the antenna units at the transmitting and receiving ends is calculated as follows: Where c is the speed of light, is the distance between the two clusters at the transmitting and receiving ends, τ link is the virtual link delay that obeys exponential distribution; Step S204: Generate the power of each sub-path according to the sub-path delay obtained in Step S203. The calculation method is expressed as: Among them, f c is the carrier frequency, DS is the delay spread, r τ and γ are factors that depend on the time delay and frequency domain, respectively; Step S205: Generate the channel transfer function according to the sub-path delay and sub-path power obtained in Steps S203 and S204. The calculation method is: Among them, f c is the carrier frequency, N qp (t) and M n (t) are the number of clusters and the number of scatterers in the cluster respectively; according to the simulation time, channel bandwidth, antenna array size, and based on the birth and death rate of the dual clusters at the transmitting and receiving ends in the space-time-frequency domain, the index of the surviving cluster is determined, and according to step S202, the position information of the newly generated cluster that changes with the antenna unit interval, time interval, and frequency interval is calculated, and then according to steps S203 and S204, the sub-path delay and sub-path power of the corresponding cluster are calculated, and then according to step S205, the non-stationary channel transfer function updated with the space-time-frequency evolution is calculated.

4. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 3 is characterized in that: The specific content of Step S3 includes: Step S301: When the transceiver and the double cluster move at a fixed speed magnitude and direction, the distance calculation in Step S203 is approximated as: in, is the fixed speed of the originating cluster, is the fixed speed of the transmitting antenna, is the initial distance from the first antenna unit at the transmitting end to the mth scatterer in the nth cluster at the transmitting end, and is a parameter that depends on the antenna array, the moving speed and the angle of the scatterer; in, The fixed speed of the receiving cluster is is the fixed speed of the receiving antenna, is the initial distance from the first antenna unit at the receiving end to the mth scatterer in the nth cluster at the receiving end, and is a parameter that depends on the antenna array, the moving speed and the angle of the scatterer; Step S302: According to the approximate calculation relationship in Step S301, substitute it into Step S204 to obtain the approximate calculation of the channel transfer function of the MIMO channel as: Where L(t) = max{N qp (t)·M n (t)},η qp,l ,ξ qp,l ,t qp,l ,a qp,l ,b qp,l are all constant factors; j is the imaginary unit; the approximate calculation method of the Doppler frequency of the channel model is expressed as: Step S303: According to the approximate relationship in Step S302, the Doppler frequency of the channel model is approximated by a linear time-varying function, and the channel transfer function is obtained by the linear combination of multiple chirp signals; thus, a K-dimensional subspace is constructed with chirp signals as the basis functions, k << L(t); the k-th basis function of the K-dimensional subspace is expressed as: Among them, α k represents the initial frequency, β k represents the Chirp rate; according to the approximate relationship in step S302, considering the case where the transmitting end is stationary and the receiving end is moving, α k and β k The value range of is expressed as: Thus α k and β k The value of is obtained from the statistical distribution of the coefficient corresponding to the approximately linear time-varying Doppler frequency in the channel model; Step S304, determine the projection coefficient of the channel transfer function on the basis function in step S303; consider any subchannel in the MIMO channel, and express the channel fading at different frequencies as: h=[H(t,f1),H(t,f2),…,H(t,f I )] Among them, H(t,f i ), i = 1, 2, ..., I is the H (t, f) of the MIMO subchannel at frequency f i The sample value on ; The basis function of the subspace in step S303 is expressed as: Ψ=[b1(t),b2(t),…,b K (t)] Among them, b k (t), k = 1, 2, ..., K is the basis function in step S303; Assume that Ψ is orthogonalized from matrix G to matrix Φ, that is, Φ = Ψ·G, where G is expressed as: in, <,> represents the inner product operation, b k (t) is the kth column vector of Ψ, g j (t) is the jth column vector of Φ; the projection X of h on Φ is expressed as: x ji is the jth row and ith column element of X; thus, the projection A of the channel fading h on the Chirp basis function Ψ is expressed as: A=G·X.

5. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 4 is characterized in that: The step S4 specifically includes: The channel matrix preprocessing results based on subspace projection in step S303 and step S304 are Send and store to FPGA, so that this set of basic parameters in the time window T w The channel transfer function is reconstructed in , and the calculation method is expressed as: in, is the reconstructed channel transfer function, that is Among them, a k (f i ) is the kth row and ith column element of A, t0 is the initial time of channel simulation, T w is the duration of the channel simulation; therefore, a low-complexity channel transfer function calculation is completed on the FPGA by superposition of time-division Chirp signals.

6. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 5, characterized in that: The step S5 specifically includes: Step S501: The channel transfer function calculated in step S4 Get the channel time sampling interval T ch The channel frequency response within H qp (f i );According to the predefined channel bandwidth and frequency sampling interval in the channel model, the actual channel frequency response is a discrete H qp (m); Step S502: From the input signal stream s p Take N from (n) s point signals p,i (n) put into cache, N s The calculation method is expressed as: N s =N H -N a Among them, N H is the number of channel frequency response points, means round down, τ m =max{τ|h qp (t,τ)≥0}, h(t,τ) represents the channel impulse response, and is the channel transfer function H qp The inverse Fourier transform of (t,f) with respect to f; Step S503: In the cache N s The input signal of point s p,i (n)N a Zeros and perform fast Fourier transform FFT to get N points H The input signal spectrum S p,i (m); then the q-th output signal spectrum of the P×QMIMO channel simulator is expressed as: Step S504: q,i (m)N H The inverse fast Fourier transform IFFT of the point is used to obtain the Q-channel output signal r q,i (n); Using the overlap-add method, the qth continuous output signal stream of the channel simulator is expressed as Step S505: According to step S4, the channel frequency response H is refreshed. qp (m), and repeat steps S501-S504; channel update time T ch , that is, the channel time sampling interval, the calculation method is expressed as: T ch =N s ·T s Among them, T s is the input signal sampling period.

7. The method for designing a non-stationary channel simulator based on subspace reconstruction according to claim 6, characterized in that: The step S6 specifically includes: Step S601: The average error between the channel transfer function based on the subspace reconstruction and the channel transfer function based on the geometric random channel model is calculated as follows: Where B represents the channel bandwidth, H(t,f) and They represent the channel transfer function based on the geometric random channel model and the channel transfer function based on subspace reconstruction respectively; Step S602: Doppler power spectrum density is a short-time Fourier transform of the time autocorrelation function, and the calculation method is expressed as: s qp (t,v)=∫R qp (t,Δt)w(t-Δt)e -j2πvΔtd Δt Among them, R qp (t, Δt) represents the time autocorrelation function of the channel transfer function at the zero frequency point, and w(t-Δt) represents the time window function; Step S603: The time delay power spectrum density is the inverse Fourier transform of the frequency correlation function, which is also the square of the channel impulse response modulus. The calculation method is expressed as: S qp (t,τ)=∫R qp (t,Δf)e j2πτΔf dΔf=||h qp (t,τ)|| 2 Among them, R qp (t,Δf) represents the frequency correlation function of the channel transfer function, h qp (t,τ) represents the time-varying impulse response of the channel.

Images