%%%%%%% MATLAB (R) script to compute TWDP %%%%%%%%%%%%%%%%%%%%%%
%% TP-2 test inputs
%% Units for all optical power values must match.

%PD The values given below for TxDataFile and MeasuredWaveformFile are
%% examples and should be replaced by actual path\filenames and values for each waveform tested.

%% Transmit data file: The transmit data sequence is based on either of the TWDP test patterns defined in
%% Table 68-5. The file format is ASCII with a single column of chronological ones and zeros with no headers
%% or footers.
TxDataFile = 'prbs9_ 950.txt';

%% Measured waveform. The waveform consists of exactly N samples per unit interval T, where N is the
%% oversampling rate. The waveform must be circularly shifted to align with the data sequence. The file
%% format for the measured waveform is ASCII with a single column of chronological numerical samples, in
%% optical power, with no headers or footers.
MeasuredWaveformFile = 'preproc-1207-01.txt';

%% OMA and steady-state ZERO power must be entered.
% MeasuredOMA         = [6.6e-004];       % Measured OMA, in optical power
% SteadyZeroPower     = [2.76e-004]; 		% Measured steady-state logic ZERO
OverSampleRate      = 16;           	% Oversampling rate

%% Simulated fiber response, modeled as a set of ideal delta functions with specified amplitudes in optical
%% power and delays in nanoseconds in columns with no headers or footers. The number of test cases is
%% determined by the TWDP requirements in Table 68-3. The vector 'PCoefs' contains the amplitudes, and the
%% vector 'Delays' contains the delays.
FiberResp = [...
    0.000000 0.65 0.88 0.51
    0.072727 0.5 0.58 0.89
    0.145455 0.91 0.89 0.29
    0.218182 0.26 0.1 0.81];
Delays          = FiberResp(:, 1);

%% Program constants %%
SymbolPeriod    = 1/(10.3125);          % Symbol period (ns)
EFilterBW       = 7.5;                  % Front end filter bandwidth (GHz)
EqNf            = 100;                  % Number of feedforward equalizer taps
EqNb            = 50;                   % Number of feedback equalizer taps
EqDel           = ceil(EqNf/2);     % Equalizer delay
PAlloc          = 6.5;              % Allocated dispersion penalty (dBo)
Q0              = 7.03;             % BER = 10^(-12)

%% STEP 1 - Process waveform through simulated fiber channel %%%%%%%%%%%%%%%%%%%
%% Load input waveforms
XmitData    = load(TxDataFile);
yout0       = load(MeasuredWaveformFile);
PtrnLength  = length(XmitData);
TotLen      = PtrnLength*OverSampleRate;
Fgrid       = [-TotLen/2:TotLen/2-1].'/(PtrnLength*SymbolPeriod);

%PD New section calculating signal strength from waveform
MeasuredSpectrum=abs(fft(yout0));
MeasuredMean=MeasuredSpectrum(1)/TotLen
% Find the "perfect" pattern, after NRZ encoding
for k=1:length(XmitData)
        TheoryWaveform(OverSampleRate*(k-1)+1:OverSampleRate*k,1)=XmitData(k);
end % k
% Find its spectrum
TheorySpectrum=abs(fft(TheoryWaveform));
% Bessel filter
y=2.114*j*Fgrid/7.5;
H_Bessel=abs(105./(105+105*y+45*y.^2+10*y.^3+y.^4));
H_Bessel=circshift(H_Bessel,TotLen/2);  % would fftshift have done?
TheorySpectrumF=abs(TheorySpectrum.*H_Bessel);
% Spectra are rough, so do some smoothing.  First replace the DC component
MeasuredSpectrum(1)=MeasuredSpectrum(2);
TheorySpectrumF(1)=TheorySpectrumF(2);
% Smoothing.  Temporarily rearrage so 0 Hz is in the middle
Swindow=19;  % an odd number
Svector=ones(Swindow,1)/Swindow;
MeasuredSpectrumS=filter(Svector,1,fftshift(MeasuredSpectrum));
MeasuredSpectrumS=circshift(MeasuredSpectrumS,-(Swindow-1)/2);
MeasuredSpectrumS=fftshift(MeasuredSpectrumS);
TheorySpectrumFS=filter(Svector,1,fftshift(TheorySpectrumF));
TheorySpectrumFS=circshift(TheorySpectrumFS,-(Swindow-1)/2);
TheorySpectrumFS=fftshift(TheorySpectrumFS);
% "MeasuredOMA" is not exactly that, should be renamed to e.g. WaveSize
% Look for max up to half the line rate - usually max would be very near DC
MeasuredOMA=max(MeasuredSpectrumS(2:256)./TheorySpectrumFS(2:256))
% Find the expected mean considering the pattern
DataMean=sum(XmitData)/length(XmitData)
SteadyZeroPower=MeasuredMean-MeasuredOMA*(1-DataMean)
% End of changes

