Алгоритм SMBPTHSHIFT (PASCAL)
Я долго искал этот алгоритм в Паскале и не нашел, я нашел его только в C ++, он был разочарован. Затем я решил перевести код C ++ для Pascal, однако были некоторые проблемы, которые я не могу решить. Оказалось сообщение об ошибке «Переполнение плавающей точки». Я хотел бы помочь сделать этот код работать!
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));
5
задан Johan - reinstate Monica 3 June 2014 в 15:39
поделиться
0 ответов
Другие вопросы по тегам: