% Post processing MATLAB code for 1000BASE-RH transmitter distortion

% Read captured data file, min 13 us long, 3.25GSample/sec, 
% high resolution capture
fid  = fopen('RawData.bin', 'r');
xcap = fread(fid, inf, 'double');
fclose(fid);
xcap = reshape(xcap, 1, []);

% TX distortion evaluation
[HD2 HD3 RPD] = txdist(xcap, 10);

% Results
fprintf('HD2: %.2f dB\n', HD2);
fprintf('HD3: %.2f dB\n', HD3);
fprintf('RPD: %.2f dB\n', RPD);


%------------------------------------------------------------------------------
% TX distortion measurement
%------------------------------------------------------------------------------
function [HD2 HD3 RPD] = txdist(xcap, ov)
%----------------
% Compensate baseline
%----------------
xcap = xcap - mean(xcap);

%----------------
% Low pass filter at 162.5 MHz
%----------------
[hb, ha] = butter(2, 1/ov, 'low');
xcap = filter(hb, ha, xcap);

%----------------
% Reference test mode 6 signal
%----------------
tm6 = tm6gen();

%----------------
% Synchronization 
%----------------
tm6_ov  = reshape(repmat(tm6, ov, 1), 1, []);
xc      = filter(tm6_ov(end:-1:1), 1, [xcap zeros(1, length(tm6_ov))]);
[mv mi] = max(abs(xc));
dly     = mi - length(tm6_ov);
xcap    = xcap(1+dly:end);
xcap    = xcap(1:length(tm6_ov));

%----------------
% Decimation and clock phase recovery
%----------------
alpha = 0.7;

min_ted = Inf;
for i = 0:ov-1,
  xcap_dec = xcap(1+i:ov:end);
  l0       = min([length(xcap_dec) length(tm6)]);
  ted      = mean((1 - alpha)*xcap_dec(2:l0).*tm6(1:l0-1) - ...
                  alpha*xcap_dec(1:l0-1).*tm6(2:l0));

  if abs(ted) < min_ted
    min_ted = abs(ted);
    dly     = i;
  end
end
xcap_dec = xcap(1+dly:ov:end);

%----------------
% Normalize
%----------------
xcap_dec = xcap_dec/max(abs(xcap_dec));


%----------------
% Volterra estimation & analysis
%----------------
[HD2 HD3 RPD] = volest(tm6, xcap_dec, 12, 3);



%------------------------------------------------------------------------------
% Generate the reference test signal pattern for test mode 6
%------------------------------------------------------------------------------
function tm6 = tm6gen()

Ns = 2^11 - 1;

scr1 = lfsrgen([1 9 11], Ns, '7FF').';

x1 = [circshift(scr1, 0), ...
      mod((circshift(scr1, 1) + circshift(scr1, 4)), 2), ...
      mod((circshift(scr1, 2) + circshift(scr1, 9)), 2), ...
      mod((circshift(scr1, 0) + circshift(scr1, 10)), 2)];

y1 = x1(:,1) + 2*x1(:,2) + 4*x1(:,3) + 8*x1(:,4);


scr2 = lfsrgen([1 7 9 10 11], Ns, '7FF').';

x2 = [circshift(scr2, 0), ...
      mod((circshift(scr2, 1) + circshift(scr2, 4)), 2), ...
      mod((circshift(scr2, 2) + circshift(scr2, 9)), 2), ...
      mod((circshift(scr2, 0) + circshift(scr2, 10)), 2)];

y2 = x2(:,1) + 2*x2(:,2) + 4*x2(:,3) + 8*x2(:,4);

tm6 = ((2*(16*y1 + y2) - 255)/256).';


%------------------------------------------------------------------------------
% Volterra's estimation and analysis
%------------------------------------------------------------------------------
function [HD2 HD3 RPD] = volest(x, d, n, dly)

% Init
R  = 0;
rD = 0;

