Cross-correlation by convolution (unexpected results)

Michel Rouzic wrote:

> I'm convolving sounds, not images. Those are the spectrograms of these
> signals. My convolution function apparently works, I just didn't know
> how to deal with the result. Confusing thing tho is why does image 1
> appear on the first part of image 3

Have you tried convolving the stimulus with the inverse before passing 
the stimulus through a system?  If the result isn't a one sample impulse 
then there is something wrong with your implementation.  This should be 
your acid test that you have things right.


Bob
-- 

"Things should be described as simply as possible, but no simpler."

                                              A. Einstein

Bob Cain wrote:
> Michel Rouzic wrote:
>
> >
> > I'm tempted to say that for some reason signal 1 was longer than signal
> > 2, don't know what to blame tho, I thought about blaming my ADC but
> > that would mean it plays and samples at about 44,085 Hz instead of
> > 44,100, idk if it's possible.
>
> Yes, I've seen that much deviation and more between generating and
> sampling devices.  You should play by the rule that the same clock must
> generate the stimulus that samples the response.

Oh yeah wait I did a major error in claiming that it had to both play
and sample at 44,085 Hz, cuz if so then I actually couldn't tell. The
input sampling rate has to be different than the output, otherwise I
can't see any explanation. It think it could only do it if the two had
about 30 Hz of difference

> > For you to see how my small chirp looks like, here is a zoom on it
> > http://www.geocities.com/michel0528/chirp.jpg
>
> I've never tried generating a chirp and an inverse slightly different in
> length nor doing any kind of formal analysis of what one would get but a
> small chirp instead of an impulse intuitively seems to be a good
> possibility.
>
> > You can see clearly that it does look as a small chirp. The only
> > explanation I can see is that the output signal was longer than the
> > input signal, and that thus I should do one further convolution by a
> > chirp to bring that small chirp vertical.
>
> Dunno how that could be made to work even if it could be.  I'd make sure
> the generating D/A clock was the same one as the sampling A/D clock.
> I.e. same sound card on the same computer with its clock being
> internally generated.  Working with a stimulus and response that don't
> use a common time base is pretty much out of bounds.
>
> There was a DX plugin from Sound Forge once upon a time called Acoustic
> Modeler which allowed asynchronous operation.  Somehow it was able to
> very accurately determine the relative clock rate of the response
> sampling clock and analytically generated the inverse stimulus at that
> rate before doing the convolution.  The implementor was closed mouth
> about how that rate determination was done (possibly by first generating
> a little chirp the way you have and adjusting the inverse rate based on
> its length for a second pass.)
>
>
> Bob
> --
>
> "Things should be described as simply as possible, but no simpler."
> 
>                                               A. Einstein

Bob Cain wrote:
> Michel Rouzic wrote:
>
> > I'm convolving sounds, not images. Those are the spectrograms of these
> > signals. My convolution function apparently works, I just didn't know
> > how to deal with the result. Confusing thing tho is why does image 1
> > appear on the first part of image 3
>
> Have you tried convolving the stimulus with the inverse before passing
> the stimulus through a system?  If the result isn't a one sample impulse
> then there is something wrong with your implementation.  This should be
> your acid test that you have things right.

Oh crap, you had a good idea, it revealed that something is wrong
(basically the result of the convolution of the chirp by it's inverse
gives the same chirp plus another more silent chirp starting at the
middle of the signal). Explains why the beginning of the result of my
convolution started with the first signal in it as it shouldn't even
have been there. I'm tired of having broken convolution functions, I've
went through many last summer and at last I thought I had made one that
worked (it convolved signals with windowed-sinc functions perfectly
well tho, go figure). here's the code, in case you can figure out
what's wrong. I can't see what can be tho, I must have misunderstood
something about FFT convolution

I'm not including the dftfunction(), I know it works perfectly well.
it's based on FFTW 3's real to real interface

double *zeropad(double *s, int32_t nin, int32_t nout)
{
	int32_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;
}

double *fftconv(double *s1, int32_t n1, double *s2, int32_t n2, int32_t
*n3)
{
	int32_t i, parity;
	double *s3, *t, a, b, c, d;

	*n3=n1+n2-1;	//output signal length
	parity=*n3 % 2;
	s3=malloc (*n3 * sizeof(double));

	s1=zeropad(s1, n1, *n3);
	s2=zeropad(s2, n2, *n3);
	dftfunction(s1, s1, *n3, 0, 0);
	dftfunction(s2, s2, *n3, 0, 1);

	s3[0]=s1[0]*s2[0];	//DC element
	for (i=1; i<*n3/2+parity; i++)
	{
		a=s1[i];
		b=s1[*n3-i];
		c=s2[i];
		d=s2[*n3-i];

		s3[i]=a*c - b*d;
		s3[*n3-i]=b*c + a*d;
	}
	if (parity==0)
		s3[*n3/2]=s1[*n3/2] * s2[*n3/2];

	free(s1);
	free(s2);
	dftfunction(s3, s3, *n3, 1, 1);
	return s3;
}