DSP Routine to perform audio wave pitch shifting

Share your advanced PureBasic knowledge/code with the community.
User avatar
Psychophanta
Always Here
Always Here
Posts: 5153
Joined: Wed Jun 11, 2003 9:33 pm
Location: Anare
Contact:

DSP Routine to perform audio wave pitch shifting

Post by Psychophanta »

Routine For doing pitch shifting While maintaining duration using the Short Time Fourier Transform.
Translated form C++ (see comments at source)

Code: Select all

;/* NAME: smbPitchShift.cpp
; * VERSION: 1.2
; * HOME URL: http://www.dspdimension.com
; * KNOWN BUGS: none
; *
; * SYNOPSIS: Routine For doing pitch shifting While maintaining
; * duration using the Short Time Fourier Transform.
; *
; * DESCRIPTION: The routine takes a pitchShift factor value which is between 0.5
; * (one octave down) And 2. (one octave up). A value of exactly 1 does Not change
; * the pitch. numSampsToProcess tells the routine how many samples in indata[0...
; * numSampsToProcess-1] should be pitch shifted And moved To outdata[0 ...
; * numSampsToProcess-1]. The two buffers can be identical (ie. it can process the
; * Data in-place). fftFrameSize defines the FFT frame size used For the
; * processing. Typical values are 1024, 2048 And 4096. It may be any value <=
; * MAX_FRAME_LENGTH but it MUST be a power of 2. osamp is the STFT
; * oversampling factor which also determines the overlap between adjacent STFT
; * frames. It should at least be 4 For moderate scaling ratios. A value of 32 is
; * recommended For best quality. sampleRate takes the sample rate For the signal 
; * in unit Hz, ie. 44100 For 44.1 kHz audio. The Data passed To the routine in 
; * indata[] should be in the range [-1.0, 1.0), which is also the output range 
; * For the Data, make sure you scale the Data accordingly (For 16bit signed integers
; * you would have To divide (And multiply) by 32768). 
; *
; * COPYRIGHT 1999-2009 Stephan M. Bernsee <smb [AT] dspdimension [DOT] com>
; *
; * 						The Wide Open License (WOL)
; *
; * Permission To use, copy, modify, distribute And sell this software And its
; * documentation For any purpose is hereby granted without fee, provided that
; * the above copyright notice And this license appear in all source copies. 
; * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS Or IMPLIED WARRANTY OF
; * ANY KIND. See http://www.dspguru.com/wol.htm For more information.
; *
; *****************************************************************************/ 
; 
; #include <string.h>
; #include <math.h>
; #include <stdio.h>
; 
; #define M_PI 3.14159265358979323846
; #define MAX_FRAME_LENGTH 8192
; 
; void smbFft(float *fftBuffer, long fftFrameSize, long sign);
; double smbAtan2(double x, double y);
; 
; 
; // -----------------------------------------------------------------------------------------------------------------
; 
; 
; void smbPitchShift(float pitchShift, long numSampsToProcess, long fftFrameSize, long osamp, float sampleRate, float *indata, float *outdata)
; /*
; 	Routine smbPitchShift(). See top of file For explanation
; 	Purpose: doing pitch shifting While maintaining duration using the Short
; 	Time Fourier Transform.
; 	Author: (c)1999-2009 Stephan M. Bernsee <smb [AT] dspdimension [DOT] com>
; */
; {
; 
; 	Static float gInFIFO[MAX_FRAME_LENGTH];
; 	Static float gOutFIFO[MAX_FRAME_LENGTH];
; 	Static float gFFTworksp[2*MAX_FRAME_LENGTH];
; 	Static float gLastPhase[MAX_FRAME_LENGTH/2+1];
; 	Static float gSumPhase[MAX_FRAME_LENGTH/2+1];
; 	Static float gOutputAccum[2*MAX_FRAME_LENGTH];
; 	Static float gAnaFreq[MAX_FRAME_LENGTH];
; 	Static float gAnaMagn[MAX_FRAME_LENGTH];
; 	Static float gSynFreq[MAX_FRAME_LENGTH];
; 	Static float gSynMagn[MAX_FRAME_LENGTH];
; 	Static long gRover = false, gInit = false;
; 	double magn, phase, tmp, window, real, imag;
; 	double freqPerBin, expct;
; 	long i,k, qpd, index, inFifoLatency, stepSize, fftFrameSize2;
; 
; 	/* set up some handy variables */
; 	fftFrameSize2 = fftFrameSize/2;
; 	stepSize = fftFrameSize/osamp;
; 	freqPerBin = sampleRate/(double)fftFrameSize;
; 	expct = 2.*M_PI*(double)stepSize/(double)fftFrameSize;
; 	inFifoLatency = fftFrameSize-stepSize;
; 	If (gRover == false) gRover = inFifoLatency;
; 
; 	/* initialize our Static arrays */
; 	If (gInit == false) {
; 		memset(gInFIFO, 0, MAX_FRAME_LENGTH*SizeOf(float));
; 		memset(gOutFIFO, 0, MAX_FRAME_LENGTH*SizeOf(float));
; 		memset(gFFTworksp, 0, 2*MAX_FRAME_LENGTH*SizeOf(float));
; 		memset(gLastPhase, 0, (MAX_FRAME_LENGTH/2+1)*SizeOf(float));
; 		memset(gSumPhase, 0, (MAX_FRAME_LENGTH/2+1)*SizeOf(float));
; 		memset(gOutputAccum, 0, 2*MAX_FRAME_LENGTH*SizeOf(float));
; 		memset(gAnaFreq, 0, MAX_FRAME_LENGTH*SizeOf(float));
; 		memset(gAnaMagn, 0, MAX_FRAME_LENGTH*SizeOf(float));
; 		gInit = true;
; 	}
; 
; 	/* main processing loop */
; 	For (i = 0; i < numSampsToProcess; i++){
; 
; 		/* As long As we have Not yet collected enough Data just Read in */
; 		gInFIFO[gRover] = indata[i];
; 		outdata[i] = gOutFIFO[gRover-inFifoLatency];
; 		gRover++;
; 
; 		/* now we have enough Data For processing */
; 		If (gRover >= fftFrameSize) {
; 			gRover = inFifoLatency;
; 
; 			/* do windowing And re,im interleave */
; 			For (k = 0; k < fftFrameSize;k++) {
; 				window = -.5*Cos(2.*M_PI*(double)k/(double)fftFrameSize)+.5;
; 				gFFTworksp[2*k] = gInFIFO[k] * window;
; 				gFFTworksp[2*k+1] = 0.;
; 			}
; 
; 
; 			/* ***************** ANALYSIS ******************* */
; 			/* do transform */
; 			smbFft(gFFTworksp, fftFrameSize, -1);
; 
; 			/* this is the analysis Step */
; 			For (k = 0; k <= fftFrameSize2; k++) {
; 
; 				/* de-interlace FFT buffer */
; 				real = gFFTworksp[2*k];
; 				imag = gFFTworksp[2*k+1];
; 
; 				/* compute magnitude And phase */
; 				magn = 2.*sqrt(real*real + imag*imag);
; 				phase = ATan2(imag,real);
; 
; 				/* compute phase difference */
; 				tmp = phase - gLastPhase[k];
; 				gLastPhase[k] = phase;
; 
; 				/* subtract expected phase difference */
; 				tmp -= (double)k*expct;
; 
; 				/* Map delta phase into +/- Pi interval */
; 				qpd = tmp/M_PI;
; 				If (qpd >= 0) qpd += qpd&1;
; 				Else qpd -= qpd&1;
; 				tmp -= M_PI*(double)qpd;
; 
; 				/* get deviation from bin frequency from the +/- Pi interval */
; 				tmp = osamp*tmp/(2.*M_PI);
; 
; 				/* compute the k-th partials' true frequency */
; 				tmp = (double)k*freqPerBin + tmp*freqPerBin;
; 
; 				/* store magnitude And true frequency in analysis arrays */
; 				gAnaMagn[k] = magn;
; 				gAnaFreq[k] = tmp;
; 
; 			}
; 
; 			/* ***************** PROCESSING ******************* */
; 			/* this does the actual pitch shifting */
; 			memset(gSynMagn, 0, fftFrameSize*SizeOf(float));
; 			memset(gSynFreq, 0, fftFrameSize*SizeOf(float));
; 			For (k = 0; k <= fftFrameSize2; k++) { 
; 				index = k*pitchShift;
; 				If (index <= fftFrameSize2) { 
; 					gSynMagn[index] += gAnaMagn[k]; 
; 					gSynFreq[index] = gAnaFreq[k] * pitchShift; 
; 				} 
; 			}
; 			
; 			/* ***************** SYNTHESIS ******************* */
; 			/* this is the synthesis Step */
; 			For (k = 0; k <= fftFrameSize2; k++) {
; 
; 				/* get magnitude And true frequency from synthesis arrays */
; 				magn = gSynMagn[k];
; 				tmp = gSynFreq[k];
; 
; 				/* subtract bin mid frequency */
; 				tmp -= (double)k*freqPerBin;
; 
; 				/* get bin deviation from freq deviation */
; 				tmp /= freqPerBin;
; 
; 				/* take osamp into account */
; 				tmp = 2.*M_PI*tmp/osamp;
; 
; 				/* add the overlap phase advance back in */
; 				tmp += (double)k*expct;
; 
; 				/* accumulate delta phase To get bin phase */
; 				gSumPhase[k] += tmp;
; 				phase = gSumPhase[k];
; 
; 				/* get real And imag part And re-interleave */
; 				gFFTworksp[2*k] = magn*Cos(phase);
; 				gFFTworksp[2*k+1] = magn*Sin(phase);
; 			} 
; 
; 			/* zero negative frequencies */
; 			For (k = fftFrameSize+2; k < 2*fftFrameSize; k++) gFFTworksp[k] = 0.;
; 
; 			/* do inverse transform */
; 			smbFft(gFFTworksp, fftFrameSize, 1);
; 
; 			/* do windowing And add To output accumulator */ 
; 			For(k=0; k < fftFrameSize; k++) {
; 				window = -.5*Cos(2.*M_PI*(double)k/(double)fftFrameSize)+.5;
; 				gOutputAccum[k] += 2.*window*gFFTworksp[2*k]/(fftFrameSize2*osamp);
; 			}
; 			For (k = 0; k < stepSize; k++) gOutFIFO[k] = gOutputAccum[k];
; 
; 			/* shift accumulator */
; 			memmove(gOutputAccum, gOutputAccum+stepSize, fftFrameSize*SizeOf(float));
; 
; 			/* move input FIFO */
; 			For (k = 0; k < inFifoLatency; k++) gInFIFO[k] = gInFIFO[k+stepSize];
; 		}
; 	}
; }
; 
; // -----------------------------------------------------------------------------------------------------------------
; 
; 
; void smbFft(float *fftBuffer, long fftFrameSize, long sign)
; /* 
; 	FFT routine, (C)1996 S.M.Bernsee. Sign = -1 is FFT, 1 is iFFT (inverse)
; 	Fills fftBuffer[0...2*fftFrameSize-1] With the Fourier transform of the
; 	time domain Data in fftBuffer[0...2*fftFrameSize-1]. The FFT Array takes
; 	And returns the cosine And sine parts in an interleaved manner, ie.
; 	fftBuffer[0] = cosPart[0], fftBuffer[1] = sinPart[0], asf. fftFrameSize
; 	must be a power of 2. It expects a complex input signal (see footnote 2),
; 	ie. when working With 'common' audio signals our input signal has To be
; 	passed As {in[0],0.,in[1],0.,in[2],0.,...} asf. In that Case, the transform
; 	of the frequencies of interest is in fftBuffer[0...fftFrameSize].
; */
; {
; 	float wr, wi, arg, *p1, *p2, temp;
; 	float tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
; 	long i, bitm, j, le, le2, k;
; 
; 	For (i = 2; i < 2*fftFrameSize-2; i += 2) {
; 		For (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1) {
; 			If (i & bitm) j++;
; 			j <<= 1;
; 		}
; 		If (i < j) {
; 			p1 = fftBuffer+i; p2 = fftBuffer+j;
; 			temp = *p1; *(p1++) = *p2;
; 			*(p2++) = temp; temp = *p1;
; 			*p1 = *p2; *p2 = temp;
; 		}
; 	}
; 	For (k = 0, le = 2; k < (long)(log(fftFrameSize)/log(2.)+.5); k++) {
; 		le <<= 1;
; 		le2 = le>>1;
; 		ur = 1.0;
; 		ui = 0.0;
; 		arg = M_PI / (le2>>1);
; 		wr = Cos(arg);
; 		wi = sign*Sin(arg);
; 		For (j = 0; j < le2; j += 2) {
; 			p1r = fftBuffer+j; p1i = p1r+1;
; 			p2r = p1r+le2; p2i = p2r+1;
; 			For (i = j; i < 2*fftFrameSize; i += le) {
; 				tr = *p2r * ur - *p2i * ui;
; 				ti = *p2r * ui + *p2i * ur;
; 				*p2r = *p1r - tr; *p2i = *p1i - ti;
; 				*p1r += tr; *p1i += ti;
; 				p1r += le; p1i += le;
; 				p2r += le; p2i += le;
; 			}
; 			tr = ur*wr - ui*wi;
; 			ui = ur*wi + ui*wr;
; 			ur = tr;
; 		}
; 	}
; }
; 
; 
; // -----------------------------------------------------------------------------------------------------------------
; 
; /*
; 
;     12/12/02, smb
;     
;     PLEASE NOTE:
;     
;     There have been some reports on domain errors when the ATan2() function was used
;     As in the above code. Usually, a domain error should Not interrupt the program flow
;     (maybe except in Debug mode) but rather be handled "silently" And a Global variable
;     should be set according To this error. However, on some occasions people ran into
;     this kind of scenario, so a replacement ATan2() function is provided here.
;     
;     If you are experiencing domain errors And your program stops, simply replace all
;     instances of ATan2() With calls To the smbAtan2() function below.
;     
; */
; 
; 
; double smbAtan2(double x, double y)
; {
;   double signx;
;   If (x > 0.) signx = 1.;  
;   Else signx = -1.;
;   
;   If (x == 0.) Return 0.;
;   If (y == 0.) Return signx * M_PI / 2.;
;   
;   Return ATan2(x, y);
; }

