close all;
format short;

OutputPort1     = InputPort1;
OutputPort2     = InputPort2;
OutputPort3     = InputPort3;
OutputPort4     = InputPort4;

%Input ports:
%Input I-X: Electrical;
%Input Q-X: Electrical;
%Input I-Y: Electrical;
%Input Q-Y: Electrical;
%Output Ports:
%Output I-X: Electrical;
%Output Q-X: Electrical;
%Output I-Y: Electrical;
%Output Q-Y: Electrical;
%



Reference_Wavelength        =  Parameter0;
RefWavelength  =  Parameter1;
Dispersion   =   Parameter2;
DispersionSlope  =  Parameter3;
ChannelWavelength =  Parameter4;
PropagationLength =  Parameter5;

%%%%%%%%%%%%%%%%%%%%%%% Parameters Definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Reference_Wavelength  =  Reference_Wavelength  *1e-9;
RefWavelength =  RefWavelength * 1e-9;
Dispersion =  Dispersion * 1e-12 / 1e-6;
DispersionSlope =  DispersionSlope * 1e-12 / 1e-15;
ChannelWavelength = ChannelWavelength* 1e-9;
PropagationLength =PropagationLength* 1e3 ;
c0 =   2.99792458e8;
NumberOfSymbols =  2048;

SamplesPerSymbol =  2;
SymbolRate=  25e9;

StepP= 0.11; %(Control gain for channel estimation algorithm)
Initial Value= 1 ;                %(Initial value for the channel estimation)
TapNumberP= 3;                    %(Order of the FIR filter)

PMDbit=0;
%JP-CPE (joint polarization carrier phase estimation)
Coupling Factor = 0; %(Coupling factor between polarization X,Y)
TapNumberC = 3; %(Filter order for MMSE phase estimation)
RotationFactorX = 0; %(Rotating I and Q for X polarization)
RotationFactorY = 0; %(Rotating I and Q for Y polarization)
%ConstellationDiagram:
IgnoreStartSymbols=100;
IgnoreEndSymbols= 100;
%InterpolationMethod
InterpolationMethod= Cubic;

%%%%%%%%%%%%%%%%%%%%%%%%% Down Sampling %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Input ports:
if(InputPort1.TypeSignal == 'Electrical')
    (InputPort2.TypeSignal == 'Electrical')
    (InputPort3.TypeSignal == 'Electrical')
    (InputPort4.TypeSignal == 'Electrical')
%In Q-X : Electrical;
%In I-Y : Electrical;
%In Q-Y : Electrical;
%Output Ports:
%Out I-X: Electrical;
%Out Q-X: Electrical;
%Out I-Y: Electrical;
%Out Q-Y: Electrical;
SamplesPerSymbolIn  =   length(InputPort1.Sampled.Signal) / NumberOfSymbols;                              % calculates number of samples per symbol

InX_Signal = ones(1, NumberOfSymbols * SamplesPerSymbol);
InX_Noise = ones(1, NumberOfSymbols * SamplesPerSymbol);
InY_Signal = ones(1, NumberOfSymbols * SamplesPerSymbol);
InY_Noise = ones(1, NumberOfSymbols * SamplesPerSymbol);

%%%%%%%%%%%%%%%%%%%%%  module activation%%%%%%%%%%%%%
EnableC =1;
EnableD =1;
EnableP =1;

for m = 1 : NumberOfSymbols * SamplesPerSymbol
    
    InX_Signal(1,m) =   InputPort1.Sampled.Signal((2*m-1) * SamplesPerSymbolIn / (2*SamplesPerSymbol)) + 1i * InputPort2.Sampled.Signal((2*m-1) * SamplesPerSymbolIn / (2*SamplesPerSymbol));
    InX_Noise(1,m)  =   InputPort1.Noise.Signal((2*m-1) * SamplesPerSymbolIn / (2*SamplesPerSymbol)) + 1i * InputPort2.Noise.Signal((2*m-1) * SamplesPerSymbolIn / (2*SamplesPerSymbol));
    
    InX             =   InX_Signal + InX_Noise;
    
    InY_Signal(1,m) =   InputPort3.Sampled.Signal((2*m-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol))+1i*InputPort4.Sampled.Signal((2*m-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol));
    InY_Noise(1,m)  =   InputPort3.Noise.Signal((2*m-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol))+1i*InputPort4.Noise.Signal((2*m-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol));
    
    InY             =   InY_Signal + InY_Noise;
