THD from Mapping Coefficients

From Wikimization

%THD from mapping coefficients.  Submeasurable Op Amp Distortion, section 3
function [thd dc] = thdxi(xi, E, precision);
   mp.Digits(precision);  %Advanpix MCT
   Nh = numel(xi);
   two = mp('2');
   harmonic = zeros(Nh+1,1,'mp');
   for n = 1:Nh
      tn = xi(n)*E^n/two^n;
      for ell = 0:n
         tell = tn*nchoosek(mp(n), mp(ell));
         idx = abs(n - 2*ell) + 1;
         harmonic(idx) = harmonic(idx) + tell;
      end
   end
   thd = sum(harmonic(3:end).^2)/harmonic(2)^2;
   dc  = harmonic(1);
end