FFT convolution's complex multiplication problem

It's difficult what to say what is wrong without knowing more
information about how your function works. What is the total input size
in samples? The filter size? The input block size?

Michel Rouzic wrote:
> When I convolve with a function that is an array of 20,000 samples set
> to 0, and the element 4,000 being set to 1, my output signal looks like
> the orginal input (with a longer length of course), unmodified, except
> that on every section of 32,768 samples, the first 4,000 are set to
> zero.

I'm assuming your filter size is 20,000 taps. If my math is correct,
the max your input block size can be is 12,769 samples (32768-20000+1).
Anyway, the result you're getting seems correct to me. As long as your
input signal is shifted in time by 4000 samples, which it sounds like
it is, then this is the correct behavior. Also, it seems your input
block isn't too big, because I'd expect to see some non-zero samples
folding back into the first 4000 samples if it was too big. Sounds like
the convolution itself is working.

> I also tried with a windowed-sinc function, i used 88.2 Hz as a cutoff
> frequency, 8.82 Hz as the bandwidth of the rolloff, and when convolving
> with the input signal, not only the output signal has regions regularly
> set to zero, but the output signal doesnt look/sound at all like it's
> been low-passed filtered

How many taps is this filter? If it's not many taps, then this may be
such a shallow filter that it there's almost no effect. What kind of
window is being used to generate the filter? According to the first
result you mentioned it sounds like your convolution is working, so
this may be related to the filter you're using. Does this same filter
work with a time-domain implementation of convolution?

> any ideas on what can be wrong?

Again, please provide some more information. However, my hunch is that
the convolution itself is working. But confirm that your input block
isn't too big. Also check the way you're piecing the output blocks back
together. Whether you're doing overlap-add or overlap-save, check that
the overlapped portions are being added together before sending them
off to output.

BobM wrote:
> It's difficult what to say what is wrong without knowing more
> information about how your function works. What is the total input size
> in samples? The filter size? The input block size?

I used two test input signals, one 2,103,154 samples long, the other,
972 samples long. as for the filter, it's 20,000 samples long in the
case of the shifted delta function. I'm not sure to know what you mean
by input block size

> I'm assuming your filter size is 20,000 taps. If my math is correct,
> the max your input block size can be is 12,769 samples (32768-20000+1).

huh? ok, we maybe got the finger on something wrong... what I do is, in
the case of the convolution by the shifted delta function, take the
20,000 samples long signal, and provided that the other input signal is
longer than that, zero-pad that signal to the power of two above twice
the length of the signal (I guess that is 65,536 here), then, I cut the
2,103,154 samples long signal into blocks of 32,768 samples and I
zero-pad each block to 65,536 (i think i made no mistake)

> Anyway, the result you're getting seems correct to me. As long as your
> input signal is shifted in time by 4000 samples, which it sounds like
> it is, then this is the correct behavior.

No, it's not getting shifted, but muted, I mean, some regions of 4000
samples are getting, each 32,768 samples, set to zero, as the other
28,768 samples are unchanged, even when using a windowed-sinc function.
This is indeed quite far from a normal behavior

> Also, it seems your input
> block isn't too big, because I'd expect to see some non-zero samples
> folding back into the first 4000 samples if it was too big. Sounds like
> the convolution itself is working.
>
> > I also tried with a windowed-sinc function, i used 88.2 Hz as a cutoff
> > frequency, 8.82 Hz as the bandwidth of the rolloff, and when convolving
> > with the input signal, not only the output signal has regions regularly
> > set to zero, but the output signal doesnt look/sound at all like it's
> > been low-passed filtered
>
> How many taps is this filter?

20,001. I calculate the length of the filter by 4/BW (BW=the bandwidth
of the roll off) and make sure the number is odd.

> If it's not many taps, then this may be
> such a shallow filter that it there's almost no effect.
> What kind of
> window is being used to generate the filter?

Blackman. I made sure it's working right. I even made sure that the
first and last bins were always non-zero.

> According to the first
> result you mentioned it sounds like your convolution is working, so
> this may be related to the filter you're using. Does this same filter
> work with a time-domain implementation of convolution?

Haven't made a time-domain convolution. i thought i'd go straight to
the goal... anyways, it looks perfectly right.

>
> > any ideas on what can be wrong?
>
> Again, please provide some more information. However, my hunch is that
> the convolution itself is working. But confirm that your input block
> isn't too big. Also check the way you're piecing the output blocks back
> together. Whether you're doing overlap-add or overlap-save, check that
> the overlapped portions are being added together before sending them
> off to output.

Maybe it's not very appropriate to paste such a long code, but I have
the feeling that it'll go much faster if you can get to see a
stripped-down-to-the-essential version of my function, so here it is
(the first functions are here because they are used in the fftconv
function)


#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>

double log_2(double x)	//log base 2 performing function
{
	double y;
	y=log(x)/log(2);
	return y;
}

double modfloat(double a, double b)	//modulo returning function
{
	double c, d, x;

	c=(uint32_t) a/b;
	d=a/b;
	x=d-c;
	return x;
}

double roundup(double in)	//function rounding the input to the upper
integer (unless the input already is an integer)
{
	double out;
	if (modfloat(in,1) == 0)
		out=in;
	else
		out=in - (modfloat(in,1)) + 1;
	return out;
}