end;

clear InX_Signal;
clear InX_Noise;
clear InY_Signal;
clear InY_Noise;


%In_X: signal and noise for polarization X
%In_Y: signal and noise for polarization Y

SampledVector       =   (2*(1:NumberOfSymbols*SamplesPerSymbol)-1)*SamplesPerSymbolIn/(2*SamplesPerSymbol);
TimeSampledVector   =   InputPort1.Sampled.Time(SampledVector);


%%%%%%%%%%%%%% Normalizing %%%%%%%%%%%%%%%%%%%%%%
OutX        = InX/sqrt(mean(abs(InX).^2));
OutY        = InY/sqrt(mean(abs(InY).^2));


StoreInX            =   OutX;
StoreInY            =   OutY;


clear SampledVector;
%%%%%%%%%%%%%%%%%%%%%% Dispersion Compensation %%%%%%%%%%%%%%%%%%%

if EnableD == 1
    InX         = OutX;
    InY         = OutY;
    L           = PropagationLength * 1e03;
    S           = DispersionSlope;
    D           = Dispersion + (ChannelWavelength - RefWavelength) * S;
    
    NumD_Eff	= 9999;                                                         % Odd number
    DTime       = (1 / (SymbolRate)) / SamplesPerSymbol;
    DFreq       = (1 / DTime) / (NumD_Eff - 1);
    VFreq       = (-(NumD_Eff-1) / 2 : (NumD_Eff-1) / 2) * DFreq;
    VOmeg       = 2 * pi * VFreq;                                               % Angular frequency 
    
    Hf          = exp(-1i * (VOmeg .^ 2) * ( (ChannelWavelength ^ 2 * D * L) / (4 * pi * c0) ));% + (1i * ( (S * L * ChannelWavelength ^ 4 * VOmeg .^ 3) / (24 * pi ^ 2 * c0 ^ 2) )) );
    Ht          = (fftshift(ifft(ifftshift((Hf)))));
    Htc         = conj(Ht);

    B = SymbolRate * 2;
    TapNumberD = 2 * floor((abs(D) * (ChannelWavelength ^ 2) * L * (B ^ 2)) / (2 * c0)) + 1;
    if mod(TapNumberD,2) == 0
        TapNumberD = TapNumberD + 1;
    end;
    
    TapD        = Htc(1,((NumD_Eff+1)/2-(TapNumberD-1)/2):((NumD_Eff+1)/2+(TapNumberD-1)/2));

    OutX_Temp   = conv(InX,TapD);           
    conL = length(OutX_Temp);
    
    OutX        = OutX_Temp(1,((TapNumberD-1)/2) + 1 : end - ((TapNumberD-1)/2));
    
    OutY_Temp   = conv(InY,TapD);           
    OutY        = OutY_Temp(1,((TapNumberD-1)/2) + 1 : end - ((TapNumberD-1)/2));

    AfterDX     = OutX;
    AfterDY     = OutY;

end

clear OutX_Temp OutY_Temp;
clear Hf Ht Htc TapD;
clear DTime DFreq VFreq VOmeg;

