Hi, first of all hello, im new here. I have been trying to develop some microphone modelling software that uses impulse responses to apply the frequency response of a mic to a sound. I'm using fast convolution and the overlap-add method. I have developed the following code to convolve a signal with an impulse of less that 1024points. I'm presently testing using white noise as my signal. The frequency respsonse of the mic looks to be partially applied to the output signal, but not well enough really and there seems to be quite a bit of distortion when i have tested it with proper audio. I think it may be a problem with the overlaps, or the length of my segments, not really sure. Should i be using some kind of windowing? I experimented using a Hamming window on the input segment, but this then also appeared in the output. Any help would be very much appreciated. The full project is available to download here: www.sycamorefarmcottages.co.uk/benkbenkbenk/hal_v2.rar #include "WavFile.h" #include "pi.h" #include "playWavFile.h" #include "fft.h" #include <iostream> using namespace std; /*This program will convolve an N point signal with a 1024 point impulse response the input signal is broken down into (N DIV 1024)+1 segments each with 1024 points. 2048point FFTs are used*/ const int FFTLENGTH = 2048; CWavFile signal, impulse; //create objects for opening wav files float fSignal[FFTLENGTH], fImpulse[FFTLENGTH], fOverlap[FFTLENGTH]; float *pSignal, *pImpulse, *pOverlap; int iNumSegs; //Number of segments input signal is broken into int iRemainingSamples;// Number of samples in last segment int iSegCount;//Counter used to keep track of which segment is being processed int iSegLength;//length of audio segments int iOverlapLength;//length of overlap void Initialise(void) { //Load the Wav Files signal.openWav("c://Temp//Signal.wav"); impulse.openWav("c://Temp//Impulse612.wav"); //initialise pointers pSignal = &fSignal[0]; pImpulse = &fImpulse[0]; //Calculate Segment Size = fftlength -ImpulseLength + 1 iSegLength = FFTLENGTH - impulse.getNumSamples() + 1; //Calculate Overlap Length = Impulse Length -1 iOverlapLength = impulse.getNumSamples() - 1; //Zero the array holding the overlapping values for(int i=0; i<FFTLENGTH; i++) { fOverlap[i]=0; } iNumSegs = (signal.getNumSamples()/(iSegLength)); iRemainingSamples = signal.getNumSamples()%(iSegLength); } void FFTImpulse(void) { //Load the impulse into the buffer and zero pad for(int i=0; i < FFTLENGTH; i++) { if(i<impulse.getNumSamples()) fImpulse[i] = impulse.getSample(i); //Load each sample into BUFFER else fImpulse[i] = 0; //Zero pad to end of buffer } //Perform FFT on Impulse Buffer realfft_split(pImpulse, FFTLENGTH); } void HammingWindow(void) { double omega = 2.0 * M_PI / (iSegLength); for (int i=0;i<iSegLength;i++) fSignal[i] = (0.54 - 0.46 * cos(omega*(i))) * fSignal[i]; } void LoadNextSegment(void)//Loads next input segment into Signal Buffer and Performs FFT { //Check to see if it is the last segment //(in which case the number of non zero samples will differ) if(iSegCount < iNumSegs) { for(int i=0; i < FFTLENGTH; i++) //if not the last segment, half fill the buffer with samples { if(i<iSegLength) fSignal[i] = signal.getSample(i+iSegLength*iSegCount); else fSignal[i] = 0; //zero pad to end of buffer } } else { for(int i=0; i < FFTLENGTH; i++) { if(i<iRemainingSamples) fSignal[i] = signal.getSample(i+iSegLength*iSegCount); //load remaining samples into input buffer else fSignal[i] = 0; //zero pad to end of buffer } } //HammingWindow(); realfft_split(pSignal,FFTLENGTH);// FFT the Input Signal Segment } void ComplexMult(void) //Perform the complex multiplication between the Signal FFT and Impulse FFT { float temp; //temporary store for real part to avoid overwriting in complex mult /* for input signal X and impulse H at frequency f: RealY[f] = RealX[f]RealH[f] - ImagX[f]ImagH[f] ImagY[f] = ImagX[f]RealH[f] + RealX[f]ImagH[f] imaginary parts of FFT are stored in reverse order starting at FFTLENGTH/2 */ for(int i=0; i<iSegLength; i++) { temp = (fSignal[i]*fImpulse[i] - fSignal[FFTLENGTH-1-i]*fImpulse[FFTLENGTH-1-i]); fSignal[FFTLENGTH-1-i] = (fSignal[i]*fImpulse[FFTLENGTH-1-i] + fSignal[FFTLENGTH-1-i]*fImpulse[i]); fSignal[i] = temp; } } void OutputSegment(void)//stores overlap and outputs the segment to the wavfile { irealfft_split(pSignal, FFTLENGTH); //Perform Inverse FFT on segment //add the last segments overlap to the this segment for(int i=0; i<iOverlapLength; i++) { fSignal[i] = fSignal[i]+fOverlap[i]; } //save the samples which will overlap the next segment for(int i=iSegLength; i<FFTLENGTH; i++) { fOverlap[i-iSegLength] = fSignal[i]; } //output segment to wav object for(int i=0; i<iSegLength; i++) { signal.setSample(fSignal[i],(i+iSegLength*iSegCount)); } iSegCount++; //increment the counter to keep track of which segment is to be processed next } void normalise(void) { float sample = signal.getSample(0); float highest = sample*sample; float hisample = 0; for(int i=0; i<signal.getNumSamples(); i++) { sample =signal.getSample(i); if(sample*sample>highest) { highest = sample*sample; hisample = sample; } } hisample=fabs(hisample); float gain=1/hisample; for(int i=0; i<signal.getNumSamples(); i++) { signal.setSample(signal.getSample(i)*gain,i); } } void main(void) { Initialise(); //initialise the variables and buffers FFTImpulse(); //FFT the impulse for(int i=0; i<iNumSegs; i++) { LoadNextSegment(); ComplexMult(); OutputSegment(); } //Output remaining samples for(int i=0; i<iRemainingSamples; i++) { signal.setSample(fSignal[i],(i+iSegLength*iSegCount)); } normalise(); signal.saveWav("c://Temp//Output.wav"); } _____________________________________ Do you know a company who employs DSP engineers? Is it already listed at http://dsprelated.com/employers.php ?

