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

Reply by benkbenkbenk●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 Andor●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 benkbenkbenk●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 Andor●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:

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 benkbenkbenk●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 Andor●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 benkbenkbenk●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 John Sampson●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 benkbenkbenk●April 9, 20072007-04-09

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 ?