%%%  BER Analyzer code is not 100% same as OptiSystem's  %%%

%%% Gathering data from input port %%%
binary = InputPort1.Sequence;
electrical = InputPort2.Sampled.Signal + InputPort2.Noise.Signal;
time = InputPort2.Sampled.Time;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use this if you wish to load the data from a file                                             %
% load('C:\Documents and Settings\bryan\Desktop\Junk Folder\BER\binary.mat');                   %
% load('C:\Documents and Settings\bryan\Desktop\Junk Folder\BER\electrical.mat');               %
% load('C:\Documents and Settings\bryan\Desktop\Junk Folder\BER\time.mat');                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% timeSpace is the difference in time between samples %%%
timeSpace = time(2) - time(1);
samples = length(electrical)/length(binary);%%samples per bit

%%% assign initial values and threshold level (currently 50%) %%%
timeArray(1) = 0;
spaceIndex = 0;
markIndex = 0;
elecIndex = 0;
threshold = 0.5;

%%% finds all 0's and 1's, and groups all their samples into separate arrays %%%
for i = 1:length(binary) 
    if (binary(i) == 0)
        spaceIndex = spaceIndex + 1;
        for j = 1:samples 
            spaceArray(spaceIndex, j) = electrical(elecIndex+j)
        end;
    else
        markIndex = markIndex + 1;
        for j = 1:samples
            markArray(markIndex, j) = electrical(elecIndex+j);
        end;  
    end;
    elecIndex = elecIndex + samples;
end;


%%% creates the timeArray or x-axis of the graphs %%%
for i = 2:samples 
    timeArray(i) = timeArray(i-1) + timeSpace;
end;
        
%%% plots the Eye Diagram %%%
figure
plot(timeArray, markArray, timeArray, spaceArray); 
title('Eye Diagram','FontSize',16);
pause(3);

%%% calculates average amplitude of 1's and 0's, their %%%
%%% standard deviations, Q factor, and eye height      %%% 
for j = 1:samples
    for i = 1:markIndex
      temp1(i) = markArray(i,j);  
    end
    u1(j) = sum(temp1)/markIndex;
    std1(j) = std(temp1);
    for i = 1:spaceIndex
      temp0(i) = spaceArray(i,j);  
    end;
    u0(j) = sum(temp0)/spaceIndex;
    std0(j) = std(temp0);
    if (std1(j)+std0(j) == 0)
        Q(j) = 0;
    else
        Q(j) = abs(u1(j)-u0(j))/(std1(j)+std0(j));
        eyeHeight(j) = (u1(j)-3*std1(j)^2) - (u0(j)+3*std0(j)^2);      
    end;
end

%%% plots the Q-factor %%%
figure
plot(timeArray, Q);
title('Q-factor','FontSize',16);
pause(3);

%%% plots the Eye Height%%%
figure
plot(timeArray, eyeHeight);
title('Eye Height','FontSize',16);
pause(3);
        
%%% calculates the threshold, probability of 1's and 0's, and the BER %%%
for j = 1:samples
    S(j) = (mean(markArray(:,j))-std1(j) + mean(spaceArray(:,j))+std0(j))/2;
    if (std0(j)==0)
       Pe0(j) = 0;
    else
       Pe0(j) = 1/2*erfc(abs(((S(j)-u0(j))/(sqrt(2)*std0(j))))); 
    end;
    if (std1(j)==0)
       Pe1(j) = 0;
    else
       Pe1(j) = 1/2*erfc(abs(((u1(j)-S(j))/(sqrt(2)*std1(j))))); 
    end;
    Pe(j) = log10(spaceIndex/(spaceIndex+markIndex)*Pe0(j) + markIndex/(markIndex+spaceIndex)*Pe1(j)); 
end

%%% plots the threshold %%%
figure
plot(timeArray, S);
title('Threshold','FontSize',16);
pause(3);

%%% plots the BER using Gaussian approx. %%%
figure
plot(timeArray, Pe);
title('Log of Min BER using Gaussian approx.','FontSize',16);
pause(3);

%%% calculates and plots the BER from the Q-factor %%%
PeWC = log10(1/2*erfc(Q/sqrt(2)));
figure
plot(timeArray, PeWC);
title('Log of Min BER from Q','FontSize',16);
pause(3);




        
        
        
        
        
        
        