lfsrc.m

来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 122 行

% lfsrc.m

function [e0_1,alpha_1,wg_1,eps_1,ece_1]=lfsrc(tc,tp,te,ta,ee)
% Convert LF timing parameters into direct synthesizing parameters.

%	tc=58.58;tp=41.41;te=55.55;ta=0.404;ee=60;

	count = 1;ue = 0.0;area_ret = 0.0;area_open = 0.0;
	dist_top_alpha_1 = 0.001; add_top = 0.1;tot_area = 0.0;
	bottom = 0.0; top = 0.1;

	ite = round(te);
	itc = round(tc);
	wg_1 = pi/tp;    % wg_1=2*pi*Fg, where Fg = 1/(2*tp) 

% if no return phase 
	if (ta == 0.0)
		index=3;
	else
%  Newton-Raphson method to solve for eps_1ilon 
	eps_1x = 1.0/ta;
	eps_1=eps_1x;
	ce = tc - te;
	index=1;
	end

while count <= 1000,

	if index==1,
	ece_1 = exp(-(eps_1)*ce);
	feps_1 = ((eps_1)*ta) - (1.0-(ece_1));
	dfeps_1 = ta - ((ece_1)*ce);
	eps_1 = eps_1 - (feps_1/dfeps_1);
	if (abs(eps_1x-(eps_1)) <= 1.0e-7)
		ece_1 = exp(-(eps_1)*ce);
		index=2;
	else
	eps_1x = eps_1;
	index=1;
	end
	end

% approximation to the area under the return phase 

	if index==2,
	area_ret=0.0;
	for i=ite+1:itc,
		area_ret=area_ret-1.0*(exp(-eps_1*(i-ite))-ece_1)/(eps_1*ta); 
		ue=area_ret;
	end
%
% in the next formula area under return phase should be positive since 
% the negative sign has been taken care of.
 
	ue = -1.0*ue;
	index=3;
	end

	if index==3,
	st = sin(wg_1*te);
	wu = (wg_1^2.0)*ue - (wg_1*cos(wg_1*te)/(st*tp));

% Wegstein's method to solve for alpha_1 
	a0 = 0.5*(wg_1);
	ea = exp(-a0*te);
	f0 = -((wg_1)*ea/st) -(((a0*a0*ue)+wu)*tp); 
	a1 = f0;
	index=4;
	end

	if index==4,
	ea = exp(-a1*te);
	f1 = -((wg_1)*ea/st) -(((a1*a1*ue)+wu)*tp); 
	alpha_1 = a1 + ((a1-a0)/(((a0-f0)/(a1-f1))-1.0));
	if (abs(a1-(alpha_1)) <= 1.0e-7)
		index=5;
	else
	f0 = f1;
	a0 = a1;
	a1 = (alpha_1);
	index=4;
	end
	end

% compute e0_1n (norm. e0_1) 
	if index==5,
	e0_1n = -1.0/(exp(alpha_1*te)*st);

% calculate the actual area under the open phase
	area_open=0.0; 
	for i=1:ite,
		area_open=area_open+e0_1n*exp(alpha_1*i)*sin(wg_1*i);% opening phase 
	end
% calculate total area under the pulse 
	tot_area = area_open + area_ret;
	if (abs(tot_area) < 0.01)
% compute e0_1
		e0_1 = -ee/(exp(alpha_1*te)*st);
		return;
	else 
		if (tot_area > 0.0) 
			if (abs(alpha_1-top) < dist_top_alpha_1)
				top = top+add_top;
			end
			bottom = alpha_1;
		else 
			top = alpha_1;
		end
		alpha_1 = (top+bottom)/2.0;
		count=count+1;
		if(count > 1000) 
% fputs("Warning: alpha_1 did not converge for the current frame\n", stderr);
% compute e0_1 
			e0_1 = -ee/(exp(alpha_1*te)*st);
			return;
		else 
			index=5;
		end
	end
	end
end