% Addition autocorrelation and cross-correlation
for k = n:length(x),
  % Volterra products
  xi = [1 ...
    x(k:-1:k-n+1) ...
    x(k:-1:k-n+1)  .*x(k:-1:k-n+1)   ...
    x(k:-1:k-n+1+1).*x(k-1:-1:k-n+1) ...
    x(k:-1:k-n+1+2).*x(k-2:-1:k-n+1) ...
    x(k:-1:k-n+1)  .*x(k:-1:k-n+1)    .*x(k:-1:k-n+1)   ...
    x(k:-1:k-n+1+1).*x(k:-1:k-n+1+1)  .*x(k-1:-1:k-n+1) ...
    x(k:-1:k-n+1+2).*x(k:-1:k-n+1+2)  .*x(k-2:-1:k-n+1) ...
    x(k:-1:k-n+1+1).*x(k-1:-1:k-n+1)  .*x(k-1:-1:k-n+1) ...
    x(k:-1:k-n+1+2).*x(k-1:-1:k-n+1+1).*x(k-2:-1:k-n+1) ...
    x(k:-1:k-n+1+2).*x(k-2:-1:k-n+1)  .*x(k-2:-1:k-n+1)];

  % Autocorrelation matrix
  R = R + xi.'*xi;

  % Cross-correlation vector
  rD = rD + d(k-dly)*xi.';
end

% Wiener's MMSE solution
hw = (R\rD).';

% Separate the Volterra kernels per channel
lw = [1 ...
      n ...
      n (n-1) (n-2) ...
      n (n-1) (n-2) (n-1) (n-2) (n-2)];

ofst = 0;
for i = 1:11,
  h{i} = hw(ofst+1:ofst+lw(i));
  ofst = ofst + lw(i);
end

% Calculate harmonic distortion
HD2 = -10*log10(1/3*axc(h{2})/(1/5*axc(h{3})  + 1/9*axc(h{4})   + 1/9*axc(h{5})));
HD3 = -10*log10(1/3*axc(h{2})/(1/7*axc(h{6})  + 1/15*axc(h{7})  + 1/15*axc(h{8}) + ...
                               1/15*axc(h{9}) + 1/27*axc(h{10}) + 1/15*axc(h{11})));

% Calculate residual peak distortion
z = h{1} + ...
    filter(h{2},  1, x(3:end)) + ...
    filter(h{3},  1, x(3:end).*x(3:end))   + ...
    filter(h{4},  1, x(3:end).*x(2:end-1)) + ...
    filter(h{5},  1, x(3:end).*x(1:end-2)) + ...
    filter(h{6},  1, x(3:end).*x(3:end)  .*x(3:end))   + ...
    filter(h{7},  1, x(3:end).*x(3:end)  .*x(2:end-1)) + ...
    filter(h{8},  1, x(3:end).*x(3:end)  .*x(1:end-2)) + ...
    filter(h{9},  1, x(3:end).*x(2:end-1).*x(2:end-1)) + ...
    filter(h{10}, 1, x(3:end).*x(2:end-1).*x(1:end-2)) + ...
    filter(h{11}, 1, x(3:end).*x(1:end-2).*x(1:end-2));

z = z(1+dly-2+n:end);
d = d(1+n:end);
l = min([length(z) length(d)]);
e = z(1:l) - d(1:l);

RPD = 20*log10(max(abs(e))/abs(sum(h{2})));


%------------------------------------------------------------------------------
% Auto-correlation, main term
%------------------------------------------------------------------------------
function xc = axc(x)
if length(x) > 0
  rhc = xcorr(x);
  xc = rhc(length(x));
else
  xc = 0;
end


%------------------------------------------------------------------------------
% LFSR
%------------------------------------------------------------------------------
function out = lfsrgen(poly, len, seed)

r = double(dec2bin(hex2dec(seed))) - 48;
r = [zeros(1,poly(end)-length(r)) r];

for i = 1:len,
  out(i) = r(1);
  fb = 0;
  for j = 2:length(poly),
    fb = mod(fb + r(poly(j)), 2);
  end
  r = [fb r(1:poly(end))];
end
