Back Substitution - Revision history

Ranjelin at 02:35, 30 December 2024

2024-12-30T02:35:20Z

←Older revision		Revision as of 02:35, 30 December 2024
Line 42:		Line 42:

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

Ranjelin at 02:25, 30 December 2024

2024-12-30T02:25:33Z

←Older revision		Revision as of 02:25, 30 December 2024
Line 29:		Line 29:
	thetai = atan2(-imag(fn(i,:)), real(fn(i,:)));		thetai = atan2(-imag(fn(i,:)), real(fn(i,:)));
	% signum(i,:) = mypsi(cos(itheta1)).mypsi(cos(thetai));		% signum(i,:) = mypsi(cos(itheta1)).mypsi(cos(thetai));
-	%%% "Perhaps better for computation..."	+	%%% "Perhaps better for computation..." %%%
	varpi = abs(alias(thetai + pi - itheta1, twopi)) < abs(alias(thetai - itheta1, twopi));		varpi = abs(alias(thetai + pi - itheta1, twopi)) < abs(alias(thetai - itheta1, twopi));
	signum(i,varpi) = -one;		signum(i,varpi) = -one;

Ranjelin: New page:

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

2024-12-30T02:24:42Z

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

New page

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