开发者

Algorithm smbPitchShift (Pascal)

I looked for a long time this algorithm in Pascal and not found, I found it only in C++, it was frustrating. Then I decided to translate the C++ code for Pascal, however there were some problems that I am not able to solve. it appeared an error message "Floating point overflow". I would like help to make this code work!

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));


OK, I usually don't do this, but I also have an interest in having a Delphi version of the code, so I translated it. Try my translation and see if that works for you.

FWIW, also take a look at the Dirac3LE library by the same author. That is a much more professional library (PSOLA, not WSOLA), available for Windows, Linux, Mac, iPhone, etc. Just tried the Mac version and it sounds good.


Note the HINTS that Delphi generates when compiling the code.

For instance I get: "[DCC Hint] Unit1.pas(65): H2077 Value assigned to 'p1' never used" on this line:

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

Because it should really be:

p1^ := p2^;

Also if you add this line at the top of your unit:

{$POINTERMATH ON}

you can do C-style pointer arithmetic so you don't have to use the TmpData workaround. So this:

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

Can be simplified into this:

gInFIFO[gRover]:= InData[i];
0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