Description of MassiveMIMOAsymptoticPerformance

0001 % Massive MIMO FBMC-OQAM performance under highly frequency
0002 % selective channel
0003 %
0004 % Example showing massive MIMO FBMC-OQAM performance under highly frequency
0005 % selective channel for a ZF equalizer. As N and K grow large, the MSE gets
0006 % flat across frequency and converge to the asymptotic one [1].
0007 %
0008 % Reference [1] F. Rottenberg, X. Mestre, F. Horlin and J. Louveaux,
0009 % "Performance Analysis of Linear Receivers for Uplink Massive MIMO
0010 % FBMC-OQAM Systems," in IEEE Transactions on Signal Processing, vol. 66,
0011 % no. 3, pp. 830-842, 1 Feb.1, 2018. doi: 10.1109/TSP.2017.2778682
0012 
0013 
0014 
0015 % This file is part of WaveComBox: www.wavecombox.com and is distributed under
0016 % the terms of the MIT license. See accompanying LICENSE file.
0017 % Original author: François Rottenberg, June 24, 2019.
0018 % Contributors:
0019 % Change log:
0020 
0021 % Definition of general Parameters
0022 Para=InitializeChainParameters( 'FBMC-OQAM' ); % initialize parameters of the FBMC-OQAM chain
0023 Para.Ns=100;
0024 K=15; % Number of users
0025 N=30; % Number of BS antennas
0026 Para.Es_N0_dB=20;
0027 Para.N_R=N;
0028 Para.N_T=K;
0029 Para.S=min(Para.N_T,Para.N_R);
0030 
0031 %% Transmitter
0032 d = GenerateData (Para); % Generate data symbols
0033 s = Modulator(d, Para ); % FBMC-OQAM modulation
0034 
0035 %% Channel
0036 
0037 % BS antenna correlation
0038 r=0.7;
0039 row=r.^[0:1:N-1];
0040 C_BS = toeplitz(row);
0041 [V,D] = eig(C_BS);
0042 C_BS_1_2=V*sqrt(D)*V';
0043 
0044 % Large-scale fading coefficients and PDP (same PDP per user, choose one of the following by uncommenting)
0045 beta=ones(K,1)-0.5*rand(K,1);
0046 beta=beta./sum(beta)*K;
0047 D_beta=diag(beta);
0048 
0049 
0050 % % % % ***************************** Vehicular A  (speed 60 km/h)
0051 % delays=1e-9*[0 300 700 1100 1700 2500];
0052 % powers=[0 -1 -9 -10 -15 -20];
0053 %
0054 % % ***************************** Vehicular B  (speed 60 km/h)
0055 delays=1e-9*[0 300 8900 12900 17100 20000];
0056 powers=[-2.5 0 -12.8 -10 -25.2 -16];
0057 %
0058 % % %****************************** model TU 50 model A
0059 % delays=1e-6*[0 0.217 0.512 0.514 0.517 0.674 0.882 1.230 1.287 1.311 1.349 1.533 1.535 1.622 1.818 1.836 1.884 1.943 2.048 2.140] ;
0060 % powers=[ -5.7 -7.6 -10.1 -10.2 -10.2 -11.5 -13.4 -16.3 -16.9 -17.1 -17.4 -19.0 -19.0 -19.8 -21.5 -21.6 -22.1 -22.6 -23.5 -24.3];
0061 
0062 % % 3GPP-HT hilly terrain
0063 % powers = [-3.6 -8.9 -10.2 -11.5 -11.8 -12.7 -13 -16.2 -17.3 ...
0064 %     -17.7 -17.6 -22.7 -24.1 -25.8 -25.8 -26.2 -29 -29.9 -30 -30.7];
0065 % delays = [0 356 441 528 546 609 625 842 916 941 15000 16172 ...
0066 %     16492 16876 16882 16978 17615 17827 17849 18016]*1e-9;
0067 
0068 Ts=Para.T/Para.nSubcarriers;
0069 delays=round(delays/Ts)+1;
0070 nbtaps=length(powers);
0071 L=max(delays);
0072 p_b=zeros(1,L);
0073 variances=10.^(powers/10);
0074 variances=variances/sum(variances);
0075 for i=1:nbtaps
0076     p_b(delays(i))=p_b(delays(i))+ variances(i);
0077 end
0078 
0079 L=length(p_b);
0080 D_beta_b=zeros(K,K,L);
0081 D_1_2=zeros(K*L,K);
0082 for k=1:K
0083     for b=1:L
0084         D_beta_b(k,k,b)=beta(k)*(p_b(b));
0085     end
0086 end
0087 
0088 for b=1:L
0089     D_1_2(1+(b-1)*K:K+(b-1)*K,:)=sqrt(D_beta_b(:,:,b));
0090 end
0091 
0092 % Channel statitics
0093 Psi=diag(0:L-1);
0094 Phi=kron(-1j*Psi,eye(K,K));
0095 D=D_1_2'*D_1_2;
0096 D_av=D\D_1_2'*1j*Phi*D_1_2;
0097 D_rms=sqrt(-inv(D)*(D_1_2'*Phi^2*D_1_2)-D_av^2);
0098 
0099 % Small-scale fading and channel generation
0100 C=zeros(N,K,L);
0101 for b=1:L
0102     C(:,:,b)=C_BS_1_2*1/sqrt(2)*(randn(N,K)+1j*randn(N,K))*sqrt(D_beta_b(:,:,b));
0103 end
0104 
0105 % C=GenerateRayleighChannelReal('ITU_VehA', Para);
0106 r = Channel_Multipath( s, C ); % Multipath channel
0107 r=Channel_AWGN( r, Para ); % AWGN channel, SNR fixed in Para
0108 
0109 %% Receiver
0110 z = Demodulator( r, Para ); % FBMC demodulation
0111 
0112 % Single-tap equalizer
0113 ParaEqualizer.criterion='ZF';
0114 [ x, B ] = Equalizer( z, C, Para, 'SingleTap', ParaEqualizer );
0115 d_hat=real(x);
0116 
0117 %% Performance
0118 % Simulated MSE
0119 MSE = MSEComputes( d, d_hat, Para );
0120 
0121 % Theoretical approximation of the MSE
0122 MSE_theo_approx = FBMC_OQAM_MSESingleTapTheorApprox( B, C, Para );
0123 
0124 % Asymptotic equivalent of the MSE
0125 [eta_tilde, zeta_tilde]=ComputeEtaTilde(Para);
0126 
0127 % delta and delta_tilde computation ZF
0128 delta_tilde_prev=N/K-1;
0129 precision=1;
0130 while abs(precision)>1e-3
0131     delta_tilde=1/K*trace(C_BS/( delta_tilde_prev^(-1)*C_BS  +eye(N,N)));
0132     precision=delta_tilde-delta_tilde_prev;
0133     delta_tilde_prev=delta_tilde;
0134 end
0135 delta=delta_tilde^(-1);
0136 alpha=1/K*trace(inv(D))*delta_tilde^(-2);
0137 gamma=delta_tilde^(-2);
0138 alpha_tilde=1/K*trace(C_BS*(inv( delta*C_BS  +eye(N,N)))^2 );
0139 gamma_tilde=1/K*trace(C_BS^2*(inv( delta*C_BS  +eye(N,N)))^2 );
0140 
0141 Para.Es_N0=10.^(Para.Es_N0_dB/10);
0142 sigma_2=Para.Es./Para.Es_N0;
0143 MSE_bar_ZF=sigma_2/2*alpha*alpha_tilde/(1-gamma*gamma_tilde)...
0144     +eta_tilde*trace(D_av^2)/K+...
0145     +eta_tilde*alpha*gamma_tilde/(1-gamma*gamma_tilde)*trace(D*D_rms^2)/K;
0146 
0147 
0148 %% plot results
0149 
0150 figure
0151 plot(Para.ActiveSubcarriers, 10*log10(MSE),'xb','markersize',8,'linewidth',1.5)
0152 hold on
0153 plot(Para.ActiveSubcarriers,10*log10(MSE_theo_approx(Para.ActiveSubcarriers)),'-b','markersize',8,'linewidth',1.5)
0154 plot(Para.ActiveSubcarriers,10*log10(MSE_bar_ZF)*ones(length(Para.ActiveSubcarriers),1),'--b','markersize',8,'linewidth',1.5)
0155 xlabel('Subcarrier index')
0156 ylabel('MSE [dB]')
0157 xlim([min(Para.ActiveSubcarriers) max(Para.ActiveSubcarriers)])
0158 legend('ZF - simu', 'ZF - approx','ZF - asymptotic')
0159 ylim([-20 0])
MassiveMIMOAsymptoticPerformance

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
