delphi - 算法 smbPitchShift (帕斯卡)

Question

我在Pascal中找了很长时间这个算法并没有找到，我只在C++中找到了它，这令人沮丧。然后我决定将 C++ 代码翻译为 Pascal，但是有一些问题我无法解决。出现错误消息“浮点溢出”。我需要帮助才能使此代码正常工作!

var
  WFX: pWaveFormatEx;

  {** Algoritmo 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

好吧，我通常不这样做，但我也对代码的 Delphi 版本感兴趣，所以我翻译了它。试试my translation看看这是否适合您。

FWIW，还可以看看 Dirac3LE library出自同一作者之手。这是一个更专业的库(PSOLA，不是 WSOLA)，可用于 Windows、Linux、Mac、iPhone 等。刚刚尝试了 Mac 版本，听起来不错。