# Mic Modelling with Fast Convolution

Started by ●April 9, 2007

Reply by ●April 9, 20072007-04-09

benkbenkbenk wrote:> Hi, first of all hello, im new here. > > I have been trying to develop some microphone modelling software that uses > impulse responses to apply the frequency response of a mic to a sound. > > I'm using fast convolution and the overlap-add method. > > I have developed the following code to convolve a signal with an impulse > of less that 1024points. I'm presently testing using white noise as my > signal. The frequency respsonse of the mic looks to be partially applied > to the output signal, but not well enough really and there seems to be > quite a bit of distortion when i have tested it with proper audio. > > I think it may be a problem with the overlaps, or the length of my > segments, not really sure. Should i be using some kind of windowing? I > experimented using a Hamming window on the input segment, but this then > also appeared in the output. > > Any help would be very much appreciated. >A window is not required for fast convolution. One item that I spotted in the function ComplexMult(void) is that only iSegLength complex multiplies are performed, but it should execute FFTLENGTH complex multiplies. John

Reply by ●April 10, 20072007-04-10

> >A window is not required for fast convolution. > >One item that I spotted in the function ComplexMult(void) is that only >iSegLength complex multiplies are performed, but it should execute >FFTLENGTH complex multiplies. > >John >Thanks for looking at the code. The function ComplexMult(void) only requires iSegLength multiplies as it fills two bins with each loop. The real part is stored at the begining of the array and the imaginary is part stored in the second half of the array in reverse order. Any other ideas where these distortions may arise from? Thanks Ben _____________________________________ Do you know a company who employs DSP engineers? Is it already listed at http://dsprelated.com/employers.php ?

Reply by ●April 11, 20072007-04-11

Ben wrote:> >A window is not required for fast convolution. > > >One item that I spotted in the function ComplexMult(void) is that only > >iSegLength complex multiplies are performed, but it should execute > >FFTLENGTH complex multiplies. > > >John > > Thanks for looking at the code. > > The function ComplexMult(void) only requires iSegLength multiplies as it > fills two bins with each loop. The real part is stored at the begining of > the array and the imaginary is part stored in the second half of the array > in reverse order. > > Any other ideas where these distortions may arise from?You can try to use time-domain convolution on the blocks (this is trivial to program), using overlap-add to splice together the output. If you still have distortion, the overlap-add isn't working correctly. If not, you must look at the frequency-domain convolution engine. Regards, Andor

Reply by ●April 11, 20072007-04-11

I have now tested this with time-domain convolution on the blocks, this (although rather slow) works fine. Therefore it must be something wrong with the complex mult or maybe my FFT, but i cant really see what. The FFT im using is this, which i found on dsparchive: // 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 Many thanks Ben>You can try to use time-domain convolution on the blocks (this is >trivial to program), using overlap-add to splice together the output. >If you still have distortion, the overlap-add isn't working correctly. >If not, you must look at the frequency-domain convolution engine. > >Regards, >Andor > >_____________________________________ Do you know a company who employs DSP engineers? Is it already listed at http://dsprelated.com/employers.php ?

Reply by ●April 11, 20072007-04-11

