thresh_gauss_irreg.m

来自「software to find the threshold of an LDP」· M 代码 · 共 44 行

function y = thresh_gauss_irreg(rho,lambda,K,delta,lpad,iter,mu_start,mu_inc,mu_term)

x_in = (-1*K*delta):delta:(K*delta);
f_in = exp(-1/(4*mu_start)*(x_in-mu_start).^2);
f_in = f_in / sum(f_in);

mu_start

[y,init_df] = feed_irreg(rho,lambda,f_in,K,delta,lpad,iter);

df = init_df;
mu = mu_start;

while (df == init_df)
	last_mu = mu;

	mu = mu + ((-1)^(1-df))*mu_inc; % adds mu_inc if decoding failure, otherwise subtracts
mu
	f_in = exp(-1/(4*mu)*(x_in-mu).^2);
	f_in = f_in / sum(f_in);
	[y,df] = feed_irreg(rho,lambda,f_in,K,delta,lpad,iter);
end

% now, the threshold is between last_mu and mu
if (last_mu > mu)
	t_limits = [mu last_mu];
else
	t_limits = [last_mu mu];
end

while ((t_limits(2)-t_limits(1)) > mu_term)
	foo = (t_limits(2)-t_limits(1))/2 + t_limits(1);
foo
	f_in = exp(-1/(4*foo)*(x_in-foo).^2);
	f_in = f_in / sum(f_in);
	[y,df] = feed_irreg(rho,lambda,f_in,K,delta,lpad,iter);
	if (df == 1)
		t_limits(1) = foo;
	else
		t_limits(2) = foo;
	end
end

y = t_limits(2);