%FUNCTION: generate the glottal pulse using the LF model.
%          gpulse=lfpuls(Tp,Te,Ta,Tc,Ee,lens,noise) returns the extimated
%          glottal wavform.
%
% INPUT:   Tp,Te,Ta,Tc = four timing parameters for the LF model
%          Ee = max negative differential waveform value
%          lens = the length of the pitch period
%          noise = superpose the white noise if 'noise~=0'
%
% OUTPUT:  
%       gpulse = an estimated differential glottal waveform
%
% SEE ALSO: lfmodel.

function gpulse=lfpuls(Tp,Te,Ta,Tc,Ee,lens,noise);

if nargin==6
   noise=0;
end

% issue warning message
if Te>Tc
   disp('Error in LF timing parameters; Te should be less than Tc');
   disp(['Te=' num2str(Te) ' Tc=' num2str(Tc) ]);
   Tc=1;
   return;
end

if Tp>Tc
   disp('Error in LF timing parameters; Tp should be less than Tc');
   disp(['Tp=' num2str(Tp) ' Tc=' numstr(Tc) ]);
   return;
end

if Tp>Te
   disp('Error in LF timing parameters; Tp should be less then Te');
   disp(['Tp=' num2str(Tp) ' Te=' num2str(Te)]);
   return;
end  

 
% 1. estimate the LF computing parameters
  T0=lens;
  Tp=Tp*T0;
  Te=floor(Te*T0);
  Ta=Ta*T0;
  Tc=floor(Tc*T0);

 [alpha,eps1,ece,wg,e0,OK]=lfsrc(Tp,Te,Ta,Tc,Ee);
 if OK==0
     return;
 end

% 2. generate the LF pulse

 gp=zeros(1,lens);
 for i=1:Te
     gp(i)=e0*exp(alpha*i)*sin(wg*i);
 end

 for i=Te:Tc
     tmp_gp=-Ee*exp((-eps1*(i-Te)-ece))/(eps1*Ta);
     if tmp_gp<0.0
        gp(i)=tmp_gp;
     else
        gp(i)=0;
     end
 end 

 gpulse=gp;
 gpulse=[gp(Te:T0) gp(1:Te-1)];

if noise==1
   
   gpulse0=gpulse;
 % 3. produce high-energy pulse at glottal closure (follow Dr.Hu's approach)

  % create a pulse swing
  gpulse(1)=gpulse(lens);
  gpulse(3)=0.6*gpulse(2)+0.4*gpulse(1);
  gpulse(6)=0.5*gpulse(6)+0.5*gpulse(3);
  gpulse=gpulse-mean(gpulse);

  % Remove the spectral tilt of the excitation pulse by inverse filtering
  ss=gpulse(2:lens);
  ss=ss(:);
  energy=ss'*ss;
  rc1=ss(1:lens-2)'*ss(2:lens-1)/energy; % first order reflection coefficient
  gpulse=filter([1 -.75*rc1],1,gpulse);  % first-order inverse filtering
  gpulse(1)=gpulse(lens);

  gpulse=0.3*gpulse+0.7*gpulse0;

end

if noise>0

% 4. Add modulated noise.

   W=cshift(gauss(-1:2/lens:1,0,.25),floor(lens/2))+.5;
   nn=(rand(1,lens+1)-.5).*W;
   nn=nn(1:lens);
   amp=sqrt(0.0032*sum(gpulse.^2)/(nn*nn'));
   gpulse=gpulse+nn*amp;

end