unit MMStretch;

(*-------------------------------------------------------------------
 * Time domain harmonic scaling by
 * Pointer Inteval Controled OverLap and ADD (PICOLA) Method
 *		C version by IKEDA Mikio
 *		original argolithm is developed by MORITA Naotaka
 *		about detail, see original paper.
 *-------------------------------------------------------------------
 * Usage
 *  PICOLA <source signal> <companded (destination) signal>
 *	   <compansion ratio>
 *         <window length> <pitch minimum> <pitch maximum>
 * Last three arguments can be abbriviated.
 *------------------------------------------------------------------*)

 // Does not work (horrible quality) leave it for now....

interface

uses
    SysUtils,
    Windows,
    Classes,
    MMSystem,
    MMRegs,
    MMObj,
    MMDSPObj,
    MMPCMSup,
    MMUtils;


type
   EMMTimeStretchError = class(Exception);

   {-- TMMTimeStretch ---------------------------------------------------------}
   TMMTimeStretch = class(TMMDSPComponent)
   private
      FEnabled       : Boolean;
      FOpen          : Boolean;
      FFirstRead     : Boolean;
//      FPitch         : Float;
      FWaveHdr       : TMMWaveHdr;
      FRealBufSize   : Longint;
      FBytesRead     : Longint;
      FMoreBuffers   : Boolean;
      FWriteBuffer   : PChar;
      FBytesWritten  : Longint;
      FDone          : Boolean;

      procedure SetEnabled(aValue: Boolean);
  //    procedure SetPitch(aValue: Float);
      function  ReadData(Buffer: PChar; dwLength: Longint; var MoreData: Boolean): Longint;
      function  WriteData(Buffer: PChar; dwLength: Longint): Longint;
      procedure ReadFromInput(lpwh: PWaveHdr; var MoreBuffers: Boolean);

   protected
      procedure SetPWaveFormat(aValue: PWaveFormatEx); override;
      procedure Opened; override;
      procedure Closed; override;
      procedure Started; override;
      procedure Reseting; override;
      procedure BufferReady(lpwh: PWaveHdr); override;
      procedure BufferLoad(lpwh: PWaveHdr; var MoreBuffers: Boolean); override;

   public
      constructor Create(aOwner: TComponent); override;
      destructor Destroy; override;

      procedure Open;
      procedure Close;
      procedure Reset;

   published
      property Input;
      property Output;
      property Enabled: Boolean read FEnabled write SetEnabled default True;
    //  property Pitch: Float read FPitch write SetPitch;
   end;

procedure StretchFile(SrcFile,DstFile: String);

implementation

{---- find maximum covariance  = pitch ----------------------------------------}
function covpitch(pitmin, pitmax, length: integer; _is: array of Smallint): integeR;
var
  i,j,pitch: integer;
  covst, covs0t, covmax, s: Float;

begin
   covmax := 0.0;

   pitch := pitmin;
   for i := pitmin to pitmax do
   begin
      covst := 0.0;
      covs0t := 0.0;
      for j := 0 to length-1 do
      begin
	 s := _is[i+j];
	 covs0t := covs0t + s * s;
	 covst  := covst + _is[j] * s;
      end;
      covst := covst / sqrt(covs0t);
      if (covst >= covmax) then
      begin
         covmax := covst;
         pitch := i;
      end;
   end;
   Result := pitch;
end;

{------------ PICOLA OverLap and Add (picOLA) stage ---------------------------}
procedure ola(pitch: integer; is1, is2: PSmallint);
var
   i: integer;
   s, w: Float;
begin
   for i := 0 to pitch-1 do
   begin
      w := i / (pitch - 1);
      s :=  is1^ * (1.0 - w) + is2^ * w;
      inc(is1);
      is2^ := Trunc(s);
      inc(is2);
   end;
end;

{------------------------------------------------------------------------------}
function amdfpitch(pitmin, pitmax, length: integer; _is: array of Smallint): integer;
var
   i, j, diff, acc, accmin, pitch: integer;

begin
   pitch := pitmin;
   accmin := 0;
   for j := 0 to length-1 do
   begin
      diff := _is[j+pitmin] - _is[j];
      if (diff > 0) then
          accmin := accmin + diff
      else
          accmin := accmin - diff;
   end;

   for i := pitmin+1 to pitmax do
   begin
      acc := 0;
      for j := 0 to length-1 do
      begin
         diff := _is[i+j] - _is[j];
	 if (diff > 0) then
             acc := acc + diff
         else
	     acc := acc - diff;
      end;

      if (acc < accmin) then
      begin
	 accmin := acc;
	 pitch := i;
      end;
   end;
   Result := pitch;
end;

