Fast Convolution

[fromthisdrive]

Ok, so im getting really stuck with this now. I am writing a VST plugin to perform microphone modelling using fast convolution. I thought i had my FFT working fine, but now im not sure. I'm using Sorensens In Place Split Radix FFT for real values which i found here http://www.musicdsp.org/files/rvfft.cpp . I presumed (before looking at how it worked) that i could perform this in the processreplacing function and it would output the FFT to another array. However it acts directly on the data and what i think im getting output doesnt look like it could be the bins of an FFT. When i perform an inverse FFT before outputing the audio is back there fine, so i presume the algorithms are working, im just not sure what i should be getting in between. Should this algorithm give me a set of bins with real values between 0 and 1 representing the magnitude of each frequency?

I have included below my processreplacing function and the sorensen fft:

Code: Select all

void VST_Plug_in::processReplacing (float** inputs, float** outputs, VstInt32 sampleFrames)
{
	

    float* in1  =  inputs[0];
    float* in2  =  inputs[1];
    float* out1 = outputs[0];
    float* out2 = outputs[1];

	float leftSample  = 0.0;
	float rightSample = 0.0;
	float impulseSample = 0.0;

	realfft_split(inputs[0] , 1024);
	realfft_split(inputs[1] , 1024);

	
	for(int i = 0; i < sampleFrames; i++)
    {
		//get samples from input buffer
		leftSample = (*in1++);
		rightSample = (*in2++);
		
		//do stuff here
		
		

		// write samples to output buffer
        (*out1++) = leftSample;
        (*out2++) = rightSample;
    }

	//irealfft_split(outputs[0], 1024);  //commented out to see what is at output when 
	//irealfft_split(outputs[1], 1024);  //inverse fft is not performed
}

Code: Select all

/////////////////////////////////////////////////////////
// Sorensen in-place split-radix FFT for real values
// data: array of floats:
// re(0),re(1),re(2),...,re(size-1)
// 
// output:
// re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1)
// normalized by array length
//
// Source: 
// Sorensen et al: Real-Valued Fast Fourier Transform Algorithms,
// IEEE Trans. ASSP, ASSP-35, No. 6, June 1987

void realfft_split(float *data,long n){

  long i,j,k,i5,i6,i7,i8,i0,id,i1,i2,i3,i4,n2,n4,n8;
  float t1,t2,t3,t4,t5,t6,a3,ss1,ss3,cc1,cc3,a,e,sqrt2;
  
  sqrt2=sqrt(2.0);
  n4=n-1;
  
  //data shuffling
      for (i=0,j=0,n2=n/2; i<n4 ; i++){
	  if (i<j){
				t1=data[j];
				data[j]=data[i];
				data[i]=t1;
				}
	  k=n2;
	  while (k<=j){
				j-=k;
				k>>=1;	
				}
	  j+=k;
      }
	
/*----------------------*/
	
	//length two butterflies	
	i0=0;
	id=4;
   do{
       for (; i0<n4; i0+=id){ 
			i1=i0+1;
			t1=data[i0];
			data[i0]=t1+data[i1];
			data[i1]=t1-data[i1];
		}
	   id<<=1;
       i0=id-2;
       id<<=1;
    } while ( i0<n4 );

   /*----------------------*/
   //L shaped butterflies
n2=2;
for(k=n;k>2;k>>=1){  
	n2<<=1;
	n4=n2>>2;
	n8=n2>>3;
	e = 2*M_PI/(n2);
	i1=0;
	id=n2<<1;
	do{ 
	    for (; i1<n; i1+=id){
			i2=i1+n4;
			i3=i2+n4;
			i4=i3+n4;
			t1=data[i4]+data[i3];
			data[i4]-=data[i3];
			data[i3]=data[i1]-t1;
			data[i1]+=t1;
			if (n4!=1){
				i0=i1+n8;
				i2+=n8;
				i3+=n8;
				i4+=n8;
				t1=(data[i3]+data[i4])/sqrt2;
				t2=(data[i3]-data[i4])/sqrt2;
				data[i4]=data[i2]-t1;
				data[i3]=-data[i2]-t1;
				data[i2]=data[i0]-t2;
				data[i0]+=t2;
			}
	     }
		 id<<=1;
	     i1=id-n2;
	     id<<=1;
	  } while ( i1<n );
	a=e;
	for (j=2; j<=n8; j++){  
	      a3=3*a;
	      cc1=cos(a);
	      ss1=sin(a);
	      cc3=cos(a3);
	      ss3=sin(a3);
	      a=j*e;
	      i=0;
	      id=n2<<1;
	      do{
		   for (; i<n; i+=id){  
			  i1=i+j-1;
			  i2=i1+n4;
			  i3=i2+n4;
			  i4=i3+n4;
			  i5=i+n4-j+1;
			  i6=i5+n4;
			  i7=i6+n4;
			  i8=i7+n4;
			  t1=data[i3]*cc1+data[i7]*ss1;
			  t2=data[i7]*cc1-data[i3]*ss1;
			  t3=data[i4]*cc3+data[i8]*ss3;
			  t4=data[i8]*cc3-data[i4]*ss3;
			  t5=t1+t3;
			  t6=t2+t4;
			  t3=t1-t3;
			  t4=t2-t4;
			  t2=data[i6]+t6;
			  data[i3]=t6-data[i6];
			  data[i8]=t2;
			  t2=data[i2]-t3;
			  data[i7]=-data[i2]-t3;
			  data[i4]=t2;
			  t1=data[i1]+t5;
			  data[i6]=data[i1]-t5;
			  data[i1]=t1;
			  t1=data[i5]+t4;
			  data[i5]-=t4;
			  data[i2]=t1;
		     }
		   id<<=1;
		   i=id-n2;
		   id<<=1;
		 } while(i<n);
	   }
      }

	//division with array length
   for(i=0;i<n;i++) data[i]/=n;
}

