Algorytm smbPitchShift (Pascal)

Question

Długo szukałem tego algorytmu w Pascalu i nie znalazłem, znalazłem go tylko w C++, to było frustrujące. Potem zdecydowałem się przetłumaczyć kod C++ dla Pascala, jednak były pewne problemy, których nie jestem w stanie rozwiązać. pojawił się komunikat o błędzie "Przepełnienie zmiennoprzecinkowe". Chciałbym pomóc, aby ten kod działał!

var 
    WFX: pWaveFormatEx; 

    {** Algoritimo Pitch Shift **} 
    gInFIFO, gOutFIFO, gLastPhase, gSumPhase, gOutputAccum: Array Of Extended; 
    gAnaMagn, gAnaFreq, gSynFreq, gSynMagn, gFFTworksp: Array Of Extended; 

Const 
    MAX_FRAME_LENGTH = 8192; 

implementation 

{$R *.dfm} 

procedure smbFft(fftBuffer: PExtended; fftFrameSize, sign: Integer); 
var 
    p1, p2, p1r, p1i, p2r, p2i: PExtended; 
    wr, wi, arg, temp: EXTENDED; 
    tr, ti, ur, ui: EXTENDED; 
    i, bitm, j, le, le2, k: Integer; 
begin 
    i:= 2; 
    WHILE (i < 2*fftFrameSize-2) DO            //for (i = 2; i < 2*fftFrameSize-2; i += 2) { 
    BEGIN 
     bitm:= 2; 
     j:= 0; 
     WHILE (bitm < (2 * fftFrameSize)) DO          //for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1) { 
     BEGIN 
     if ((i and bitm) <> 0) then           //if (i & bitm) j++; 
      inc(j); 
      // 
      j:= j shl 1;               //j <<= 1; 
      bitm:= bitm shl 1;              //bitm <<= 1 
    END; 
     // 
     if (i < j) then 
     begin 
     p1:= fftBuffer;              //^ 
      Inc(p1, i);               //p1 = fftBuffer+i; 
      p2:= fftBuffer;              //^ 
      Inc(p2, j);               //p2 = fftBuffer+j; 
     temp:= p1^;               //temp = *p1; 
      inc(p1, 1);               //*(p1++) 
      p1:= p2;                //p1 = *p2; 
      inc(p2, 1);               //*(p2++) 
      p2^:= temp;               //p2 = temp; 
      temp:= p1^;               //temp = *p1; 
     p1:= p2;                //*p1 = *p2; 
      p2^:= temp;               //*p2 = temp; 
    end; 
     INC(I, 2); 
    END; 
    // 
    le:= 2; 
    k:= 0; 
    WHILE (k < (ln(fftFrameSize)/ln(2.0)+0.5)) DO         //for (k = 0, le = 2; k < (long)(log(fftFrameSize)/log(2.)+.5); k++) { 
    BEGIN 
     le:= le shl 1;                //le <<= 1; 
     le2:= le shr 1;               //le2 = le>>1; 
     ur:= 1.0;                 //ur = 1.0; 
     ui:= 0.0;                 //ui = 0.0; 
     arg:= PI/(le2 shr 1);             //arg = M_PI/(le2>>1); 
     wr:= cos(arg);                //wr = cos(arg); 
     wi:= sign * sin(arg);              //wi = sign*sin(arg); 
     j:=0; 
     WHILE (j < le2) DO               //for (j = 0; j < le2; j += 2) { 
     BEGIN 
      p1r:= fftBuffer;              //^ 
      INC(p1r, j);               //p1r = fftBuffer+j; 
      p1i:= p1r;                //^ 
      INC(p1i, 1);               //p1i = p1r+1; 
      p2r:= p1r;                //^ 
      INC(p2r, le2);               //p2r = p1r+le2; 
      p2i:= p2r;                //^ 
      INC(p2i, 1);               //p2i = p2r+1; 
      i:= j; 
     WHILE (i < 2*fftFrameSize) DO           //for (i = j; i < 2*fftFrameSize; i += le) { 
      BEGIN 
      tr:= p2r^ * ur - p2i^ * ui;          //tr = *p2r * ur - *p2i * ui; 
      ti:= p2r^ * ui + p2i^ * ur;          //ti = *p2r * ui + *p2i * ur; 
      p2r^:= p1r^ - tr;             //*p2r = *p1r - tr; 
       p2i^:= p1i^ - ti;             //*p2i = *p1i - ti; 
      p1r^:= p1r^ + tr;             //*p1r += tr; 
       p1i^:= p1i^ + ti;             //*p1i += ti; 
      INC(p1r, le);              //p1r += le; 
       INC(p1i, le);              //p1i += le; 
       INC(p2r, le);              //p2r += le; 
       INC(p2i, le);              //p2i += le; 
       INC(i, le); 
     END; 
      // 
      tr:= ur * wr - ui * wi;            //tr = ur*wr - ui*wi; 
     ui:= ur * wi + ui * wr;            //ui = ur*wi + ui*wr; 
      ur:= tr;                //ur = tr; 
      INC(J, 2); 
      END; 
     inc(k); 
    END; 
