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

>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 ?

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

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 ?

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

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 ?

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

>
>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 ?

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

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 ?