%%%%%%%%%%%%%%%%%%%%%% Polarization Demultiplexing %%%%%%%%%%%%%%%%%%%%%%%%%%
if EnableP==1
    
     InX     = circshift(OutX',0)';%OutX;
     InY     = circshift(OutY',0)';
    
    clear TempInX TempInY;
     
    Hxx     = zeros(1,TapNumberP);
    Hxy     = zeros(1,TapNumberP);
    Hyx     = zeros(1,TapNumberP);
    Hyy     = zeros(1,TapNumberP); 
    
% Hxx = [0 0 0 0 0 0 1 0 0 0 0 0 0];
% Hyy = [0 0 0 0 0 0 1 0 0 0 0 0 0];

    
    if InitialValue == 1
    Hxx((TapNumberP+1)/2) = 1;
   Hyy((TapNumberP+1)/2) = 1;
   Hxy((TapNumberP+1)/2) = 0;
   Hyx((TapNumberP+1)/2) = 0;
   elseif InitialValue == 2
        Hxx((TapNumberP+1)/2) = 0;
   Hyy((TapNumberP+1)/2) = 0;
   Hxy((TapNumberP+1)/2) = 1;
   Hyx((TapNumberP+1)/2) = 1; 
      elseif InitialValue == 3
       Hxx((TapNumberP+1)/2) = 1;
   Hyy((TapNumberP+1)/2) = 0;
   Hxy((TapNumberP+1)/2) = 0;
   Hyx((TapNumberP+1)/2) = 1; 
   elseif InitialValue == 4
       Hxx((TapNumberP+1)/2) = 0;
   Hyy((TapNumberP+1)/2) = 1;
   Hxy((TapNumberP+1)/2) = 1;
   Hyx((TapNumberP+1)/2) = 0; 
   end
    
    OutX = zeros(1,length(InX));
    OutY = zeros(1,length(InY));
   
    for m = 1 : length(InX) - TapNumberP + 1
        
           if mod(m,2) == PMDBit    
        OutX(m)  = sum(Hxx .* InX(m : m + TapNumberP - 1)) + sum(Hxy .* InY(m : m + TapNumberP - 1));       
        OutY(m)  = sum(Hyx .* InX(m : m + TapNumberP - 1)) + sum(Hyy .* InY(m : m + TapNumberP - 1));

        Ex(m) = (1 - (abs(OutX(m)) ^ 2));
        Ey(m) = (1 - (abs(OutY(m)) ^ 2));
              
        Hxx     = Hxx + StepP *  (1 - (abs(OutX(m)) ^ 2)) * (OutX(m) * conj(InX(m : m + TapNumberP - 1)));
        
        Hxy     = Hxy + StepP *  (1 - (abs(OutX(m)) ^ 2)) * (OutX(m) * conj(InY(m : m + TapNumberP - 1)));

        Hyx     = Hyx + StepP *  (1 - (abs(OutY(m)) ^ 2)) * (OutY(m) * conj(InX(m : m + TapNumberP - 1)));

        Hyy     = Hyy + StepP *  (1 - (abs(OutY(m)) ^ 2)) * (OutY(m) * conj(InY(m : m + TapNumberP - 1)));  
           end
    end;        
       
if PMDBit == 0
   PMDBit = 1; 
end;
    
%     
%     OutX_Temp    = conv(Hxx,InX(PMDBit : 2 : end)) +  conv(Hxy,InY(PMDBit : 2 : end));
%     OutY_Temp    = conv(Hyx,InX(PMDBit : 2 : end)) +  conv(Hyy,InY(PMDBit : 2 : end)); 
    
    OutX_Temp    = conv(Hxx,InX) +  conv(Hxy,InY);
    OutY_Temp    = conv(Hyx,InX) +  conv(Hyy,InY); 
        
%     filterCenterTap_X = find(abs(Hxx) == max(abs(Hxx)));
%     filterCenterTap_Y = find(abs(Hyy) == max(abs(Hyy)));

% OutX_Temp = OutX(PMDBit : 2 : end);
% OutY_Temp = OutY(PMDBit : 2 : end);
    
%     OutX = OutX_Temp(filterCenterTap_X : end - (TapNumberP - filterCenterTap_X + 0));
%  OutX = OutX(PMDBit:2:end); 

%     OutY = OutY_Temp(filterCenterTap_Y : end - (TapNumberP - filterCenterTap_Y + 0));
%  OutY = OutY(PMDBit:2:end); 
    
   OutX        = OutX_Temp(1,((TapNumberP-1)/2) + 1 : end - ((TapNumberP-1)/2));
   OutX        = OutX(1 : 2 : end);
%    OutX        = circshift(OutX',-1)';
    
   OutY        = OutY_Temp(1,((TapNumberP-1)/2) + 1 : end - ((TapNumberP-1)/2));
   OutY        = OutY(1 : 2 : end);
%    OutY        = circshift(OutY',0)';
%     
% OutX = OutX_Temp;
% OutY = OutY_Temp;    

AfterPX     = OutX;
AfterPY     = OutY;
    
         SamplesPerSymbol = SamplesPerSymbol / 2;
        TimeSampledVector = TimeSampledVector(PMDBit : 2 : end);


end

% clear Hxx Hxy Hyx Hyy;

%%%%%%%%%%%%%%%%%%%%%%%%% Joint Polarization Carrier Phase Estimation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if EnableC==1
    InX     = OutX;
    InY     = OutY;
    
    clear TempInX;
    clear TempInY;
    
    TempInX = [InX zeros(1,TapNumberC - 1)] * exp(-1i*pi/4);
    TempInY = [InY zeros(1,TapNumberC - 1)] * exp(-1i*pi/4);
    
    
    PhaseTemp = zeros(1,NumberOfSymbols * SamplesPerSymbol);
    PhaseErr = zeros(1,NumberOfSymbols * SamplesPerSymbol + 1);
   
    for m = 1 : NumberOfSymbols * SamplesPerSymbol
        
        FlagX       = sum(TempInX(m : (m + TapNumberC - 1)) .^ 4);
        FlagY       = sum(TempInY(m : (m + TapNumberC - 1)) .^ 4);
        
        Flag = 0.5 * (FlagX + FlagY);
              
        PhaseTemp(m) = (1/4) * angle(Flag);
        
        PhaseErr(m+1) =  PhaseTemp(m) + floor(0.5 + ((PhaseErr(m) -  PhaseTemp(m))) * (2 / pi)) * (pi / 2);
       
    end



Phase = PhaseErr(2 : end);
       
    OutX    = exp(-1i * Phase) .* InX * exp(1i * RotationFactorX * pi);
    OutY    = exp(-1i * Phase) .* InY * exp(1i * RotationFactorY * pi);
    
  
    AfterCX     = OutX;
    AfterCY     = OutY;
end

clear TempInX TempInY;
clear PhaseX PhaseY FlagX FlagY;
clear Flag PhaseTemp PhaseErr Phase;
%%%%%%%%%%%%%%%%%%%%%%%%% Plot constellation diagrams %%%%%%%%%%%%%%%%%%%%%%%%%
NumCount = 1+EnableD+EnableP+EnableC;
figure;
h1=subplot(2,NumCount,1);
scatter(h1,real(StoreInX((IgnoreStartSymbols*2+1):2:(length(OutX)-IgnoreEndSymbols*2))),imag(StoreInX((IgnoreStartSymbols*2+1):2:(length(OutX)-IgnoreEndSymbols*2))),'r*')
xlabel('In-Phase'); ylabel('Quadrature'); grid on; axis([-1.2 1.2 -1.2 1.2]); title('Before DSP - X','FontSize',14)
h5=subplot(2,NumCount,1+NumCount);
scatter(h5,real(StoreInY((IgnoreStartSymbols*2+1):2:(length(OutX)-IgnoreEndSymbols*2))),imag(StoreInY((IgnoreStartSymbols*2+1):2:(length(OutX)-IgnoreEndSymbols*2))),'r*')
xlabel('In-Phase'); ylabel('Quadrature'); grid on; axis([-1.2 1.2 -1.2 1.2]); title('Before DSP - Y','FontSize',14)

if EnableD==1
    h2=subplot(2,NumCount,1+EnableD);
    scatter(h2,real(AfterDX((IgnoreStartSymbols*2+1):2:(length(AfterDX)-IgnoreEndSymbols*2))),imag(AfterDX((IgnoreStartSymbols*2+1):2:(length(AfterDX)-IgnoreEndSymbols*2))),'b*')
    xlabel('In-Phase'); ylabel('Quadrature'); grid on; axis([-1.2 1.2 -1.2 1.2]); title('After Dispersion Compensation - X','FontSize',14)
    h6=subplot(2,NumCount,1+EnableD+NumCount);
    scatter(h6,real(AfterDY((IgnoreStartSymbols*2+1):2:(length(AfterDY)-IgnoreEndSymbols*2))),imag(AfterDY((IgnoreStartSymbols*2+1):2:(length(AfterDY)-IgnoreEndSymbols*2))),'b*')
    xlabel('In-Phase'); ylabel('Quadrature'); grid on; axis([-1.2 1.2 -1.2 1.2]); title('After Dispersion Compensation - Y','FontSize',14)

end

if EnableP==1
    h3=subplot(2,NumCount,1+EnableD+EnableP);
    scatter(h3,real(AfterPX((IgnoreStartSymbols*2+1):2:(length(AfterPX)-IgnoreEndSymbols*2))),imag(AfterPX((IgnoreStartSymbols*2+1):2:(length(AfterPX)-IgnoreEndSymbols*2))),'b*')
    xlabel('In-Phase'); ylabel('Quadrature'); grid on; axis([-1.2 1.2 -1.2 1.2]); title('After Polarization Demultiplexing - X','FontSize',14)
    h7=subplot(2,NumCount,1+EnableD+EnableP+NumCount);
    scatter(h7,real(AfterPY((IgnoreStartSymbols*2+1):2:(length(AfterPY)-IgnoreEndSymbols*2))),imag(AfterPY((IgnoreStartSymbols*2+1):2:(length(AfterPY)-IgnoreEndSymbols*2))),'b*')
    xlabel('In-Phase'); ylabel('Quadrature'); grid on; axis([-1.2 1.2 -1.2 1.2]); title('After Polarization Demultiplexing - Y','FontSize',14)
end

if EnableC==1
    h4=subplot(2,NumCount,1+EnableD+EnableP+EnableC);
    scatter(h4,real(AfterCX((IgnoreStartSymbols*2+1):2:(length(AfterCX)-IgnoreEndSymbols*2))),imag(AfterCX((IgnoreStartSymbols*2+1):2:(length(AfterCX)-IgnoreEndSymbols*2))),'b*')
    xlabel('In-Phase'); ylabel('Quadrature'); grid on; axis([-1.2 1.2 -1.2 1.2]); title('After Carrier Phase Estimation - X','FontSize',14)
    h8=subplot(2,NumCount,1+EnableD+EnableP+EnableC+NumCount);
    scatter(h8,real(AfterCY((IgnoreStartSymbols*2+1):2:(length(AfterCY)-IgnoreEndSymbols*2))),imag(AfterCY((IgnoreStartSymbols*2+1):2:(length(AfterCY)-IgnoreEndSymbols*2))),'b*')
    xlabel('In-Phase'); ylabel('Quadrature'); grid on; axis([-1.2 1.2 -1.2 1.2]); title('After Carrier Phase Estimation - Y','FontSize',14)

end

clear StoreInX StoreInY;
clear AfterDX AfterDY;
clear AfterPX AfterPY;
clear AfterCX AfterCY;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Interpolation and Output to OptiSystem %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OutXT = zeros(1, NumberOfSymbols * SamplesPerSymbol * SamplesPerSymbolIn);
% OutYT = zeros(1, NumberOfSymbols * SamplesPerSymbol * SamplesPerSymbolIn);
% 
% 
% for i = 1 : NumberOfSymbols * SamplesPerSymbol
%     
%     for m = SamplesPerSymbolIn * (i - 1) + 1 : SamplesPerSymbolIn * i
%         OutXT(m) = OutX(i);
%         OutYT(m) = OutY(i);
%     end;
%     
% end;
% 
% OutX = OutXT;
% OutY = OutYT;
% 
% OutputPort1.Sampled.Signal  = real(OutX); 
% OutputPort2.Sampled.Signal  = imag(OutX);
% OutputPort3.Sampled.Signal  = real(OutY);
% OutputPort4.Sampled.Signal  = imag(OutY);
% 
OutputPort1.Sampled.Signal  = interp1(TimeSampledVector,real(OutX),InputPort1.Sampled.Time,InterpolationMethod); 
OutputPort2.Sampled.Signal  = interp1(TimeSampledVector,imag(OutX),InputPort1.Sampled.Time,InterpolationMethod);
OutputPort3.Sampled.Signal  = interp1(TimeSampledVector,real(OutY),InputPort1.Sampled.Time,InterpolationMethod);
OutputPort4.Sampled.Signal  = interp1(TimeSampledVector,imag(OutY),InputPort1.Sampled.Time,InterpolationMethod);

clear TimeSampledVector;
clear InX InY;
clear OutX OutY;
clear InputPort1 InputPort2 InputPort3 InputPort4;
end