Ben wrote:> I have now tested this with time-domain convolution on the blocks, this > (although rather slow) works fine. Therefore it must be something wrong > with the complex mult or maybe my FFT, but i cant really see what.I'm not going to debug your code. However, a quick look suggests that your complex multiplication routine does not fit with the data format of the FFT output:> // output: > // re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1)for(int i=0; i<iSegLength; i++) ... There is no im(0) in your output (it is implicitely assumed to be 0). Similarly, there is no im(size/2) (also implicitely assumed to be 0). Anycase, debugging means that you take each part of your code and test it. First test if the complex multiplication works. Then generate two complex arrays with the above format and test if the for loop works. And so on. Regards, Andor

Reply by ●April 11, 20072007-04-11

Thanks again andor for your help. Sorted it out now with this as the complexMult() function: void ComplexMult(void) //Perform the complex multiplication between the Signal FFT and Impulse FFT { realfft_split(pSignal,FFTLENGTH);// FFT the Input Signal Segment float temp; //temporary store for real part to avoid overwriting in complex mult /* for input signal X and impulse H at frequency f: RealY[f] = RealX[f]RealH[f] - ImagX[f]ImagH[f] ImagY[f] = ImagX[f]RealH[f] + RealX[f]ImagH[f] imaginary parts of FFT are stored in reverse order starting at FFTLENGTH/2 */ temp=(fSignal[0]*fImpulse[0] - fSignal[FFTLENGTH/2]*fImpulse[FFTLENGTH/2]); fSignal[FFTLENGTH/2]=(fSignal[0]*fImpulse[FFTLENGTH/2] + fSignal[FFTLENGTH/2]*fImpulse[0]); fSignal[0]=temp; for(int i=1; i<(FFTLENGTH/2); i++) { temp = (fSignal[i]*fImpulse[i] - fSignal[FFTLENGTH-i]*fImpulse[FFTLENGTH-i]); fSignal[FFTLENGTH-i] = (fSignal[i]*fImpulse[FFTLENGTH-i] + fSignal[FFTLENGTH-i]*fImpulse[i]); fSignal[i] = temp; } irealfft_split(pSignal, FFTLENGTH); //Perform Inverse FFT on segment }>I'm not going to debug your code. However, a quick look suggests that >your complex multiplication routine does not fit with the data format >of the FFT output: > >> // output: >> // re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1) > >for(int i=0; i<iSegLength; i++) > >... > >There is no im(0) in your output (it is implicitely assumed to be 0). >Similarly, there is no im(size/2) (also implicitely assumed to be 0). > >Anycase, debugging means that you take each part of your code and test >it. First test if the complex multiplication works. Then generate two >complex arrays with the above format and test if the for loop works. >And so on. > >Regards, >Andor > >_____________________________________ Do you know a company who employs DSP engineers? Is it already listed at http://dsprelated.com/employers.php ?

Reply by ●April 11, 20072007-04-11

Ben wrote:> Thanks again andor for your help. > > Sorted it out now with this as the complexMult() function:I think there is still a bug in your initial multplication (with index 0). You have to convert the implicit 0 imaginary parts of DC and Nyquist into an explicit computation. Also, using half of the FFT size for the impulse response is far from optimal. You can compute the convolution even faster if you use larger FFT sizes. I computed the optimal FFT size for a given impulse response for a certain class of DSPs here: http://groups.google.ch/group/comp.dsp/browse_frm/thread/0985ad0ca9d10309/e416169ba33926a0?&hl=de#e416169ba33926a0 On your system, the break-even points will be different, but probably of the same order of magnitude. This means that for impulse response size of 1024, you should consider 8k or 16k FFT sizes. Regards, Andor

Reply by ●April 22, 20072007-04-22

>I think there is still a bug in your initial multplication (with index >0). You have to convert the implicit 0 imaginary parts of DC and >Nyquist into an explicit computation.I'm not quite sure what you mean by this. What do i have to do to convert the implicit to the explicit? Cheers, Ben _____________________________________ Do you know a company who employs DSP engineers? Is it already listed at http://dsprelated.com/employers.php ?

Reply by ●April 22, 20072007-04-22

On 22 Apr, 12:45, "benkbenkbenk" <n0363...@hud.ac.uk> wrote:> >I think there is still a bug in your initial multplication (with index > >0). You have to convert the implicit 0 imaginary parts of DC and > >Nyquist into an explicit computation. > > I'm not quite sure what you mean by this.If you have eight real-valued samples in your data set, you will get five different coefficients in your spectrum: - The DC coefficient - The coefficients with index 1 to 3 - The Fs/2 coefficient The DC coefficient and Fs/2 coeffcient belong in different ends of your spectrum, but both are real-valued, which means that the two fit inside one complex-valued variable. So if you use a FFT implementation which returns N/2 coefficients for the real-valued DFT, the DC coefficient and Fs/2 coefficient are stored together in one complex variable. Andor said that you will need to be very careful to convert this complex-valued variable to two real-valued variables yourself. Rune