%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% TP2_TEST.M - MATLAB (R) script to compute TP-2 penalty
%%%%%%%
%%%%%%% Written and submitted to IEEE 802.3aq TP-2 by ClariPhy Communications, Inc.
%%%%%%% Refer to ClariPhy document "Proposed TP-2 Test Methodology" for further 
%%%%%%% description of the processing steps contained in this program.
%%%%%%%
%%%%%%% VERSION HISTORY
%%%%%%% 1/12/05 v0.1 - Initial
%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TP-2 test inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Transmit data file. The current sequence is for demonstration purposes
%% and will be replaced by a data sequence TBD by the LRM committee 
%% (In the case used here, the transmit data sequence 
%% is derived from a PRBS9 data sequence generated by polynomial x^9+x^4+1.
%% (That is, the data sequence d(n) is given by d(n)=d(n-9)+d(n-4), mod 2.)
%% The sequence is initially aligned so that it starts with 9 ones.  A zero is 
%% inserted immediately after the string of eight zeros in the sequence. 2 bits 
%% are inverted at location 391 and 392. The resulting data sequence is
%% inverted. Then the entire sequence is circularly right-shifted 476 bits
%% to align with the measured waveform specified below.)
TxDataFile           = 'TxData.txt';

%% Measured waveform. The current waveform is for demonstration purposes
%% only. The waveform consists of exactly N samples per bit
%% period T, where N is the oversampling rate. 
%% (For the PRBS9 case used here, the oversampling rate is 16 and the waveform 
%% consists of 512*16 = 8192 samples. The data sequence described above has been 
%% aligned with the waveform.) OMA and steady-state ZERO power must also be specified.
MeasuredWaveformFile = 'G05.txt';   % Measured waveform samples, in optical power
MeasuredOMA          = 3.8e-004; % Measured OMA, in optical power
SteadyZeroPower      = 3.2e-004; % Measured optical power, steady-state logic ZERO 
OverSampleRate       = 16;          % Oversampling rate

%% Simulated fiber response, modeled as a set of ideal delta functions with
%% specified amplitudes and delays.  The vector 'PCoefs' contains the amplitudes, 
%% and the vector 'Delays' contains the delays.  (This particular case is
%% taken from Cambridge model release 2.1, Fiber number = 1, Offset = 17
%% microns. This fiber is for demonstration purposes and will be replaced
%% by a set of channels TBD by the LRM committee)
FiberResp   = load('cam2p1_1_17.txt');
PCoefs      = FiberResp(:, 1);
Delays      = FiberResp(:, 2); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Static parameters that should not change once test is specified
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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);
yout       = load(MeasuredWaveformFile);
PRBSlength = length(XmitData);
TotLen     = PRBSlength*OverSampleRate;
Fgrid      = [-TotLen/2:TotLen/2-1].'/(PRBSlength*SymbolPeriod);

%% Process through fiber model.  Fiber frequency response is normalized 
%% to 1 at DC
ExpArg    = -j*2*pi*Fgrid;
Hsys      = exp(ExpArg * Delays') * PCoefs; 
Hx        = fftshift(Hsys/abs(Hsys(find(Fgrid==0)))); 
yout      = real(ifft(fft(yout).*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 E[||Z-XIn||^2], where XIn is the input data sequence, 
%% and the expectation operator refers to the Gaussian noise component of Z.

%% Compute the noise autocorrelation sequence at the output of the front-end 
%% filter and rate-2/T sampler.  Constuct a toeplitz autocorrelation matrix.
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);
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+PRBSlength*C);
P    = X'*(eye(PRBSlength) - 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 Matlab 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;
fprintf(1,'TP2 penalty equals %5.4f dB\n', RefSNR-10*log10(Q));