FFT convolution's complex multiplication problem

I've made an overlap-add/FFT convolution function, the problem with it
is that it mostly outputs samples which value is infinite. I isolated
the problem in the function and it appears to be at the level of my
complex multiplication

context : i have cut the long input signal into nblocks and put each
part into an array block[], i of course zeropadded each block, then
FFTed each block, and here i'm attempting to perform a complex
multiplication with the FFTed zeropadded short input signal :

for (i=0; i<nblocks; i++)
		for (j=0; j<blocksize; j++)
		{
			a=block[i][j];
			b=block[i][j-blocksize-1];
			c=smin[j];
			d=smin[j-blocksize-1];

			if (i<blocksize/2)
				block[i][j]=a*c - b*d;
			else
				block[i][j]=b*c + a*d;
		}

this is where i performed the complex multiplication, then i IFFT each
block and put them back together into a signal.

everything is ok when i take off the multiplication part, thus i know
that it's really that part that causes my infinite value results.

i use FFTW's real-to-real interface to perform FFTs.

anyone knows whats wrong with my complex multiplication?

in article 1123869584.368669.306040@g43g2000cwa.googlegroups.com, Michel
Rouzic at Michel0528@yahoo.fr wrote on 08/12/2005 13:59:

> I've made an overlap-add/FFT convolution function, the problem with it
> is that it mostly outputs samples which value is infinite. I isolated
> the problem in the function and it appears to be at the level of my
> complex multiplication
> 
> context : i have cut the long input signal into nblocks and put each
> part into an array block[], i of course zeropadded each block, then
> FFTed each block, and here i'm attempting to perform a complex
> multiplication with the FFTed zeropadded short input signal :
> 
> for (i=0; i<nblocks; i++)
>         for (j=0; j<blocksize; j++)
>         {
>             a=block[i][j];
>             b=block[i][j-blocksize-1];
>             c=smin[j];
>             d=smin[j-blocksize-1];
> 
>             if (i<blocksize/2)
>                 block[i][j]=a*c - b*d;
>             else
>                 block[i][j]=b*c + a*d;
>         }

i presume that "a" and "c" have your real parts and "b" and "d" are the imag
parts.  is "smin[]" your transfer function?  is it a DFT of a finite impulse
response that is at least as short as your zero padding?  (that's really
important in overlap-add.)

what about what you're storing at block[i][j-blocksize-1]?


> this is where i performed the complex multiplication, then i IFFT each
> block and put them back together into a signal.
> 
> everything is ok when i take off the multiplication part, thus i know
> that it's really that part that causes my infinite value results.

do you know for sure if your transfer function has only decent, finite
values in it?

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

robert bristow-johnson wrote:
> i presume that "a" and "c" have your real parts and "b" and "d" are the imag
> parts.

that's right.

> is "smin[]" your transfer function?

transfer function? i dont know what you mean by that, but it's the
short input signal, quite often its the filter kernel

> is it a DFT of a finite impulse
> response that is at least as short as your zero padding?  (that's really
> important in overlap-add.)

not sure i understood this one, but there is as much or less non-zero
samples (in both smin[] or block[i][]) than zero samples, and then
that's DFTed.

>
> what about what you're storing at block[i][j-blocksize-1]?

i'm storing the imaginary part of the DFTed block.


> do you know for sure if your transfer function has only decent, finite
> values in it?

if your askin for what my smin[] is made of, it's created in the main()
function by :

double *kernel;
kernelsize=3;
kernel=malloc (sizeof(double) * kernelsize);
kernel[0]=0;
kernel[1]=1;
kernel[2]=0;

and then zeropadded/realloc'ed inside the fftconvolution function to
the calculated length of 8 (thus being the FFT of (0,1,0,0,0,0,0,0))

