FBMC_OQAM_EqualizerParallelMultistage

 Parallel multi-stage equalization of received FBMC-OQAM signal for
 time-invariant frequency selective channels.

 function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para )
 function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para, B )
 function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para, B, B1 )
 function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para, B, B1, B2 )

 The function works for MIMO and SISO systems.

 Input arguments:

   r: received signal. Size: vector [para.N_R, nFrameSamples].

   R: number of additional parallel analysis filterbanks (total = R+1),
   cannot be bigger than 2.

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

   criterion: determines the 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. 

       'Specific': the specific equalizing matrices given in the function
       argument B, B1, B2 should be provided and will be used for
       equalization. Argument C is not regarded.

   Para: structure containing the modulation parameters.

   B, B1, B2: the equalizing matrices and their derivatives at each
   subcarrier should be provided if 'specific' equalizer criterion is chosen.
   If not chosen, the matrices are disregarded. Size: vector [Para.nSubcarriers ,1]
   if Para.N_T == 1 and multi- dimensional array [Para.S, Para.N_R,
   Para.nSubcarriers] otherwise.

 Ouput arguments:

   x: equalized symbols. Size: matrix [Para.nSubcarriers, 2*Para.Ns] if Para.S == 1,
   multidimensional array [para.S, Para.nSubcarriers, 2*Para.Ns] if Para.S > 1.   


 References: [1]    X. Mestre and D. Gregoratti, "Parallelized Structures
   for MIMO FBMC under Strong Channel Frequency Selectivity," IEEE Trans.
   Signal Process., vol. PP, no. 99, pp. 1–1, 2015.

