Back Substitution

From Wikimization

%mapping coefficients by back substitution.
%This routine is for thd+n measurement only.
%fn is function of time in optional second dimension. 
function [zeta signum flip] = backsubz2(E, fn, Nh, precision)  %fn is complex 2*fft/N
   if nargin < 4, precision = 34; end
   mp.Digits(precision);  %requires Advanpix Multiprecision Toolbox
   if size(fn,1) ~= Nh, disp(' ');disp('backsubz2(): fn length error');disp(' '); return; end 
   zeta   = zeros(Nh,size(fn,2),'mp');      
   signum = ones (Nh,size(fn,2),'mp'); 
   flip   = 0;
   zero   = mp('0');
   one    = mp('1');
   two    = mp('2');
   pi     = mp('pi');
   twopi  = mp('2*pi');
   theta1 = atan2(-imag(fn(1,:)), real(fn(1,:)));
   for i=Nh:-1:1  
      hi1 = zero;
      if mod(i,2)  %odd
         for n=(i+3)/2:ceil(Nh/2)
            hi1 = hi1 + nchoosek(two*n-one, n-(i+one)/two)*zeta(2*n-1,:);
         end
      else         %even
         for n=(i+2)/2:floor(Nh/2)
            hi1 = hi1 + nchoosek(two*n,     n-i/two)      *zeta(2*n,:);
         end   
      end
      thetai = atan2(-imag(fn(i,:)), real(fn(i,:)));
%     signum(i,:) = mypsi(cos(i*theta1)).*mypsi(cos(thetai));
      %%% "Perhaps better for computation..." %%%
      varpi = abs(alias(thetai + pi - i*theta1, twopi)) < abs(alias(thetai - i*theta1, twopi));
      signum(i,varpi) = -one;
      zeta(i,:) = signum(i,:).*abs(fn(i,:)) - hi1;
   end
   zeta = two.^((1:Nh)'-1).*zeta./abs(E.^(1:Nh)');
   t = find(zeta(1,:) < 0);  %back substitution assumes positive zeta(1,:)
   if any(t), flip = 1; end
   zeta  (1:2:Nh,t) = -zeta  (1:2:Nh,t);
   signum(1:2:Nh,t) = -signum(1:2:Nh,t);
end

%nonplatonic alias frequency function
function fa = alias(f,Fs)  % Fs is sample rate [Hz] or 2*pi. 
   fa = f - Fs*round(f/Fs);
end