;***********************************************************************
;************** PB TRANSLATION by Psychophanta 20130101 ****************
;***********************************************************************
; * NAME: smbPitchShift.cpp
; * VERSION: 1.2
; * HOME URL: http://www.dspdimension.com
; * KNOWN BUGS: none
; *
; * SYNOPSIS: Routine For doing pitch shifting While maintaining
; * duration using the Short Time Fourier Transform.
; *
; * DESCRIPTION: The routine takes a pitchShift factor value which is between 0.5
; * (one octave down) And 2. (one octave up). A value of exactly 1 does Not change
; * the pitch. numSampsToProcess tells the routine how many samples in indata[0...
; * numSampsToProcess-1] should be pitch shifted And moved To outdata[0 ...
; * numSampsToProcess-1]. The two buffers can be identical (ie. it can process the
; * Data in-place). fftFrameSize defines the FFT frame size used For the
; * processing. Typical values are 1024, 2048 And 4096. It may be any value <=
; * MAX_FRAME_LENGTH but it MUST be a power of 2. osamp is the STFT
; * oversampling factor which also determines the overlap between adjacent STFT
; * frames. It should at least be 4 For moderate scaling ratios. A value of 32 is
; * recommended For best quality. sampleRate takes the sample rate For the signal 
; * in unit Hz, ie. 44100 For 44.1 kHz audio. The Data passed To the routine in 
; * indata[] should be in the range [-1.0, 1.0), which is also the output range 
; * For the Data, make sure you scale the Data accordingly (For 16bit signed integers
; * you would have To divide (And multiply) by 32768). 
; *
; * COPYRIGHT 1999-2009 Stephan M. Bernsee <smb [AT] dspdimension [DOT] com>
; *
; * 						The Wide Open License (WOL)
; *
; * Permission To use, copy, modify, distribute And sell this software And its
; * documentation For any purpose is hereby granted without fee, provided that
; * the above copyright notice And this license appear in all source copies. 
; * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS Or IMPLIED WARRANTY OF
; * ANY KIND. See http://www.dspguru.com/wol.htm For more information.
; *
; *****************************************************************************/ 
;#PI=3.14159265358979323846
#MAX_FRAME_LENGTH=8192
Procedure.f float(a.l)
  ProcedureReturn a
