Frequency Domain Inversion

From Wikimization

%computation of mapping coefficients by frequency domain inversion. August 9 2023
%this subroutine is used for thd+n measurements only; single test tone.
function [zeta lp_failure] = backsubzfreq(Ep, E, fn, Nh, theta, precision) %fn is complex conjugate data 2*conj(fft)/N.
   version -blas; version -lapack;  %load AMD AOCL libraries
   maximallyflat        = ~true;                                           %theta is vector or scalar [rads] modifying twiddle W
   closed               =  true;  %closed form solution. nonclosed is for backsubz3Test (thd); linprog prone to failure for small amplitude test tone or high Nh
      matrix_Inversion  = ~true;  %choose closed form method: matrix inversion or analytical inverse
   unitygainconstraint  = ~true;  %amplitude Ep used for unity gain xi determination. E used only to make twiddle matrix.
   lorenzo              = ~true;  %Lorenzo David closed form twiddle matrix.
   if size(fn,1) ~= Nh, disp(' ');disp('backsubzfreq(): fn length error');disp(' '); return; end 
   if nargin < 6
      precision = 34;
   end
   mp.Digits(precision);
   if numel(theta) == 1
      theta = (1:Nh)'*theta;  %linear phase
   elseif isrow(theta)
      theta = theta';
   end
   zero = mp('0');
   one  = mp('1');
   two  = mp('2');
   W = zeros(Nh+1,Nh,'mp');
   C = zeros(1,Nh,'mp');
   %%% make real twiddle matrix
   if lorenzo
      for i=1:Nh
         for j=i:Nh
            if ~mod(i+j,2)
               W(i,j) = two*(E/two)^j*nchoosek(mp(j),(j-i)/two); %second binomial coefficient (i+j)/2 is equivalent to (-i+j)/2 because of numerator j
            end 
         end
      end
   else
      for i=1:Nh  
         k = 0;
         idx = i:2:Nh;
         for j=idx  %binomial coefs
            C(j) = nchoosek(mp(j),mp(k));
            k = k + 1;
         end
         W(i,idx) = two*(E/two).^idx.*C(idx);  %E corresponds to measured signal amplitude in frequency domain fn
      end
   end
   if maximallyflat
      for n=2:ceil(Nh/2) 
         W(Nh+1,2*n-1) = nchoosek(two*n-one,n-one)/two^(2*n-2);
      end
      theta = [theta; zero];
      fn    = [fn;    zero];
   else
      W(Nh+1,:) = [];
   end   
   %%% mapping coef solution by closed form projection OR by linear algebraic inversion (linprog, simplex)
   lp_failure = false;
   if closed
      if matrix_Inversion
         zeta = real((exp(1j*theta).*W)\fn);
      else
         zeta = zeros(Nh,1,'mp');  %from "monotone frequency domain inversion" in chapter 8
         for j=1:Nh
          % for k=j:floor((Nh+j)/two)
          %    zeta(j) = zeta(j) + real((-one)^(k-j)*(two*k-j)*nchoosek(k-one,mp(k-j))*exp(-1j*theta(2*k-j))*fn(2*k-j));
          % end
            for k=0:floor((Nh-j)/two)  %Nov.24 2024 changing real to abs fails backsubz3Test.m
             % zeta(j) = zeta(j) +  abs((two*k+j)*binomial(mp(-j),mp(k),precision)*exp(-1j*theta(2*k+j))*fn(2*k+j));
               zeta(j) = zeta(j) + real((two*k+j)*binomial(mp(-j),mp(k),precision)*exp(-1j*theta(2*k+j))*fn(2*k+j));
            end
            zeta(j) = two^(j-1)*zeta(j)/(j*E^j);  %two^(j-1) because DFT already multiplied by 2
         end
      end
      if unitygainconstraint
         zeta = ug(zeta, Ep, precision);       %unity gain and negativity check
      end
   else
      Aeq = double([real(exp(1j*theta).*W) 
                    imag(exp(1j*theta).*W)]);
      beq = double([real(fn) 
                    imag(fn)]);
      if unitygainconstraint
         D = zeros(1,Nh);
         for n=1:ceil(Nh/2)
            D(1,2*n-1) = nchoosek(two*n-one,n-one)*(Ep/two)^(two*n-two);
         end
         Aeq = [Aeq; D];  beq = [beq; 1];  %unity gain constraint.
      end
      [zeta,~,flag] = linprog(ones(Nh,1), [], [], Aeq, beq, [0; -inf(Nh-1,1)]); %prone to failure for small amplitude test tone or high Nh
      if flag < 0
         zeta = zeros(Nh,1,'mp');
         lp_failure = true;
      end
   end
end