%% Process through fiber model. Fiber frequency response is normalized to 1 at DC.
for i=1:3
    PCoefs      = FiberResp(:,i+1);
    ExpArg      = -j*2*pi*Fgrid;
    Hsys        = exp(ExpArg * Delays') * PCoefs;
    Hx          = fftshift(Hsys/sum(PCoefs));
    yout        = real(ifft(fft(yout0).*Hx));

    %% STEP 2 - Normalize OMA%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    yout        = (yout - SteadyZeroPower)/MeasuredOMA;

    %% STEP 3 - Process signal through front-end antialiasing filter %%%%%%%%%%%%%%%%%%
    %% Compute frequency response of front-end Butterworth filter
    [b,a]       = butter(4, 2*pi*EFilterBW,'s');
    H_r         = freqs(b,a,2*pi*Fgrid);
    %% Process signal through front-end filter
    yout        = real(ifft(fft(yout) .* fftshift(H_r)));

    %% STEP 4 - Sample at rate 2/T %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    yout        = yout(1:OverSampleRate/2:end);

    %% STEP 5 - Compute MMSE-DFE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% The MMSE-DFE filter coefficients computed below minimize mean-squared error at the slicer input.
    %% The derivation follows from the fact that the slicer input over one period (which is the same as the
    %% period of the input data sequence) can be expressed as Z = (R+N)*W - X*[0 B]', where R and N are
    %% Toeplitz matrices constructed from the signal and noise components, respectively, at the sampled output
    %% of the antialiasing filter, W is the feedforward filter, X is a Toeplitz matrix constructed from the
    %% input data sequence, and B is the feedback filter. The optimal W and B minimize the mean square error
    %% between the input to the slicer and the transmitted sequence due to residual ISI and Gaussian noise.
    %% Compute the noise autocorrelation sequence at the output of the front-end filter and rate-2/T sampler.
    N0   	= SymbolPeriod/(2 * Q0^2 * 10^(2*PAlloc/10));
    Snn 	= N0/2 * fftshift(abs(H_r).^2) * 1/SymbolPeriod * OverSampleRate;
    Rnn 	= real(ifft(Snn));
    Corr 	= Rnn(1:OverSampleRate/2:end);
    %% Construct a Toeplitz autocorrelation matrix.
    C     	= toeplitz(Corr(1:EqNf));
    %% Construct Toeplitz matrix from input data sequence
    X     	= toeplitz(XmitData, [XmitData(1); XmitData(end:-1:end-EqNb+1)]);
    %% Construct Toeplitz matrix from signal at output of 2/T sampler.
    %% This sequence gets wrapped by equalizer delay
    R     	= toeplitz(yout, [yout(1); yout(end:-1:end-EqNf+2)]);
    R   	= [R(EqDel+1:end,:); R(1:EqDel,:)];
    R   	= R(1:2:end, :);
    %% Compute least-squares solution for filter coefficients
    RINV 	= inv(R'*R+PtrnLength*C);
    P     	= X'*(eye(PtrnLength) - R*RINV*R')*X;
    P01 	= P(1,2:EqNb+1);
    P11 	= P(2:EqNb+1,2:EqNb+1);
    B   	= -inv(P11)*P01'; 	% Feedback filter
    W     	= RINV*R'*X*[1;B]; 	% Feedforward filter
    Z     	= R*W - X*[0;B]; 	% Input to slicer

    %% STEP 6 - Compute BER using semi-analytic method %%%%%%%%%%%%%%%%%%%%%%
    MseGaussian 	= W'*C*W;
    Ber             = sum(0.5*erfc((abs(Z-0.5)/sqrt(MseGaussian))/sqrt(2)))/length(Z);

    %% STEP 7 - Compute equivalent SNR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% This function computes the inverse of the Gaussian error probability function. The built-in
    %% function erfcinv() is not sensitive enough for the low probability of error case.
    Q 	= inf;
    if Ber>10^(-12)         Q 	= sqrt(2)*erfinv(1-2*Ber);
    elseif Ber>10^(-300) 	Q 	= 2.1143*(-1.0658-log10(Ber)).^0.5024;
    end

    %% STEP 8 - Compute penalty %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    RefSNR       	= 10 * log10(Q0) + PAlloc;
    TrialTWDP(i) 	= RefSNR-10*log10(Q);
end

%% Pick highest value due to the multiple responses from TrialTWDP.
TWDP    = max(TrialTWDP)	% End of program