lfsrc.m

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

% This program was used to transform the LF timing parameters: tp,te,ta,tc
% to its direct paramaters: alpha, ece, eps1, wg, e0
% The whole program is translated from C program written by Ajit Lalwani
function [alpha,eps1,ece,wg,e0,OK] = lfsrc(tp,te,ta,tc,ee)

ite=floor(te);
itc=floor(tc);
OK=1;

% Given initial condition

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

% define wg

wg=pi/tp;

% 
if ta ~= 0.0
	ra=ta/(tc-te);
	if ra < 0.1
		ak = 2.0;
	elseif ( ra >=0.1 & ra <= 0.5 )
		ak = 2.0 - 2.34*(ra^2) + 1.34*(ra^4);
	else
		ak = 2.16 - 1.32*ra + 0.64*((ra-0.5)^2);
	end
	epsx = 1.0/ta;
	eps1 = epsx;
	ce = tc - te;
	ece = exp(-eps1*ce);
	feps = eps1*ta - (1-ece);
	dfeps = ta - (ece*ce);
	eps1 = eps1 - (feps/dfeps);
	while ((abs(epsx-eps1)) > 1.0e-7)
		epsx=eps1;
		ece = exp(-eps1*ce);		% level 1
		feps = eps1*ta - (1-ece);
		dfeps = ta - (ece*ce);
		eps1 = eps1 - (feps/dfeps);
	end

	% calculate the actual area under the return phase
	area_ret=0.0;
	for i = ite+1 : itc
		area_ret = area_ret - 1.0*(exp(-eps1*(i-ite))-ece)/(eps1*ta);
	end
	ue = -1.0*area_ret;
end

% level 3 

st = sin(pi*te/tp);
wu = (wg^2)*ue - wg*cos(pi*te/tp)/(st*tp);

% Wegstein's method to solve for alpha

a0 = 0.5*wg;
ea = exp(-a0*te);
f0 = -wg*ea/st - (a0*a0*ue+wu)*tp;
a1 = f0;

% level 4

ea = exp(-a1*te);
f1 = -wg*ea/st - (a1*a1*ue+wu)*tp;
alpha = a1 + ((a1-a0)/(((a0-f0)/(a1-f1))-1.0));
while (abs(a1-alpha) > 1e-7)
	f0 = f1;
	a0 = a1;
	a1 = alpha;
	ea = exp(-a1*te);
	f1 = -wg*ea/st - (a1*a1*ue+wu)*tp;
	alpha = a1 + ((a1-a0)/(((a0-f0)/(a1-f1))-1.0));
end

% level 5
% compute e0n , normal e0

e0n = -1.0/(exp(alpha*te)*st);

% calculate the actual area under the open phase
area_open = 0.0;
for i = 1 : ite
	area_open = area_open + e0n*exp(alpha*i)*sin(wg*i);
end

% calculate total area under the pulse
tot_area = area_open + area_ret;
while abs(tot_area) >= 0.01
	if tot_area > 0.0
		if abs(alpha-top) < dist_top_alpha
			top = top + add_top;
		end
		bottom = alpha;
	else
		top = alpha;
	end
	alpha = (top+bottom)/2;
	e0n = -1.0/(exp(alpha*te)*st);
	area_open = 0.0;
	for i = 1 : ite
		area_open = area_open + e0n*exp(alpha*i)*sin(wg*i);
	end
	tot_area = area_open + area_ret;
	count = count + 1;
	if count > 1000
		tot_area = 0;
		fprintf('Warning: alpha did not converge \n');
                OK=0;
	end
end

% level 6
% compute e0

e0 = -ee/(exp(alpha*te)*st);