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];
加载中,请稍侯......
精彩评论