EndProcedure
Procedure smbFft(*fftBuffer.float,fftFrameSize.l,sign.l)
; /* 
; 	FFT routine, (C)1996 S.M.Bernsee. Sign=-1 is FFT, 1 is iFFT (inverse)
; 	Fills fftBuffer[0...2*fftFrameSize-1] With the Fourier transform of the
; 	time domain Data in fftBuffer[0...2*fftFrameSize-1]. The FFT Array takes
; 	And returns the cosine And sine parts in an interleaved manner, ie.
; 	fftBuffer[0]=cosPart[0], fftBuffer[1]=sinPart[0], asf. fftFrameSize
; 	must be a power of 2. It expects a complex input signal (see footnote 2),
; 	ie. when working With 'common' audio signals our input signal has To be
; 	passed As {in[0],0.,in[1],0.,in[2],0.,...} asf. In that Case, the transform
; 	of the frequencies of interest is in fftBuffer[0...fftFrameSize].
; */
	Protected wr.f,wi.f,arg.f,*p1.float,*p2.float,temp.f
	Protected tr.f,ti.f,ur.f,ui.f,*p1r.float,*p1i.float,*p2r.float,*p2i.float
	Protected i.l,bitm.l,j.l,le.l,le2.l,k.l
	i=2
  While i<2*fftFrameSize-2
    j=0:bitm=2
		While bitm<2*fftFrameSize
			If i&bitm:j+1:EndIf
			j<<1
			bitm<<1
		Wend
		If i<j
			p1=*fftBuffer+i:p2=*fftBuffer+j
			temp=*p1:*p1+SizeOf(float):*p1=*p2
			*p2+SizeOf(float):*p2=temp:temp=*p1
			*p1=*p2:*p2=temp
		EndIf
		i+2
	Wend
  le=2:k=0
	While k<IntQ(Log(fftFrameSize)/Log(2.0)+0.5)
		le<<1
		le2=le>>1
		ur=1.0
		ui=0.0
		arg=#PI/float(le2>>1)
		wr=Cos(arg)
		wi=sign*Sin(arg)
		j=0
		While j<le2
			p1r=*fftBuffer+j:p1i=p1r+1
			p2r=p1r+le2:p2i=p2r+1
			i=j
			While i<2*fftFrameSize
				tr=*p2r*ur-*p2i*ui
				ti=*p2r*ui+*p2i*ur
				*p2r=*p1r-tr:*p2i=*p1i-ti
				*p1r+tr:*p1i+ti
				p1r+le:p1i+le
				p2r+le:p2i+le
				i+le
			Wend
			tr=ur*wr-ui*wi
			ui=ur*wi+ui*wr
			ur=tr
			j+2
		Wend
		k+1
	Wend