0001 function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para, B, B1, B2 )
0002 % Parallel multi-stage equalization of received FBMC-OQAM signal for
0003 % time-invariant frequency selective channels.
0004 %
0005 % function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para )
0006 % function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para, B )
0007 % function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para, B, B1 )
0008 % function [ x ] = FBMC_OQAM_EqualizerParallelMultistage( r, C, R, criterion, Para, B, B1, B2 )
0009 %
0010 % The function works for MIMO and SISO systems.
0011 %
0012 % Input arguments:
0013 %
0014 %   r: received signal. Size: vector [para.N_R, nFrameSamples].
0015 %
0016 %   R: number of additional parallel analysis filterbanks (total = R+1),
0017 %   cannot be bigger than 2.
0018 %
0019 %   C: tme-invariant channel impulse response. Multidimensional array
0020 %   [Para.N_R, Para.N_T, ~].
0021 %
0022 %   criterion: determines the criterion to use. Several
0023 % options are possible:
0024 %
0025 %       'ZF': zero forcing. Error if Para.N_R < Para.N_T.
0026 %
0027 %       'MMSE': minimum mean squared error. Error if Para.N_R < Para.N_T.
0028 %
0029 %       'Specific': the specific equalizing matrices given in the function
0030 %       argument B, B1, B2 should be provided and will be used for
0031 %       equalization. Argument C is not regarded.
0032 %
0033 %   Para: structure containing the modulation parameters.
0034 %
0035 %   B, B1, B2: the equalizing matrices and their derivatives at each
0036 %   subcarrier should be provided if 'specific' equalizer criterion is chosen.
0037 %   If not chosen, the matrices are disregarded. Size: vector [Para.nSubcarriers ,1]
0038 %   if Para.N_T == 1 and multi- dimensional array [Para.S, Para.N_R,
0039 %   Para.nSubcarriers] otherwise.
0040 %
0041 % Ouput arguments:
0042 %
0043 %   x: equalized symbols. Size: matrix [Para.nSubcarriers, 2*Para.Ns] if Para.S == 1,
0044 %   multidimensional array [para.S, Para.nSubcarriers, 2*Para.Ns] if Para.S > 1.
0045 %
0046 %
0047 % References: [1]    X. Mestre and D. Gregoratti, "Parallelized Structures
0048 %   for MIMO FBMC under Strong Channel Frequency Selectivity," IEEE Trans.
0049 %   Signal Process., vol. PP, no. 99, pp. 1–1, 2015.
0050 %
0051 
0052 
0053 % This file is part of WaveComBox: www.wavecombox.com and is distributed under the terms of the MIT license. See accompanying LICENSE file.
0054 % Original author: François Rottenberg, May 4, 2018.
0055 % Contributors:
0056 % Change log:
0057 
0058 
0059 if R>2
0060     error('Case R>2 not implemented yet')
0061 end
0062 M=Para.nSubcarriers/2;
0063 
0064 theta=ones(2*M,2*(Para.Ns));
0065 theta(2:2:end,1:2:end)=1j;
0066 theta(1:2:end,2:2:end)=1j;
0067 [pulse, pulse10, pulse20, ~] = GeneratePrototypePhydyas(M, Para.kappa);
0068 pulse=pulse.*sqrt(1/M);
0069 pulse10=pulse10.*sqrt(1/M);
0070 pulse20=pulse20.*sqrt(1/M);
0071 L=length(C(1,1,:));
0072 
0073 if Para.N_R==1 && Para.N_T==1
0074     C=squeeze(C);
0075     H=fft(C,2*M);
0076     H1=-sqrt(-1)*fft((0:L-1)'.*C,2*M);
0077     H2=-fft((0:L-1).^2'.*C,2*M);
0078     
0079     switch criterion
0080         case 'ZF'
0081             B=1./H;
0082             B1=-H1./(H.^2);
0083             B2=-H2./(H.^2)+2*H1.^2./(H.^3);
0084         case 'MMSE'
0085             % Noise power
0086             Es_N0=10.^(Para.Es_N0_dB/10);
0087             B=conj(H)./(abs(H).^2+1/Es_N0);
0088             B1=(conj(H1).*1/Es_N0-H1.*conj(H).^2)./((abs(H).^2+1/Es_N0).^2);
0089             if R>1
0090                 error('Case R>1 not implemented yet for MMSE equalizer criterion')
0091             end
0092         case 'Specific'
0093             if exist('B','var')==0
0094                 error('B should be provided as argument for "Specific" criterion')
0095             end
0096             if R>1
0097                 if exist('B1','var')==0
0098                     error('B1 should be provided as argument for "Specific" criterion')
0099                 end
0100                 if R>2
0101                     if exist('B1','var')==0
0102                         error('B2 should be provided as argument for "Specific" criterion')
0103                     end
0104                 end
0105             end
0106         otherwise
0107             error('Equalizer criterion not existing')
0108     end
0109     z=AFB(r,pulse,Para);
0110     z_10=AFB(r,fliplr(pulse10),Para);
0111     z_20=AFB(r,fliplr(pulse20),Para);
0112     x=zeros(2*M,2*Para.Ns);
0113     for m=Para.ActiveSubcarriers
0114         x(m,:)=B(m)*z(m,:);
0115         if R>0
0116             x(m,:)=x(m,:)+(1j/(2*M))*B1(m)*z_10(m,:);
0117             if R>1
0118                 x(m,:)=x(m,:)-1/(2*(2*M)^2)*B2(m)*z_20(m,:);
0119             end
0120         end
0121     end
0122     x=x(:,2*Para.PreambleLength+1:end);
0123     x=x.*conj(theta);
0124     
0125 else
0126     if Para.N_T>Para.N_R
0127         error('Error: N_T>N_R')
0128     else       
0129         H=zeros(Para.N_R,Para.N_T,2*M);
0130         H_1=zeros(Para.N_R,Para.N_T,2*M);
0131         H_2=zeros(Para.N_R,Para.N_T,2*M);
0132         for i=1:Para.N_R
0133             for j=1:Para.N_T
0134                 H(i,j,:)=fft(squeeze(C(i,j,:)),2*M).';
0135                 H_1(i,j,:)=-sqrt(-1)*fft((0:L-1).*squeeze(C(i,j,:)).',2*M).';
0136                 H_2(i,j,:)=-fft((0:L-1).^2.*squeeze(C(i,j,:)).',2*M).';
0137             end
0138         end
0139         B=zeros(Para.N_T, Para.N_R, 2*M);
0140         B1=zeros(Para.N_T, Para.N_R, 2*M);
0141         B2=zeros(Para.N_T, Para.N_R, 2*M);
0142         switch criterion
0143             case 'ZF'
0144                 for m=Para.ActiveSubcarriers
0145                     H_m=H(:,:,m);
0146                     temp_inv=(H_m'*H_m)\H_m';
0147                     B(:,:,m)=temp_inv;
0148                     if R>0
0149                         B1(:,:,m)=-temp_inv*H_1(:,:,m)*temp_inv;
0150                         if R>1
0151                             B2(:,:,m)=2*temp_inv*H_1(:,:,m)*temp_inv*H_1(:,:,m)*temp_inv-temp_inv*H_2(:,:,m)*temp_inv;
0152                         end
0153                     end
0154                 end
0155             case 'MMSE'
0156                 % Noise power
0157                 Es_N0=10.^(Para.Es_N0_dB/10);
0158                 for m=Para.ActiveSubcarriers
0159                     H_m=H(:,:,m);
0160                     H1_m=H_1(:,:,m);
0161                     temp_inv=eye(Para.N_T,Para.N_T)/(H_m'*H_m+1/Es_N0*Para.N_T*eye(Para.N_T,Para.N_T));
0162                     B(:,:,m)=temp_inv*H_m';
0163                     if R>0
0164                         if R>0
0165                             B1(:,:,m)=-temp_inv*(H1_m'*H_m+H_m'*H1_m )*temp_inv*H_m'...
0166                                 +temp_inv*H1_m';
0167                             if R>1
0168                                 error('MMSE equalizer criterion not yet implemented for MIMO and R>0')
0169                             end
0170                         end
0171                     end
0172                 end
0173             case 'Specific'
0174                 error('Specific option not implemented for MIMO')
0175             otherwise
0176                 error('Equalizer criterion not existing')
0177         end
0178         z=zeros(Para.N_R,2*M,2*(Para.Ns+Para.PreambleLength));
0179         z_10=zeros(Para.N_R,2*M,2*(Para.Ns+Para.PreambleLength));
0180         z_20=zeros(Para.N_R,2*M,2*(Para.Ns+Para.PreambleLength));
0181         for index_N_R=1:Para.N_R
0182             z(index_N_R,:,:)=AFB(r(index_N_R,:),pulse,Para);
0183             z_10(index_N_R,:,:)=AFB(r(index_N_R,:),fliplr(pulse10),Para);
0184             z_20(index_N_R,:,:)=AFB(r(index_N_R,:),fliplr(pulse20),Para);
0185         end
0186         z=z(:,:,2*Para.PreambleLength+1:end);
0187         z_10=z_10(:,:,2*Para.PreambleLength+1:end);
0188         z_20=z_20(:,:,2*Para.PreambleLength+1:end);
0189         x=zeros(Para.N_T,2*M,2*Para.Ns);
0190         for m=Para.ActiveSubcarriers
0191             for l=1:2*Para.Ns
0192                 x(:,m,l)=B(:,:,m)*z(:,m,l);
0193             end
0194             if R>0
0195                 for l=1:2*Para.Ns
0196                     x(:,m,l)=x(:,m,l)+(1j/(2*M))*B1(:,:,m)*z_10(:,m,l);
0197                 end
0198                 if R>1
0199                     for l=1:2*Para.Ns
0200                         x(:,m,l)=x(:,m,l)-1/(2*(2*M)^2)*B2(:,:,m)*z_20(:,m,l);
0201                     end
0202                 end
0203             end
0204             
0205         end
0206         for index_N_T=1:Para.N_T
0207             x(index_N_T,:,:)=squeeze(x(index_N_T,:,:)).*conj(theta);
0208         end
0209         if Para.N_T==1
0210             x=squeeze(x);
0211         end
0212     end
0213     
0214 end
0215 
0216 end
0217 
0218 
0219 
0220 
0221 
0222 
0223 
0224 
0225 
0226 
0227 
0228 
0229 
0230 
0231 
0232 
0233