Matlab code follows:

%%% 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;

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

%%% 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);