var
   _is : array[0..4096] of Smallint; // signal buffer
   rate: Float = 1.1;	// compansion rate
			//case of less than 1.0 compression,
			//case of greater than 1.0 expansion

   rcomp: Float;	// internal compansion ratio
   sl: Float;
   err: float = 0.0;	// compansion rate error estimate
   acclen: Float = 0.0;

   pitmin: integer = 32;  // minimal pitch period //
   pitmax: integer = 1024; // maximal pitch period //
   pitch : integer; 	  // detected pitch period */
   length: integer = 1024;

   total: integer;
   nread: integer;	// number of read samples (from file) */
   wantread: integer;	// desired number of read samples */
   lcp: integer;	// number of copy samples */
   point: integer;	// PICOLA's pointer */
//   i: integer;		// loop counter */
   lproc: integer = 0;	// processed speech samples */

   Src,Dst: THandle;

procedure StretchFile(SrcFile,DstFile: String);
var
   i: integer;
begin
   // length := atoi(argv[4]); option
   // pitmin := atoi(argv[5]); option
   // pitmax := atoi(argv[6]); option

   //-------------- error check and initialize ---------------------

   {
    if (rate <= 0.0 || rate == 1.0)
    begin
       printf("illeagal compansion rate !!\n");
	exit(0);
    end;

    if (pitmin < 16)
    begin
        printf("pitch detection minimum threshold modified !!\n");
	pitmin = 16;
    end;

    if (pitmax > 256)
    begin
        printf("pitch detection maximum threshold modified !!\n");
	pitmax = 256;
    end;

    if (length <= 64 || length + pitmax >= 1024)
    begin
	printf("frame length out of range !!\n");
	exit(0);
    end;
}
    total := length + pitmax;
    if (rate >= 1.0) then
    begin
       // TODO:rate darf nicht 1.0 sein, also bei rate := 1.0 skippen
       rcomp := 1.0  / (rate - 1.0);
    end
    else if (rate > 0) then
    begin
       rcomp := rate / (1.0 - rate);
    end
    else
    begin
       //	fprintf(stderr, "Error from %s: illeagal compansion rate!\n", argv[0]);
       //	exit(0);
    end;

    Src := FileOpen(SrcFile,fmOpenRead);
    Dst := FileCreate(DstFile);

    //------------------- body ---------------
    wantread := total;  // total muss gesetzt werden !!!
    nread := FileRead(Src,_is, 2*wantread) div 2;

    while (nread = wantread) do
    begin
       //---- pitch extraction ----
       pitch := amdfpitch(pitmin, pitmax, length, _is);

       //---- PICOLA OverLap and ADD stage ----//

       if (rate < 1.0) then
       begin
           ola(pitch, @_is, @_is[pitch]);
           point := pitch;
       end
       else
       begin
          FileWrite(Dst,_is, 2*pitch);
          ola(pitch, @_is[pitch], @_is);
	  point := 0;
       end;

       //---- compensate compansion rate ----*/

       sl := pitch * rcomp;
       lcp := trunc(sl);
       err := err + lcp - sl;

       if (err >= 0.5) then
       begin
          dec(lcp);
          err := err - 1.0;
       end
       else if (err <= -0.5) then
       begin
       	  inc(lcp);
          err := err + 1.0;
       end;
       lproc := lproc + pitch;

       //---- PICOLA Pointer Interval Control (PIC) stage ----*/

	wantread := point + lcp;
	if (wantread > total) then
        begin
           wantread := total - point;
	   FileWrite(Dst,_is[point], 2*wantread);
	   lcp := lcp - wantread;
	   wantread := total;
	   while (lcp > 0) do
           begin
              if (lcp <= total) then
              begin
                 wantread := lcp;
	         nread := FileRead(Src,_is, 2*wantread)div 2;
	         FileWrite(Dst,_is, 2*nread);
	         if (nread <> wantread) then
                     break;
	    	 wantread := total;
	    	 nread := FileRead(Src,_is, 2*wantread)div 2;
              end
              else
              begin
                 nread := FileRead(Src,_is, 2*wantread)div 2;
	         FileWrite(Dst,_is, 2*nread);
	         if (nread <> wantread) then
	         break;
              end;
              lcp := lcp - total;
           end;
        end
        else
        begin
           FileWrite(Dst,_is[point], 2*lcp);
	   point := total - wantread;
	   // shift to next pitch period

           for i := 0 to point-1 do
           begin
              _is[i] := _is[i+wantread];
           end;
           nread := FileRead(Src,_is[point], 2*wantread)div 2;
        end;
    end;

    // write rest */
    FileWrite(Dst,_is, 2*(total - wantread + nread));

    FileClose(Src);
    FileClose(Dst);
end;

{== TMMTimeStretch ============================================================}
constructor TMMTimeStretch.Create(aOwner: TComponent);
begin
   inherited Create(aOwner);