Ranjelin at 05:33, 17 January 2025

2025-01-17T05:33:54Z

←Older revision		Revision as of 05:33, 17 January 2025
Line 14:		Line 14:
	end		end
	mp.Digits(precision);		mp.Digits(precision);
-	if numel(theta) == 1	+	if numel(theta) == 1 %preferably nonlinearity output fundamental phase
	theta = (1:Nh)'*theta; %linear phase		theta = (1:Nh)'*theta; %linear phase
	elseif isrow(theta)		elseif isrow(theta)

Ranjelin at 05:19, 11 January 2025

2025-01-11T05:19:36Z

←Older revision		Revision as of 05:19, 11 January 2025
Line 91:		Line 91:
	lp_failure = true;		lp_failure = true;
	end		end
		+	end
		+	end
		+
		+	%adjust mapping coefficients for unity gain
		+	function xi = ug(xi, Ep, precision)
		+	version -blas; version -lapack; %load AMD AOCL libraries
		+	if nargin < 3
		+	precision = 34;
		+	end
		+	mp.Digits(precision);
		+	zero = mp('0');
		+	one = mp('1');
		+	two = mp('2');
		+	Nh = numel(xi);
		+	t = zero;
		+	for n=2:ceil(Nh/2)
		+	t = t + nchoosek(twon-one,n-one)(Ep/two)^(twon-two)xi(2*n-1);
		+	end
		+	xi(1) = one - t; %adjust for unity gain.
		+	if xi(1) < 0
		+	xi(3:2:Nh) = -xi(3:2:Nh); %negativity check
		+	xi(1) = one + t;
	end		end
	end		end
	</pre>		</pre>

Ranjelin: New page:

 %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] =...

2025-01-11T05:17:51Z

New page: <pre> %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] =...

New page

<pre>
%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
</pre>

Frequency Domain Inversion - Revision history

Ranjelin at 05:33, 17 January 2025

Ranjelin at 05:19, 11 January 2025

Ranjelin: New page:
%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] =...