[/code]

[pj geerlings]

there is no gaurentee that your input buffers will contain 1024 valid samples

[fromthisdrive]

what if i used sampleframes as the length of my fft window?

[MadBrain]

Some hosts (buzz, psycle) can give you a different number of samples at each call. In general, it's better to record samples to a buffer, and FFT when that buffer is full. Also,do you do windowing?

[fromthisdrive]

not got as far as windowing yet, i just want to get things working on a basic level before i take it further. but to be honest ive got myself a bit lost with it all. Should i do some kind of overlap add windowing?

If i do store the incoming audio samples into a buffer and wait till its full before outputting, wont i get some kind of latency issues?

[Borogove]

If i do store the incoming audio samples into a buffer and wait till its full before outputting, wont i get some kind of latency issues?

Yes, FFT always has latency issues.

[fromthisdrive]

Ok, cheers for all the help so far. Might just go for making it a process rather than a real time plug-in.

Am i right in thinking the memory pointed to by *data from this fft algorithm should be consecutive real numbers representing the magnitude of the frequency at each point?

[Borogove]

Ok, cheers for all the help so far. Might just go for making it a process rather than a real time plug-in.

Am i right in thinking the memory pointed to by *data from this fft algorithm should be consecutive real numbers representing the magnitude of the frequency at each point?

Depends on the FFT implementation. Most give you real-and-imaginary pairs for each frequency bin, which you have to hit with the Pythagoras stick to get absolute magnitude for.

[fromthisdrive]

The algorithm im using (shown above) says:

output: re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1)
normalized by array length

what does this mean?

[Borogove]

re(0) through re(size/2) means the first half of the output is the real components, including the DC (0) component and the Nyquist (size/2) component.

im(size/2-1) through im(1) means that the imaginary components follow the reals in reverse order (I think) from the just-below-nyquist bin down to the just-above-DC bin. (there are 2 more real bins than imaginary, because it doesn't make sense to talk about different phases at DC or at Nyquist). I assume this weird ordering made the implementation easier or faster.

Normalized by array length means that the magnitude in the bins is scaled up by the size of the array; this is because it's easier to calculate the FFT that way, and some applications of the FFT only care what the relative magnitudes are, not the absolute ones.

So to get the actual magnitude of a certain freq component, you'll have to take the real part for that bin, the imaginary part for that bin, apply pythagoras, and multiply by (1/array length).

[fromthisdrive]

cheers borogove and others, im sure there will be more questions to come though but thanks for the help so far.

Ben

[Borogove]

Since you're doing convolution, of course, you don't really need the absolute combined magnitude -- I think you just FFT the impulse response you're trying to apply, and multiply bin-by-bin between the input signal's FFT and the IR's FFT (time domain convolution = frequency domain multiplication), then IFFT the result. Or maybe you need to do complex multiplication? I've never really gone down this road.

[Z3R0T0N1N]

i thought the whole idea of fast convolution had been patented (god bless america). i could be wrong, but i was under impression that even thinking about it without paying a license fee was forbidden

[gnjp]

i thought the whole idea of fast convolution had been patented (god bless america). i could be wrong, but i was under impression that even thinking about it without paying a license fee was forbidden

There's well known patents from Gardner, Lake (dynamic) and Microsoft on various convolution techniques, dealing with stuff like calculating optimum non-uniform partitioning sizes etc for near zero latency, if you stick to uniform partition sizes i believe your ok - see bruteFIR etc, and straight fast (direct/unpartitioned) convolution is not covered by patents in general from what ive seen other wise you couldnt even implement a FIR filter.

[fromthisdrive]

multiply bin-by-bin between the input signal's FFT and the IR's FFT (time domain convolution = frequency domain multiplication), then IFFT the result.

Yeah this is true, there's a great resource for this at www.dspguide.com which has taught me most of the stuff i know for this. You can legally download a whole book for free!!

I was hoping to get the absolute magnitudes to produce some kind of graphic display of the frequency response of the mic.