FBMC_OQAM_PreEqualizerSingleTap

 Pre-equalizes data symbols d and returns corresponding pre-equalizing and
 required equalizing matrices.

 function [d_tilde, A, B ]=FBMC_OQAM_PreEqualizerSingleTap( d, C,
 criterion, Para )

 The function works for SISO and MIMO systems.

 Input arguments:

   d: data symbols to pre-equalize. Size: matrix [Para.nSubcarriers, 2*Para.Ns] if
   Para.S == 1, multidimensional array [para.S, Para.nSubcarriers, 2*Para.Ns] if
   Para.S > 1.

   C: channel impulse response. Multidimensional array [Para.N_R, Para.N_T, ~].

   criterion: determines the design criterion to use.
   Several options are possible: 

       'ZF': zero forcing. Error if Para.N_R > Para.N_T.

       'MMSE': minimum mean squared error. Error if Para.N_R > Para.N_T.

       'MF': matched filter. Error if Para.N_R > Para.N_T.

       'SVD': singular value decomposition. Precode by S right singular
       vector of the channel frequency response corresponding to largest
       singular values.

       'ZF_opt_freq_selec': as opposed to conventional 'ZF', does not
       assume channel frequency flatness at the subcarrier level and takes
       into account derivatives of channel frequency response. Follows the
       design of [1].

       'MMSE_opt_freq_selec': as opposed to conventional 'MMSE', does not
       assume channel frequency flatness at the subcarrier level and takes
       into account derivatives of channel frequency response. Follows the
       design of [1].

   Para: structure containing the modulation parameters.

 Outputs arguments:

   d_tilde: pre-equalized data symbols. Size: matrix [Para.nSubcarriers, 2*Para.Ns] if
   Para.N_T == 1, multidimensional array [para.N_T, Para.nSubcarriers, 2*Para.Ns] if
   Para.N_T > 1.

   A: used pre-equalizing matrices at each subcarrier. Size: vector [Para.nSubcarriers, 1]
   if Para.N_T == Para.N_R == 1 and multidimensional array
   [Para.T,Para.S,Para.nSubcarriers] otherwise

   B: required equalizing matrices at each subcarrier. Size: vector [Para.nSubcarriers, 1]
   if Para.N_T == Para.N_R == 1 and multidimensional array [Para.S,Para.N_R,Para.nSubcarriers] otherwise


   References: [1]    F. Rottenberg, X. Mestre, F. Horlin and J. Louveaux,
   "Single-Tap Precoders and Decoders for Multi-User MIMO FBMC-OQAM under
   Strong Channel Frequency Selectivity," IEEE Transactions on Signal
   Processing, 2017, vol. 65, no. 3, pp. 587-600.

