function markerarray=create_monopuls(pitcharray,fs,markerfile,n);
% DOEL: berekent uit de gegeven pitchwaarden een rij van markers 
%       die de pitch perioden aanduiden.
% 
% INPUT: pitcharray     array met pitchwaarden van geluid G
%        fs             sampling frequentie van G (in Hz)
%        markerfile     naam van de .txt-file waarin de outputs worden geschreven
%        n              aantal samples tussen de pitchwaarden in pitcharray plus 1,
%                        ( n/fs  geeft dus de tijd tussen 2 pitchwaarden in)
%
% OUTPUT: markerarray   rij van markers die de pitch perioden aanduiden.
%
% OPMERKING : als markerfile = 'nowrite' dan wordt er geen file aangemaakt.
%
% VOORBEELD: markerarray=create_monopuls(pitcharray,16000,'nowrite',320)

a=size(pitcharray); numberofsamples=n*a(1);
markerarray=1;
markerarray(2,1)=1;
current=1;
i=2;
while current<numberofsamples
   current=nextmark3(pitcharray,current,fs,n);
   if current<=numberofsamples
      markerarray(i)=round(current);
      i=i+1;
   end;
end;
if strcmp(markerfile,'nowrite')==0;
   fid=fopen(markerfile,'w+');
   fprintf(fid,'%40.20f',markerarray);
   fclose(fid);
end;

% ---------------------------------------------------------
% Hulpfunctie
% ---------------------------------------------------------

function K=nextmark3(pitcharray,previousmark,fs,n);
% Deze functie berekent, gegeven een pitchwaarde om de n samples
% en de (niet-gehele) index van de vorige marker 'previousmark',
% de (niet-gehele) index van de volgende marker.

% Initialisaties
i=0; % Dit is een teller, zie verder.
xold=0;
x=0; % 
boundary1=previousmark; 

remain=rem(floor(previousmark)-1,n);
boundary2=floor(previousmark)-remain+n;

a=size(pitcharray); numberofsamples=n*a(1); 

pitch1=fs/pitcharray((boundary2-1)/n); % De functie 'pitch' leest uit het bestand 'pitchbestand' de pitchwaarde in op een bepaald sample.
if ((boundary2-1)/n+1)<=length(pitcharray)
   pitch2=fs/pitcharray((boundary2-1)/n+1);
else
   pitch2=pitch1;
end;
m=(pitch2-pitch1)/n; % Hellingscoëfficiënt.
c=pitch1-m*(boundary2-n); % Constante term van de rechte die de pitch beschrijft
% tussen twee samples in (lineaire interpolatie).

while x<1;
   xold=x;
   if (c+m*boundary1)==0
      error2='c+ma=0'
      error2
      i
   end;
   if m==0;
      x=x+(boundary2-boundary1)/c;
   else;
      x=x+log((c+m*boundary2)/(c+m*boundary1))/m; % Berekening deelintegraal.
   end;
   i=i+1; % Teller updaten, de teller telt hoeveel deelintegralen
   % al berekend zijn.
   % Updaten van de variabelen die de deelintegralen bepalen:
   if x<1;
      boundary1=boundary2;
      pitch1=pitch2;
      boundary2=boundary2+n;
      if boundary2>numberofsamples;
         K=numberofsamples+100;
         return;
      else;
      pitch2=fs/pitcharray((boundary2-1)/n+1);
      m=(pitch2-pitch1)/n;
      c=pitch1-m*boundary1;
 end;
end;
if m==0;
   K=boundary1+c*(1-xold);
else;
   K=((c+m*boundary1)*exp(m*(1-xold))-c)/m; % Berekening waar x juist
   % de waarde 1 overschrijdt.
end;
end;