double *zeropad(double *s, uint32_t nin, uint32_t nout)
{
	uint32_t i;
	double *h;

	h=malloc (sizeof(double) * nout);
	for (i=0; i<nin; i++)
		h[i]=s[i];
	for (i=nin; i<nout; i++)
		h[i]=0;
	return h;
}

void dftfunction(double *in, uint32_t N, uint8_t method, uint8_t root)
{
	uint32_t i;
	fftw_plan p;

	p = fftw_plan_r2r_1d(N, in, in, method, FFTW_ESTIMATE);
	fftw_execute(p);
	fftw_destroy_plan(p);

	if (root==1)
		for (i=0;i<N;i++)
			in[i]=in[i]/sqrt(N);
}

double *fftconv(double *s1, uint32_t n1, double *s2, uint32_t n2,
uint32_t *n3)
//s1 and s2 the two input signals, n1 and n2 their respective length
{
	uint32_t i, j, nmin, nmax, blocksize, nblocks;
	double *s3, *smin, *smax, **block, a, b, c, d;

	printf("fftconv...\n");

	*n3=n1+n2-1;	//output signal length
	s3=malloc (*n3 * sizeof(double));
	if (n1<n2)
	{
		smin=s1;
		smax=s2;
		nmin=n1;
		nmax=n2;
	}
	else
	{
		smin=s2;
		smax=s1;
		nmin=n2;
		nmax=n1;
	}

	blocksize=pow(2, (uint32_t) roundup(log_2((double) (nmin)*2)));	//size
of blocks for the overlap-add of the long signal, size of the
zeropadded short signal
	smin=zeropad(smin, nmin, blocksize);
	dftfunction(smin, blocksize, 0, 0);	//FFT of the shortest signal

	nblocks=*n3/(blocksize/2);	//number of blocks for the o-a of the long
signal
	block=malloc (nblocks * sizeof(double *));
	for (i=0; i<nblocks; i++)	//loop filling the contents of the long
signal in an array of blocks
	{
		block[i]=malloc (blocksize * sizeof(double));
		for (j=0; j<blocksize/2; j++)
			block[i][j]=smax[i*(blocksize/2) + j];
		for (j=blocksize/2; j<blocksize; j++)
			block[i][j]=0;
		dftfunction(block[i], blocksize, 0, 1);	//FFT of each block
	}

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

			block[i][j]=a*c - b*d;
			block[i][blocksize-j]=b*c + a*d;
		}
		a=block[i][blocksize/2];  // Nyquist
		c=smin[blocksize/2];
		block[i][blocksize/2]=a*c;
		dftfunction(block[i], blocksize, 1, 1);	//IFFT of the result of the
multiplication in each block
	}

	for (i=0; i<*n3; i++)	//initialisation of the output signal s3
		s3[i]=0;

	for (i=0; i<nblocks; i++)	//filling of s3 from the content of the
blocks
		for (j=0; j<blocksize/2; j++)
			s3[i*(blocksize/2) + j]=s3[i*(blocksize/2) + j]+block[i][j];
	return s3;
	free(block);
}

Michel Rouzic wrote:
> 	for (i=0; i<nblocks; i++)	//filling of s3 from the content of the
> blocks
> 		for (j=0; j<blocksize/2; j++)
> 			s3[i*(blocksize/2) + j]=s3[i*(blocksize/2) + j]+block[i][j];

The problem looks to be here. It looks like only the first half of each
block is being sent to output, and that there's no overlap.

Be sure you are overlapping each block and adding the overlapped parts.
If your input block size (the amount of input samples you filter each
time) is always blocksize/2, then an overlap size of blocksize/2
samples is ok.

However, while it looks like you intend to overlap and add the blocks,
the above code doesn't do that. You'll need to output the entire block,
overlapped by blocksize/2 samples (assuming the input block size is
blocksize/2).

I still recommend getting familiar with the time-domain implementation
of an FIR. It's more intuitive for understanding time-domain effects
like group delay.

BobM wrote:

OK, i turned this part into this :

	for (i=0; i<*n3-(blocksize/2); i++)	//filling of s3 from the content
of the blocks
		s3[i]=block[i/(blocksize/2)][i%(blocksize/2)];	//first half of each
block
	for (i=blocksize/2; i<*n3; i++)
		s3[i]+=block[i/(blocksize/2) - 1][i%(blocksize/2) +
(blocksize/2)];	//second half of each block

Now I dont have the problems I had before, and convolving the shifted
delta function gives the expected result.

As for the windowed-sinc function, forget anything I could say about
it, the code looks right, but all it outputs is a shifted delta
function... I tried a convolution with a little kernel (2 bins,
first=1, second = -1) and the result sounds right, so I can say my
convolution works right.

thank you for help

Michel Rouzic wrote:
> Now I dont have the problems I had before, and convolving the shifted
> delta function gives the expected result.

Good to hear it's working now.

> As for the windowed-sinc function, forget anything I could say about
> it, the code looks right, but all it outputs is a shifted delta
> function...

Read up on group delay. It's expected that there will be a time shift
(delay) with any FIR filter. If it doesn't seem filtered from the
frequency domain point of view, it seems like there may be a problem
with the filter coefficients.