0001 function [d_tilde, A, B ]=FBMC_OQAM_PreEqualizerSingleTap( d, C, criterion, Para )
0002 % Pre-equalizes data symbols d and returns corresponding pre-equalizing and
0003 % required equalizing matrices.
0004 %
0005 % function [d_tilde, A, B ]=FBMC_OQAM_PreEqualizerSingleTap( d, C,
0006 % criterion, Para )
0007 %
0008 % The function works for SISO and MIMO systems.
0009 %
0010 % Input arguments:
0011 %
0012 %   d: data symbols to pre-equalize. Size: matrix [Para.nSubcarriers, 2*Para.Ns] if
0013 %   Para.S == 1, multidimensional array [para.S, Para.nSubcarriers, 2*Para.Ns] if
0014 %   Para.S > 1.
0015 %
0016 %   C: channel impulse response. Multidimensional array [Para.N_R, Para.N_T, ~].
0017 %
0018 %   criterion: determines the design criterion to use.
0019 %   Several options are possible:
0020 %
0021 %       'ZF': zero forcing. Error if Para.N_R > Para.N_T.
0022 %
0023 %       'MMSE': minimum mean squared error. Error if Para.N_R > Para.N_T.
0024 %
0025 %       'MF': matched filter. Error if Para.N_R > Para.N_T.
0026 %
0027 %       'SVD': singular value decomposition. Precode by S right singular
0028 %       vector of the channel frequency response corresponding to largest
0029 %       singular values.
0030 %
0031 %       'ZF_opt_freq_selec': as opposed to conventional 'ZF', does not
0032 %       assume channel frequency flatness at the subcarrier level and takes
0033 %       into account derivatives of channel frequency response. Follows the
0034 %       design of [1].
0035 %
0036 %       'MMSE_opt_freq_selec': as opposed to conventional 'MMSE', does not
0037 %       assume channel frequency flatness at the subcarrier level and takes
0038 %       into account derivatives of channel frequency response. Follows the
0039 %       design of [1].
0040 %
0041 %   Para: structure containing the modulation parameters.
0042 %
0043 % Outputs arguments:
0044 %
0045 %   d_tilde: pre-equalized data symbols. Size: matrix [Para.nSubcarriers, 2*Para.Ns] if
0046 %   Para.N_T == 1, multidimensional array [para.N_T, Para.nSubcarriers, 2*Para.Ns] if
0047 %   Para.N_T > 1.
0048 %
0049 %   A: used pre-equalizing matrices at each subcarrier. Size: vector [Para.nSubcarriers, 1]
0050 %   if Para.N_T == Para.N_R == 1 and multidimensional array
0051 %   [Para.T,Para.S,Para.nSubcarriers] otherwise
0052 %
0053 %   B: required equalizing matrices at each subcarrier. Size: vector [Para.nSubcarriers, 1]
0054 %   if Para.N_T == Para.N_R == 1 and multidimensional array [Para.S,Para.N_R,Para.nSubcarriers] otherwise
0055 %
0056 %
0057 %   References: [1]    F. Rottenberg, X. Mestre, F. Horlin and J. Louveaux,
0058 %   "Single-Tap Precoders and Decoders for Multi-User MIMO FBMC-OQAM under
0059 %   Strong Channel Frequency Selectivity," IEEE Transactions on Signal
0060 %   Processing, 2017, vol. 65, no. 3, pp. 587-600.
0061 
0062 % This file is part of WaveComBox: www.wavecombox.com and is distributed under the terms of the MIT license. See accompanying LICENSE file.
0063 % Original author: François Rottenberg, May 3, 2018.
0064 % Contributors:
0065 % Change log:
0066 
0067 if Para.N_R==1 && Para.N_T==1
0068     H=fft(squeeze(C),Para.nSubcarriers);
0069     switch criterion
0070         case 'ZF'
0071             A_tilde=1./H;
0072             xi=abs(A_tilde);
0073             A=A_tilde./xi;
0074             B=xi;
0075         case 'MMSE'
0076             % Noise power
0077             Es_N0=10.^(Para.Es_N0_dB/10);
0078             A_tilde=conj(H)./(abs(H).^2+1/Es_N0);
0079             xi=abs(A_tilde);
0080             A=A_tilde./xi;
0081             B=xi;
0082         case 'MF'
0083             A_tilde=1./H;
0084             xi=abs(A_tilde);
0085             A=A_tilde./xi;
0086             B=xi;
0087         case 'SVD'
0088             A_tilde=1./H;
0089             xi=abs(A_tilde);
0090             A=A_tilde./xi;
0091             B=xi;
0092         case 'ZF_opt_freq_selec'
0093             A_tilde=1./H;
0094             xi=abs(A_tilde);
0095             A=A_tilde./xi;
0096             B=xi;
0097         case 'MMSE_opt_freq_selec'
0098             L=length(C(1,1,:));
0099             H_1=-sqrt(-1)*fft((0:L-1).*squeeze(C).',Para.nSubcarriers).';
0100             H_2=-fft((0:L-1).^2.*squeeze(C).',Para.nSubcarriers).';
0101             Es_N0=10.^(Para.Es_N0_dB/10);
0102             eta_tilde=ComputeEtaTilde(Para);
0103             
0104             X_inv=1./(abs(H).^2+eta_tilde*(conj(H).*H_2+conj(H_2).*H)+2*eta_tilde*abs(H_1).^2+1/Es_N0);
0105             psi=-1./(real(abs(H).^2.*X_inv)).*imag(H.*X_inv.*(conj(H)+eta_tilde.*conj(H_2)));
0106             A_tilde=X_inv.*(1j*conj(H).*psi+conj(H)+eta_tilde.*conj(H_2));
0107             xi=abs(A_tilde);
0108             A=A_tilde./xi;
0109             B=xi;            
0110         otherwise
0111             error('Equalizer type not existing')
0112     end
0113     d_tilde=zeros(Para.nSubcarriers, 2*Para.Ns);
0114     for m=Para.ActiveSubcarriers
0115         d_tilde(m,:)=A(m).*d(m,:);
0116     end
0117 else    
0118     H=zeros(Para.N_R,Para.N_T,Para.nSubcarriers);
0119     for index_N_R=1:Para.N_R
0120         for index_N_T=1:Para.N_T
0121             H(index_N_R,index_N_T,:)=fft(squeeze(C(index_N_R,index_N_T,:)),Para.nSubcarriers).';
0122         end
0123     end
0124     B=zeros(Para.S, Para.N_R, Para.nSubcarriers);
0125     A=zeros(Para.N_T, Para.S, Para.nSubcarriers);
0126     switch criterion
0127         case 'ZF'
0128             if Para.N_R>Para.N_T
0129                 error('Error: N_R>N_T')
0130             end
0131             for m=Para.ActiveSubcarriers
0132                 H_m=H(:,:,m);
0133                 Gram_inv=eye(Para.N_R,Para.N_R)/(H_m*H_m');
0134                 xi=sqrt(trace(Gram_inv)/Para.S);
0135                 A(:,:,m)=H_m'*Gram_inv.*(1/xi);
0136                 B(:,:,m)=xi*eye(Para.N_R,Para.N_R);
0137             end
0138             
0139         case 'MMSE'
0140             if Para.N_R>Para.N_T
0141                 error('Error: N_R>N_T')
0142             end
0143             for m=Para.ActiveSubcarriers
0144                 % Noise power
0145                 Es_N0=10.^(Para.Es_N0_dB/10);
0146                 H_m=H(:,:,m);
0147                 A_tilde=(H_m'*H_m+1/Es_N0*Para.N_R*eye(Para.N_T,Para.N_T))\H_m';
0148                 xi=sqrt(trace(A_tilde*A_tilde')/Para.S);
0149                 A(:,:,m)=A_tilde.*(1/xi);
0150                 B(:,:,m)=xi*eye(Para.N_R,Para.N_R);
0151             end            
0152         case 'MF'
0153             if Para.N_R>Para.N_T
0154                 error('Error: N_R>N_T')
0155             end
0156             for m=Para.ActiveSubcarriers
0157                 H_m=H(:,:,m);
0158                 A_tilde=H_m'/(diag(diag(H_m*H_m')));                
0159                 xi=sqrt(trace(A_tilde*A_tilde')/Para.S);
0160                 A(:,:,m)=A_tilde.*(1/xi);
0161                 B(:,:,m)=xi*eye(Para.N_R,Para.N_R);
0162             end
0163         case 'SVD'
0164             for m=Para.ActiveSubcarriers
0165                 if Para.S>min(Para.N_T, Para.N_R)
0166                     error('Error: S>min(N_T,N_R)')
0167                 end
0168                 H_m=H(:,:,m);
0169                 [U,S,V] = svd(H_m);
0170                 A_tilde=V(:,1:Para.S);          
0171                 xi=sqrt(trace(A_tilde*A_tilde')/Para.S);
0172                 A(:,:,m)=A_tilde.*(1/xi);
0173                 B(:,:,m)=((H_m*A(:,:,m))'*(H_m*A(:,:,m)) )\((H_m*A(:,:,m))');
0174             end
0175         case 'ZF_opt_freq_selec'
0176             if Para.N_R>Para.N_T
0177                 error('Error: N_R>N_T')
0178             end
0179             Es_N0=10.^(Para.Es_N0_dB/10);
0180             eta_tilde=ComputeEtaTilde(Para);
0181             H_1=zeros(Para.N_R,Para.N_T,Para.nSubcarriers);
0182             L=length(C(1,1,:));
0183             for index_N_R=1:Para.N_R
0184                 for index_N_T=1:Para.N_T
0185                     H_1(index_N_R,index_N_T,:)=-sqrt(-1)*fft((0:L-1).*squeeze(C(index_N_R,index_N_T,:)).',Para.nSubcarriers).';
0186                 end
0187             end
0188             for m=Para.ActiveSubcarriers
0189                 H_m=H(:,:,m);
0190                 H1_m=H_1(:,:,m);
0191                 Gram_inv=eye(Para.N_R,Para.N_R)/(H_m*H_m');
0192                 P_plus=eye(Para.N_T,Para.N_T)-H_m'*Gram_inv*H_m;
0193                 H_plus=H_m'*Gram_inv;
0194                 P_plus_A2=-P_plus*H1_m'/(H1_m*P_plus*H1_m'+1/Es_N0*Para.N_R/eta_tilde*eye(Para.N_R,Para.N_R))*H1_m*H_plus; % Matrix lemma inversion verification
0195                 A_tilde=( H_plus+P_plus_A2 );
0196                 xi=sqrt(trace(A_tilde*A_tilde')/Para.S);
0197                 A(:,:,m)=A_tilde.*(1/xi);
0198                 B(:,:,m)=xi*eye(Para.N_R,Para.N_R);
0199             end
0200         case 'MMSE_opt_freq_selec'
0201             if Para.N_R>Para.N_T
0202                 error('Error: N_R>N_T')
0203             end
0204             Es_N0=10.^(Para.Es_N0_dB/10);
0205             eta_tilde=ComputeEtaTilde(Para);
0206             H_1=zeros(Para.N_R,Para.N_T,Para.nSubcarriers);
0207             H_2=zeros(Para.N_R,Para.N_T,Para.nSubcarriers);
0208             L=length(C(1,1,:));
0209             for index_N_R=1:Para.N_R
0210                 for index_N_T=1:Para.N_T
0211                     H_1(index_N_R,index_N_T,:)=-sqrt(-1)*fft((0:L-1).*squeeze(C(index_N_R,index_N_T,:)).',Para.nSubcarriers).';
0212                     H_2(index_N_R,index_N_T,:)=-fft((0:L-1).^2.*squeeze(C(index_N_R,index_N_T,:)).',Para.nSubcarriers).';
0213                 end
0214             end
0215             for m=Para.ActiveSubcarriers
0216                 H_m=H(:,:,m);
0217                 H1=H_1(:,:,m);
0218                 H2=H_2(:,:,m);                
0219                 X_inv=eye(Para.N_T,Para.N_T)/(H_m'*H_m+eta_tilde*(H_m'*H2+H2'*H_m)+2*eta_tilde*(H1'*H1)+1/Es_N0*Para.N_R*eye(Para.N_T,Para.N_T));
0220                 psi=-inv(real(H_m*X_inv*H_m'))*imag(H_m*X_inv*(H_m'+eta_tilde*H2'));
0221                 A_tilde=X_inv*(1j*H_m'*psi+H_m'+eta_tilde*H2');
0222                 xi=sqrt(trace(A_tilde*A_tilde')/Para.S);
0223                 A(:,:,m)=A_tilde.*(1/xi);
0224                 B(:,:,m)=xi*eye(Para.N_R,Para.N_R);
0225             end
0226             
0227         otherwise
0228             error('Equalizer type not existing')
0229     end
0230     d_tilde=zeros(Para.N_T,Para.nSubcarriers, 2*Para.Ns);
0231     if Para.S==1
0232         for m=Para.ActiveSubcarriers
0233             for l=1:2*Para.Ns
0234                 d_tilde(:,m,l)=A(:,:,m)*d(m,l);
0235             end
0236         end
0237     else
0238         for m=Para.ActiveSubcarriers
0239             for l=1:2*Para.Ns
0240                 d_tilde(:,m,l)=A(:,:,m)*d(:,m,l);
0241             end
0242         end
0243     end
0244     if Para.N_T==1
0245         d_tilde=squeeze(d_tilde);
0246     end
0247     
0248     
0249 end
0250 
0251 
0252 
0253 
0254 
0255 
0256 
0257 
0258 
0259 
0260 
0261 
0262 
0263 
0264 
0265 end
0266