end; 

Procedure smbPitchShift(pitchShift: Double; numSampsToProcess, fftFrameSize, osamp, sampleRate: Integer; indata, outdata: PExtended); 

    function atan2 (y, x : Extended) : Extended; Assembler; 
    asm 
    fld [y] 
    fld [x] 
    fpatan 
    end; 
var magn, phase, tmp, window, xreal, imag: Extended; 
    freqPerBin, expct, CC: Extended; 
    i, k, qpd, index, inFifoLatency, stepSize, fftFrameSize2: Integer; 
    gRover: Integer; 
    TmpData: PExtended; 
begin 
    gRover:= 0; 
    {* set up some handy variables *} 
    fftFrameSize2:= Round(fftFrameSize/2);          //fftFrameSize2 = fftFrameSize/2; 
    stepSize:= Round(fftFrameSize/osamp);          //stepSize = fftFrameSize/osamp; 
    freqPerBin:= sampleRate/fftFrameSize;          //freqPerBin = sampleRate/(double)fftFrameSize; 
    expct:= 2.0 * PI * stepSize/fftFrameSize;         //expct = 2.*M_PI*(double)stepSize/(double)fftFrameSize; 
    inFifoLatency:= fftFrameSize - stepSize;          //inFifoLatency = fftFrameSize-stepSize; 
    if (gRover = 0) then gRover:= inFifoLatency;         //if (gRover == false) gRover = inFifoLatency; 
    // 
    {* main processing loop *} 
    for i:=0 to numSampsToProcess-1 do            //for (i = 0; i < numSampsToProcess; i++){ 
    begin 
     {* As long as we have not yet collected enough data just read in *} 
     TmpData:= indata;               //^ 
     inc(TmpData, i);               // [i] 
     gInFIFO[gRover]:= TmpData^;            //gInFIFO[gRover] = indata[i]; 
     TmpData:= outdata;               //^ 
     inc(TmpData, i);               // [i] 
     TmpData^:= gOutFIFO[gRover - inFifoLatency];        //outdata[i] = gOutFIFO[gRover-inFifoLatency]; 
     Inc(gRover);                //gRover++; 
     {* now we have enough data for processing *} 
     if (gRover >= fftFrameSize) then           //if (gRover >= fftFrameSize) { 
      begin 
      gRover:= inFifoLatency;            //gRover = inFifoLatency; 
      {* do windowing and re,im interleave *} 
      for k:=0 to fftFrameSize-1 do          //for (k = 0; k < fftFrameSize;k+ 
       begin 
       window:= -0.5 * Cos(2.0 * PI * k/fftFrameSize) + 0.5;   //window = -.5*cos(2.*M_PI*(double)k/(double)fftFrameSize)+.5; 
       gFFTworksp[2 * k]:= gInFIFO[k] * window;       //gFFTworksp[2*k] = gInFIFO[k] * window; 
       gFFTworksp[2 * k + 1]:= 0.0;          //gFFTworksp[2 * k + 1]:= 0.0F; 
       end; 
      {****************** ANALYSIS ********************} 
      {* do transform *} 
      SmbFft(Ptr(DWORD(gFFTworksp)), fftFrameSize, -1);     //smbFft(gFFTworksp, fftFrameSize, -1); 
      {* this is the analysis step *} 
      for k:= 0 to fftFrameSize2 do          //for (k = 0; k <= fftFrameSize2; k++) { 
       begin 
       {* de-interlace FFT buffer *} 
       xreal:= gFFTworksp[2 * k];          //real = gFFTworksp[2*k]; 
       imag:= gFFTworksp[2 * k + 1];         //imag = gFFTworksp[2*k+1]; 
       {* compute magnitude and phase *} 
       magn:= 2.0 * Sqrt(xreal * xreal + imag * imag);     //magn = 2.*sqrt(real*real + imag*imag); 
       phase:= Atan2(imag, xreal);          //phase = atan2(imag,real); 
       {* compute phase difference *} 
       tmp:= phase - gLastPhase[k];          //tmp = phase - gLastPhase[k]; 
       gLastPhase[k]:= phase;           //gLastPhase[k] = phase; 
       {* subtract expected phase difference *} 
       tmp:= tmp - k * expct;           //tmp -= (double)k*expct; 
       {* map delta phase into +/- Pi interval *} 
       qpd:= Round(tmp/PI);           //qpd = tmp/M_PI; 
       if (qpd >= 0) then 
        qpd:= qpd + qpd and 1           // if (qpd >= 0) qpd += qpd&1; 
       else 
        qpd:= qpd - qpd and 1;           // else qpd -= qpd&1; 
       // 
       tmp:= tmp - (PI * qpd);           //tmp -= M_PI*(double)qpd; 
       {* get deviation from bin frequency from the +/- Pi interval *} 
       tmp:= osamp * tmp/(2.0 * PI);         //tmp = osamp*tmp/(2.*M_PI); 
       {* compute the k-th partials' true frequency *} 
       tmp:= k * freqPerBin + tmp * freqPerBin;       //tmp = (double)k*freqPerBin + tmp*freqPerBin; 
       {* store magnitude and true frequency in analysis arrays *} 
       gAnaMagn[k]:= magn;            //gAnaMagn[k] = magn; 
       gAnaFreq[k]:= tmp;            //gAnaFreq[k] = tmp; 
       end; 
      {****************** PROCESSING ********************} 
      {* this does the actual pitch shifting *} 
      for k:=0 to fftFrameSize2 do           //for (k = 0; k <= fftFrameSize2; k++) { 
       begin 
       index:= Round(k * pitchShift);         //index = (long)(k*pitchShift); 
       if (index <= fftFrameSize2) then         //if (index <= fftFrameSize2) { 
        begin 
        IF K >= LENGTH(gSynFreq) THEN 
         SetLength(gSynFreq , LENGTH(gSynFreq)+1);     //memset(gSynFreq, 0, fftFrameSize*sizeof(float)); 
        IF K >= LENGTH(gSynMagn) THEN 
         SetLength(gSynMagn , LENGTH(gSynMagn)+1);     //memset(gSynMagn, 0, fftFrameSize*sizeof(float)); 
        // 
        gSynMagn[index]:= gSynMagn[index] + gAnaMagn[k];    //gSynMagn[index] += gAnaMagn[k]; 
        gSynFreq[index]:= gAnaFreq[k] * pitchShift;     //gSynFreq[index] = gAnaFreq[k] * pitchShift; 
        end; 
       end; 
      {****************** SYNTHESIS ********************} 
      {* this is the synthesis step *} 
      for k:=0 to fftFrameSize2 do           //for (k = 0; k <= fftFrameSize2; k++) { 
       begin 
       {* get magnitude and true frequency from synthesis arrays *} 
       magn:= gSynMagn[k];            // magn = gSynMagn[k]; 
       tmp:= gSynFreq[k];            //tmp = gSynFreq[k] 
       {* subtract bin mid frequency *} 
       tmp:= tmp - (k * freqPerBin);         //tmp -= (double)k*freqPerBin; 
       {* get bin deviation from freq deviation *} 
       tmp:= tmp/freqPerBin;           //tmp /= freqPerBin; 
       {* take osamp into account *} 
       tmp:= 2.0 * PI * tmp/osamp;         //tmp = 2.*M_PI*tmp/osamp; 
       {* add the overlap phase advance back in *} 
       tmp:= tmp + (k * expct);           //tmp += (double)k*expct; 
       {* accumulate delta phase to get bin phase *} 
       gSumPhase[k]:= gSumPhase[k] + tmp;        //gSumPhase[k] += tmp; 
       phase:= gSumPhase[k];           //phase = gSumPhase[k]; 
       {* get real and imag part and re-interleave *} 
       gFFTworksp[2 * k]:= (magn * Cos(phase));       //gFFTworksp[2*k] = magn*cos(phase); 
       gFFTworksp[2 * k + 1]:= (magn * Sin(phase));      //gFFTworksp[2*k+1] = magn*sin(phase); 
       end; 
      {* zero negative frequencies *} 
      k:= fftFrameSize + 2; 
      WHILE (k < 2 * fftFrameSize) DO          //for (k = fftFrameSize+2; k < 2*fftFrameSize; k++) 
       BEGIN 
       gFFTworksp[k]:= 0.0;            //gFFTworksp[k] = 0.0F; 
       inc(k); 
       END; 
      {* do inverse transform *} 
      SmbFft(Ptr(DWORD(gFFTworksp)), fftFrameSize, 1);      //smbFft(gFFTworksp, fftFrameSize, 1); 
      {* do windowing and add to output accumulator *} 
      for k:=0 to fftFrameSize-1 do          // for(k=0; k < fftFrameSize; k++) { 
       begin 
       window:= -0.5 * Cos(2.0 * PI * k/fftFrameSize) + 0.5;   //window = -.5*cos(2.*M_PI*(double)k/(double)fftFrameSize)+.5; 
       gOutputAccum[k]:= gOutputAccum[k] + (2.0 * window * gFFTworksp[2 * k]/(fftFrameSize2 * osamp)); 
       end;                //gOutputAccum[k] += 2.*window*gFFTworksp[2*k]/(fftFrameSize2*osamp); 
      // 
      for k:=0 to stepSize-1 do gOutFIFO[k]:= gOutputAccum[k];    //for (k = 0; k < stepSize; k++) gOutFIFO[k] = gOutputAccum[k]; 
      {* shift accumulator *} 
      // 
      TmpData:= PTR(DWORD(gOutputAccum));         //^ 
      Inc(TmpData, StepSize);            //gOutputAccum + stepSize 
      MoveMemory(TmpData, PTR(DWORD(gOutputAccum)), fftFrameSize * sizeof(Extended)); 
                      //memmove(gOutputAccum, gOutputAccum + stepSize, fftFrameSize * sizeof(float)); 
      // 
      {* move input FIFO *} 
      for k:=0 to inFifoLatency-1 do          //for (k = 0; k < inFifoLatency; k++) 
       gInFIFO[k]:= gInFIFO[k + stepSize];        //gInFIFO[k] = gInFIFO[k+stepSize]; 
      end; 
    end; 
end; 

procedure TWavAnalize.FormCreate(Sender: TObject); 
begin 
    {** algoritimo pitchshift **} 
    SetLength(gInFIFO ,MAX_FRAME_LENGTH); 
    SetLength(gOutFIFO ,MAX_FRAME_LENGTH); 
    SetLength(gSynFreq ,MAX_FRAME_LENGTH); 
    SetLength(gSynMagn ,MAX_FRAME_LENGTH); 
    SetLength(gAnaFreq ,MAX_FRAME_LENGTH); 
    SetLength(gAnaMagn ,MAX_FRAME_LENGTH); 
    SetLength(gFFTworksp ,2 * MAX_FRAME_LENGTH); 
    SetLength(gLastPhase , Round(MAX_FRAME_LENGTH/2) + 1); 
    SetLength(gSumPhase , Round(MAX_FRAME_LENGTH/2) + 1); 
    SetLength(gOutputAccum , 2 * MAX_FRAME_LENGTH); 
    {** algoritimo pitchshift **} 
end; 

procedure TWavAnalize.Button8Click(Sender: TObject); 
VAR T: TMEMORYSTREAM; 
    DSize, DataOffset, i: cARDINAL; 
    AIN, AOUT: ARRAY OF EXTENDED; 
begin 
    T:= TMEMORYSTREAM.CREATE; 
    T.LoadFromFile(PATH); 
    GetStreamWaveAudioInfo(T, WFX, DSize, DataOffset); 
    T.Position:= DataOffset; 
    SETLENGTH(AIN, DSIZE); 
    SETLENGTH(AOUT, DSIZE); 
    T.READ(AIN[0], DSIZE); 
    smbPitchShift(0.5, DSize, 2048, 10, WFX.nSamplesPerSec, Ptr(DWORD(AIN)), Ptr(DWORD(AOUT))); 
    T.Clear; 
    T.WRITE(AOUT[0], LENGTH(AOUT)); 

Answer 1

Zwróć uwagę na HINTS generowany przez Delphi podczas kompilowania kodu.

Na przykład mam: "[DCC Hint] Unit1.pas (65): H2077 Wartość przypisana do 'P1' nigdy nie używane" na tej linii:

p1:= p2;                //*p1 = *p2; 

ponieważ powinien być naprawdę:

p1^ := p2^; 

także jeśli dodać tę linię na górnej części urządzenia:

{$POINTERMATH ON} 

można zrobić w stylu C pointe r arytmetyczna, więc nie musisz używać obejścia TmpData. Więc tak:

TmpData:= indata;             
inc(TmpData, i);              
gInFIFO[gRover]:= TmpData^;  

można uprościć w ten sposób:

gInFIFO[gRover]:= InData[i]; 

Answer 2

OK, ja zwykle nie rób tego, ale również mają interes w wersji Delphi kodu, więc przetłumaczyć to. Wypróbuj my translation i sprawdź, czy to działa.

FWIW, również spójrz na Dirac3LE library tego samego autora. To jest o wiele bardziej profesjonalna biblioteka (PSOLA, nie WSOLA) dostępna dla systemów Windows, Linux, Mac, iPhone itp. Po prostu wypróbował wersję na Maca i brzmi dobrze.