FBMC_OQAM_PhaseTrackerMLMAP

 Tracks and corrects phase noise following MAP/ML algorithm described in [1].

 function [x, phi]=FBMC_OQAM_PhaseTrackerMLMAP( z, linewidth, type, Para )

 The function works for SISO systems.

 Input arguments:

   z: demodulated FBMC-OQAM symbols. Size: matrix [Para.nSubcarriers,2*Para.Ns]

   linewidth: phase noise linewidth [Hz]

   type: determines the type of tracker to use [1]. Several
 options are possible:

       'ML': maximum likelihood.

       'MAP': maximum a posteriori.

   Para: structure containing the modulation parameters.

 Ouput arguments:

   x: equalized symbols. Size: matrix [Para.nSubcarriers, 2*Para.Ns].

   phi: estimated phase. Size: vector [2*Para.Ns, 1].

 References: [1] F. Rottenberg, T.-H. Nguyen, S.-P. Gorza, J. Louveaux
   and F. Horlin, "ML and MAP phase noise estimators for optical fiber
   FBMC-OQAM systems." In IEEE International Conference on Communications
   (ICC), 2017, Paris.

0001 function [x, phi]=FBMC_OQAM_PhaseTrackerMLMAP( z, linewidth, type, Para )
0002 % Tracks and corrects phase noise following MAP/ML algorithm described in [1].
0003 %
0004 % function [x, phi]=FBMC_OQAM_PhaseTrackerMLMAP( z, linewidth, type, Para )
0005 %
0006 % The function works for SISO systems.
0007 %
0008 % Input arguments:
0009 %
0010 %   z: demodulated FBMC-OQAM symbols. Size: matrix [Para.nSubcarriers,2*Para.Ns]
0011 %
0012 %   linewidth: phase noise linewidth [Hz]
0013 %
0014 %   type: determines the type of tracker to use [1]. Several
0015 % options are possible:
0016 %
0017 %       'ML': maximum likelihood.
0018 %
0019 %       'MAP': maximum a posteriori.
0020 %
0021 %   Para: structure containing the modulation parameters.
0022 %
0023 % Ouput arguments:
0024 %
0025 %   x: equalized symbols. Size: matrix [Para.nSubcarriers, 2*Para.Ns].
0026 %
0027 %   phi: estimated phase. Size: vector [2*Para.Ns, 1].
0028 %
0029 % References: [1] F. Rottenberg, T.-H. Nguyen, S.-P. Gorza, J. Louveaux
0030 %   and F. Horlin, "ML and MAP phase noise estimators for optical fiber
0031 %   FBMC-OQAM systems." In IEEE International Conference on Communications
0032 %   (ICC), 2017, Paris.
0033 %
0034 
0035 % This file is part of WaveComBox: www.wavecombox.com and is distributed under the terms of the MIT license. See accompanying LICENSE file.
0036 % Original author: François Rottenberg, May 4, 2018.
0037 % Contributors:
0038 % Change log:
0039 
0040 M=Para.nSubcarriers/2;
0041 Ns=Para.Ns;
0042 M_PAM=Para.M_PAM;
0043 mod_norm=sqrt(Para.Es/2/(1./Para.M_PAM.*sum(abs(pammod([0:Para.M_PAM-1],Para.M_PAM)).^2)))*1./sqrt(Para.S);
0044 Es_N0=10.^(Para.Es_N0_dB/10);
0045 N0=Para.Es./Es_N0;
0046 sigma_epsilon_2=2*pi* linewidth *Para.T;
0047 sigma_u_2 =Para.Es/2;
0048 
0049 phi=zeros(2*Ns,1);
0050 phi_l_1=0;
0051 if Para.N_R==1 && Para.N_T==1
0052     x=zeros(2*M,2*Para.Ns);
0053     for l=1:2*Ns
0054         
0055         %% MAP without u_1
0056         z_tilde=exp(-1j*phi_l_1).*z(:,l);
0057         temp=real(z_tilde);
0058         d_hard_decision=pammod(pamdemod(temp/mod_norm,M_PAM),M_PAM).*mod_norm;
0059         
0060         switch type
0061             case 'MAP'
0062                 epsilon=sum((N0/2+sigma_u_2).*d_hard_decision.*imag(z_tilde)-sigma_u_2.*imag(z_tilde).*real(z_tilde))/(sum(sigma_u_2.*imag(z_tilde).^2)+(N0/2*(N0/2+sigma_u_2))/(sigma_epsilon_2));
0063             case 'ML'
0064                 epsilon=sum((N0/2+sigma_u_2).*d_hard_decision.*imag(z_tilde)-sigma_u_2.*imag(z_tilde).*real(z_tilde))/(sum(sigma_u_2.*imag(z_tilde).^2));
0065             otherwise
0066                 error('Type not corrected')
0067         end
0068         phi_l=phi_l_1+epsilon;
0069         phi_l_1=phi_l;
0070         phi(l)=phi_l;
0071         
0072         %% Re-estimation of the phase
0073         x(:,l)=exp(-1j*phi(l)).*z(:,l);
0074     end
0075 else
0076     if Para.S>2
0077         error('Function not implemented for S>2')
0078     end
0079     x=zeros(Para.S,2*M,2*Para.Ns);
0080     z_phase_track=[squeeze(z(1,:,:));squeeze(z(2,:,:))];
0081     for l=1:2*Ns
0082         z_tilde=exp(-1j*phi_l_1).*z_phase_track(:,l);
0083         temp=real(z_tilde);
0084         d_hard_decision=pammod(pamdemod(temp/mod_norm,M_PAM),M_PAM).*mod_norm;
0085         switch type
0086             case 'MAP'
0087                 epsilon=sum((N0/2+sigma_u_2).*d_hard_decision.*imag(z_tilde)-sigma_u_2.*imag(z_tilde).*real(z_tilde))/(sum(sigma_u_2.*imag(z_tilde).^2)+(N0/2*(N0/2+sigma_u_2))/(sigma_epsilon_2));
0088             case 'ML'
0089                 epsilon=sum((N0/2+sigma_u_2).*d_hard_decision.*imag(z_tilde)-sigma_u_2.*imag(z_tilde).*real(z_tilde))/(sum(sigma_u_2.*imag(z_tilde).^2));
0090             otherwise
0091                 error('Type not corrected')
0092         end
0093         phi_l=phi_l_1+epsilon;
0094         phi_l_1=phi_l;
0095         phi(l)=phi_l;
0096         % Correction with re-estimated current phase
0097         x(1,:,l)=exp(-1j*phi_l).*z(1,:,l);
0098         x(2,:,l)=exp(-1j*phi_l).*z(2,:,l);
0099     end
0100 end
0101 
0102 end
0103