Michel Rouzic wrote:
> 		for (j=0; j<blocksize; j++)
> 		{
> 			a=block[i][j];
> 			b=block[i][j-blocksize-1];
> 			c=smin[j];
> 			d=smin[j-blocksize-1];

Perhaps you're inadvertantly indexing a memory location outside the
array. In this loop j starts at 0, so doesn't [j-blocksize-1] result in
a value less than 0? This would cause a problem in the above lines:

> 			b=block[i][j-blocksize-1];
and
> 			d=smin[j-blocksize-1]; 

Regards,
Bob

in article 1123971289.474099.169600@g44g2000cwa.googlegroups.com, Michel
Rouzic at Michel0528@yahoo.fr wrote on 08/13/2005 18:14:

> 
> robert bristow-johnson wrote:
>> i presume that "a" and "c" have your real parts and "b" and "d" are the imag
>> parts.
> 
> that's right.
> 
>> is "smin[]" your transfer function?
> 
> transfer function? i dont know what you mean by that, but it's the
> short input signal,

i hope you mean the DFT or FFT of the short finite impulse response (the
impulse response is the filter kernel).

> quite often its the filter kernel

tomaato vs. tomahto

potaato     potahto
 

>> is it a DFT of a finite impulse
>> response that is at least as short as your zero padding?  (that's really
>> important in overlap-add.)
> 
> not sure i understood this one, but there is as much or less non-zero
> samples (in both smin[] or block[i][]) than zero samples, and then
> that's DFTed.

it's the impulse response (that would be used directly for an FIR filter if
you weren't doing "fast convolution") whose length that must not exceed the
zero-padding done in each frame (by more than one sample at least).

>> what about what you're storing at block[i][j-blocksize-1]?
> 
> i'm storing the imaginary part of the DFTed block.

your code snippet doesn't show it.

>> do you know for sure if your transfer function has only decent, finite
>> values in it?
> 
> if your askin for what my smin[] is made of, it's created in the main()
> function by :
> 
> double *kernel;
> kernelsize=3;
> kernel=malloc (sizeof(double) * kernelsize);
> kernel[0]=0;
> kernel[1]=1;
> kernel[2]=0;

looks finite to me.

> and then zeropadded/realloc'ed inside the fftconvolution function to
> the calculated length of 8 (thus being the FFT of (0,1,0,0,0,0,0,0))

a good test kernal.  you should get your input copied to the output, one
sample delayed.

in article 1123976350.825104.161260@g14g2000cwa.googlegroups.com, BobM at
BobM.DSP@gmail.com wrote on 08/13/2005 19:39:

> Perhaps you're inadvertantly indexing a memory location outside the
> array. In this loop j starts at 0, so doesn't [j-blocksize-1] result in
> a value less than 0? This would cause a problem in the above lines:

boy, i really agree with BobM.  your indexing scheme seems very
"convoluted".  (and it's this convolution you DON'T want.)


-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

oh yeah you're right, i knew it had to be a stupid error. it's
blocksize-j-1, not j-blocksize-1. never got a segfault with that tho..
maybe it's cuz i was dealing with small arrays

BobM wrote:
> Michel Rouzic wrote:
> > 		for (j=0; j<blocksize; j++)
> > 		{
> > 			a=block[i][j];
> > 			b=block[i][j-blocksize-1];
> > 			c=smin[j];
> > 			d=smin[j-blocksize-1];
>
> Perhaps you're inadvertantly indexing a memory location outside the
> array. In this loop j starts at 0, so doesn't [j-blocksize-1] result in
> a value less than 0? This would cause a problem in the above lines:
>
> > 			b=block[i][j-blocksize-1];
> and
> > 			d=smin[j-blocksize-1]; 
> 
> Regards,
> Bob

> i hope you mean the DFT or FFT of the short finite impulse response (the
> impulse response is the filter kernel).

you know, it's a convolution, doesnt matter whats the FIR kernel or
whats the input signal, to me it's just two input signals, one shorter,
one larger, and gotta be treated as such.

> > i'm storing the imaginary part of the DFTed block.
>
> your code snippet doesn't show it.

i wasnt gonna show the whole function, i did hat before and nobody
replied to my message, so i cut the code to the essential

> in article 1123976350.825104.161260@g14g2000cwa.googlegroups.com, BobM at
> BobM.DSP@gmail.com wrote on 08/13/2005 19:39:
>
> > Perhaps you're inadvertantly indexing a memory location outside the
> > array. In this loop j starts at 0, so doesn't [j-blocksize-1] result in
> > a value less than 0? This would cause a problem in the above lines:
>
> boy, i really agree with BobM.  your indexing scheme seems very
> "convoluted".  (and it's this convolution you DON'T want.)

yeah it was a stupid error, now i fixed it and it seems to work fine.
the suprising thing is that it didnt give me any segfault

> yeah it was a stupid error, now i fixed it and it seems to work fine.

no nevermind, it's giving me weird results (like convolution of a
shifted delta function with another shifted delta function gives me 8
non-zero samples). at least, it's not giving me infinite values
anymore..

in article 1123987189.813993.322130@g44g2000cwa.googlegroups.com, Michel
Rouzic at Michel0528@yahoo.fr wrote on 08/13/2005 22:39:

>> i hope you mean the DFT or FFT of the short finite impulse response (the
>> impulse response is the filter kernel).
> 
> you know, it's a convolution, doesnt matter whats the FIR kernel or
> whats the input signal, to me it's just two input signals, one shorter,
> one larger, and gotta be treated as such.

just make sure that the integer lengths of the two input signals add to a
number no more than 1 sample longer than the size of the FFT buffer.


-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

BobM wrote:
> Michel Rouzic wrote:
> > 		for (j=0; j<blocksize; j++)
> > 		{
> > 			a=block[i][j];
> > 			b=block[i][j-blocksize-1];
> > 			c=smin[j];
> > 			d=smin[j-blocksize-1];
>
> Perhaps you're inadvertantly indexing a memory location outside the
> array. In this loop j starts at 0, so doesn't [j-blocksize-1] result in
> a value less than 0? This would cause a problem in the above lines:
>
> > 			b=block[i][j-blocksize-1];
> and
> > 			d=smin[j-blocksize-1];


I also found out another mistake. once j gets bigger than blocksize/2,
b and d are out of bouds too, my mistake is making j go to blocksize
instead of blocksize/2.

I think it should be alright like that :

	for (i=0; i<nblocks; i++)
		for (j=0; j<blocksize/2; j++)
		{
			a=block[i][j];
			b=block[i][blocksize-j-1];
			c=smin[j];
			d=smin[blocksize-j-1];

			block[i][j]=a*c - b*d;
			block[i][blocksize-j-1]=b*c + a*d;
		}

but it's still giving me weird results, so there must be something
wrong still, do you see it?