EndProcedure
Procedure.d smbAtan2(x.d,y.d)
; 
; // -----------------------------------------------------------------------------------------------------------------
; 
; /*
; 
;     12/12/02, smb
;     
;     PLEASE NOTE:
;     
;     There have been some reports on domain errors when the ATan2() function was used
;     As in the above code. Usually, a domain error should Not interrupt the program flow
;     (maybe except in Debug mode) but rather be handled "silently" And a Global variable
;     should be set according To this error. However, on some occasions people ran into
;     this kind of scenario, so a replacement ATan2() function is provided here.
;     
;     If you are experiencing domain errors And your program stops, simply replace all
;     instances of ATan2() With calls To the smbAtan2() function below.
;     
; */
; 
; 
  Protected signx.d
  If x>0.0:signx=1.0:Else:signx=-1.0:EndIf
  If x=0.0:ProcedureReturn 0.0:EndIf
  If y=0.0:ProcedureReturn signx*#PI/2.0:EndIf
  ProcedureReturn ATan2(x,y)
EndProcedure
Procedure smbPitchShift(pitchShift.f,numSampsToProcess.l,fftFrameSize.l,osamp.l,sampleRate.f,Array indata.f(1),Array outdata.f(1))
; /*
; 	Routine smbPitchShift(). See top of file For explanation
; 	Purpose: doing pitch shifting While maintaining duration using the Short
; 	Time Fourier Transform.
; 	Author: (c)1999-2009 Stephan M. Bernsee <smb [AT] dspdimension [DOT] com>
; */
	Static Dim gInFIFO.f(#MAX_FRAME_LENGTH)
	Static Dim gOutFIFO.f(#MAX_FRAME_LENGTH)
	Static Dim gFFTworksp.f(2*#MAX_FRAME_LENGTH)
	Static Dim gLastPhase.f(#MAX_FRAME_LENGTH/2+1)
	Static Dim gSumPhase.f(#MAX_FRAME_LENGTH/2+1)
	Static Dim gOutputAccum.f(2*#MAX_FRAME_LENGTH)
	Static Dim gAnaFreq.f(#MAX_FRAME_LENGTH)
	Static Dim gAnaMagn.f(#MAX_FRAME_LENGTH)
	Static Dim gSynFreq.f(#MAX_FRAME_LENGTH)
	Static Dim gSynMagn.f(#MAX_FRAME_LENGTH)
	Static gRover.l=#False,gInit.l=#False
	Protected magn.d,phase.d,tmp.d,window.d,real.d,imag.d
	Protected freqPerBin.d,expct.d
	Protected i.l,k.l,qpd.l,index.l,inFifoLatency.l,stepSize.l,fftFrameSize2.l
; 	/* set up some handy variables */
	fftFrameSize2=fftFrameSize/2
	stepSize=fftFrameSize/osamp
	freqPerBin=sampleRate/fftFrameSize
	expct=2.0*#PI*stepSize/fftFrameSize
	inFifoLatency=fftFrameSize-stepSize
	If gRover=#False:gRover=inFifoLatency:EndIf
; 	/* initialize our Static arrays */
  If gInit=#False
    FillMemory(gInFIFO(),#MAX_FRAME_LENGTH*SizeOf(float))
		FillMemory(gOutFIFO(),#MAX_FRAME_LENGTH*SizeOf(float))
		FillMemory(gFFTworksp(),2*#MAX_FRAME_LENGTH*SizeOf(float))
		FillMemory(gLastPhase(),(#MAX_FRAME_LENGTH/2+1)*SizeOf(float))
		FillMemory(gSumPhase(),(#MAX_FRAME_LENGTH/2+1)*SizeOf(float))
		FillMemory(gOutputAccum(),2*#MAX_FRAME_LENGTH*SizeOf(float))
		FillMemory(gAnaFreq(),#MAX_FRAME_LENGTH*SizeOf(float))
		FillMemory(gAnaMagn(),#MAX_FRAME_LENGTH*SizeOf(float))
		gInit=#True
	EndIf
; 	/* main processing loop */
  i=0
	While i<numSampsToProcess
; 	/* As long As we have Not yet collected enough Data just Read in */
		gInFIFO(gRover)=indata(i)
		outdata(i)=gOutFIFO(gRover-inFifoLatency)
		gRover+1
; 		/* now we have enough Data For processing */
		If gRover>=fftFrameSize
			gRover=inFifoLatency
; 			/* do windowing And re,im interleave */
			k=0
      While k<fftFrameSize
 				window=-0.5*Cos(2.0*#PI*k/fftFrameSize)+0.5
				gFFTworksp(2*k)=gInFIFO(k)*window
				gFFTworksp(2*k+1)=0.0
      k+1
      Wend
    EndIf
;		/* ***************** ANALYSIS ******************* */
; 	/* do transform */
		smbFft(gFFTworksp,fftFrameSize,-1)
;		/* this is the analysis Step */
    k=0
    While k<=fftFrameSize2
;			/* de-interlace FFT buffer */
			real=gFFTworksp(2*k)
			imag=gFFTworksp(2*k+1)
;			/* compute magnitude And phase */
			magn=2.0*Sqr(real*real+imag*imag)
			phase=ATan2(imag,real)
;			/* compute phase difference */
			tmp=phase-gLastPhase(k)
			gLastPhase(k)=phase
;			/* subtract expected phase difference */
			tmp-k*expct
;			/* Map delta phase into +/- Pi interval */
			qpd=tmp/#PI
			If qpd>=0:qpd+qpd&1:Else:qpd-qpd&1:EndIf
			tmp-#PI*qpd
;			/* get deviation from bin frequency from the +/- Pi interval */
			tmp=osamp*tmp/(2.0*#PI)
;			/* compute the k-th partials' true frequency */
			tmp=k*freqPerBin+tmp*freqPerBin
;			/* store magnitude And true frequency in analysis arrays */
			gAnaMagn(k)=magn
			gAnaFreq(k)=tmp
			k+1
		Wend
; 	/* ***************** PROCESSING ******************* */
; 	/* this does the actual pitch shifting */
		FillMemory(gSynMagn,fftFrameSize*SizeOf(float))
		FillMemory(gSynFreq,fftFrameSize*SizeOf(float))
		k=0
		While k<=fftFrameSize2
			index=k*pitchShift
			If index<=fftFrameSize2
				gSynMagn(index)+gAnaMagn(k)
				gSynFreq(index)=gAnaFreq(k)*pitchShift
			EndIf
		  k+1
		Wend
; 	/* ***************** SYNTHESIS ******************* */
; 	/* this is the synthesis Step */
    k=0
    While k<=fftFrameSize2
; 			/* get magnitude And true frequency from synthesis arrays */
			magn=gSynMagn(k)
			tmp=gSynFreq(k)
; 			/* subtract bin mid frequency */
			tmp-k*freqPerBin
; 			/* get bin deviation from freq deviation */
			tmp/freqPerBin
; 			/* take osamp into account */
			tmp=2.0*#PI*tmp/osamp
; 			/* add the overlap phase advance back in */
			tmp+k*expct
; 			/* accumulate delta phase To get bin phase */
			gSumPhase(k)+tmp
			phase=gSumPhase(k)
; 			/* get real And imag part And re-interleave */
			gFFTworksp(2*k)=magn*Cos(phase)
			gFFTworksp(2*k+1)=magn*Sin(phase)
			k+1
		Wend
; 		/* zero negative frequencies */
		k=fftFrameSize+2:While k<2*fftFrameSize:gFFTworksp(k)=0.0:k+1:Wend
; 		/* do inverse transform */
		smbFft(gFFTworksp,fftFrameSize,1)
; 		/* do windowing And add To output accumulator */ 
    k=0
    While k<fftFrameSize
			window=-0.5*Cos(2.0*#PI*k/fftFrameSize)+0.5
			gOutputAccum(k)+2.0*window*gFFTworksp(2*k)/(fftFrameSize2*osamp)
			k+1
		Wend
		k=0:While k<stepSize:gOutFIFO(k)=gOutputAccum(k):k+1:Wend
; 		/* shift accumulator */
		MoveMemory(gOutputAccum+stepSize,gOutputAccum,fftFrameSize*SizeOf(float))
; 		/* move input FIFO */
    k=0
    While k<inFifoLatency:gInFIFO(k)=gInFIFO(k+stepSize):k+1:Wend
  i+1
  Wend
EndProcedure 
;\
http://www.zeitgeistmovie.com

while (world==business) world+=mafia;
User avatar
einander
Enthusiast
Enthusiast
Posts: 744
Joined: Thu Jun 26, 2003 2:09 am
Location: Spain (Galicia)

Re: DSP Routine to perform audio wave pitch shifting

Post by einander »

Thanks for this!
Exactly what I need, as I'm experimenting with sound samples.
Cheers!
User avatar
Psychophanta
Always Here
Always Here
Posts: 5153
Joined: Wed Jun 11, 2003 9:33 pm
Location: Anare
Contact:

Re: DSP Routine to perform audio wave pitch shifting

Post by Psychophanta »

Waiting for some good tip from you... :)
http://www.zeitgeistmovie.com

while (world==business) world+=mafia;
User avatar
oryaaaaa
Addict
Addict
Posts: 825
Joined: Mon Jan 12, 2004 11:40 pm
Location: Okazaki, JAPAN

Re: DSP Routine to perform audio wave pitch shifting

Post by oryaaaaa »

Very nice. I feel your sound processing code.
But, more speed and more preicision, then I want to need 32-bit integer mode.

When the USB-DAC 32-bit any one appeared, I hate being overtaken by others.
Post Reply