Fwd: Potential serious problem with xcorr

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Fwd: Potential serious problem with xcorr

JuanPi
Hi Mike,

Just wanted to ask you something regarding xcorr in signal.
First of all, I get the same apparent mistake with matlab online, so
nothing to worry on that level.
However, have a look at the output of this code:
N    = 256;
t    = linspace (0, 1, N).' * 2 * pi;
dt   = t(2) - t(1);
freq = 3;
x    = sin (freq * t);
y    = x.^2;

[C lag] = xcorr (x, y);
Cw      = arrayfun (@(i)sum (x .* circshift (y, i)), lag.');

plot(lag, C, lag, Cw)

indeed the equivalence is obtained with

plot(lag, C, lag, cumtrapz(lag, Cw.*sign(lag.')))

So the C obtained with xcorr is sort of the integral of the actual
result for periodic signals.

I find this gross difference whenever I use periodic signals (both x and y).
I think the issue comes from the fact that xcorr is padding with
zeros, but here the signals are actually periodic.

Don't you think it would be nice to offer a warning about this in the
doc (which needs some fixing, will send patch), and maybe offer
xcorr_periodic or something like that for the case in which the
signals are periodic?

Regards,


--
JuanPi Carbajal
https://goo.gl/ayiJzi

-----
“An article about computational result is advertising, not
scholarship. The actual scholarship is the full software environment,
code and data, that produced  the  result.”
- Buckheit and Donoho


--
JuanPi Carbajal
https://goo.gl/ayiJzi

-----
“An article about computational result is advertising, not
scholarship. The actual scholarship is the full software environment,
code and data, that produced  the  result.”
- Buckheit and Donoho

Reply | Threaded
Open this post in threaded view
|

Re: Fwd: Potential serious problem with xcorr

Mike Miller-4
On Tue, Apr 02, 2019 at 00:42:56 +0200, JuanPi wrote:

> Just wanted to ask you something regarding xcorr in signal.
> First of all, I get the same apparent mistake with matlab online, so
> nothing to worry on that level.
> However, have a look at the output of this code:
> N    = 256;
> t    = linspace (0, 1, N).' * 2 * pi;
> dt   = t(2) - t(1);
> freq = 3;
> x    = sin (freq * t);
> y    = x.^2;
>
> [C lag] = xcorr (x, y);
> Cw      = arrayfun (@(i)sum (x .* circshift (y, i)), lag.');
>
> plot(lag, C, lag, Cw)
>
> indeed the equivalence is obtained with
>
> plot(lag, C, lag, cumtrapz(lag, Cw.*sign(lag.')))
>
> So the C obtained with xcorr is sort of the integral of the actual
> result for periodic signals.
>
> I find this gross difference whenever I use periodic signals (both x and y).
> I think the issue comes from the fact that xcorr is padding with
> zeros, but here the signals are actually periodic.

Right, the 'xcorr' function does not compute a circular correlation. I
wouldn't call this a serious problem or a mistake.

Try comparing the above against

    Cw = arrayfun (@(i)sum (x .* zeroshift (y, i)), lag.');

instead, where 'zeroshift' is defined as

    function y = zeroshift (x, i)
      if (i < 0)
        y = [x(1-i:end); zeros(-i, 1)];
      else
        y = [zeros(i, 1); x(1:end-i)];
      endif
    endfunction

> Don't you think it would be nice to offer a warning about this in the
> doc (which needs some fixing, will send patch), and maybe offer
> xcorr_periodic or something like that for the case in which the
> signals are periodic?

Sure, the docs can always be improved. It is only very briefly mentioned
currently with

    where data not provided (for example x(-1), y(N+1)) is zero.

Thanks,

--
mike

Reply | Threaded
Open this post in threaded view
|

Re: Fwd: Potential serious problem with xcorr

Mike Miller-4
In reply to this post by JuanPi
On Tue, Apr 02, 2019 at 00:42:56 +0200, JuanPi wrote:
> maybe offer
> xcorr_periodic or something like that for the case in which the
> signals are periodic?

I forgot about this last part of your message. The simplest and fastest
way to compute the circular cross-correlation today is with 'fft'.

Continuing with your example

    corr_fft = ifft (fft (x) .* conj (fft (y)));
    corr_circ = real ([corr_fft(2:end); corr_fft]);

This agrees with the circular cross-correlation you computed using
'arrayfun' and 'circshift'.

You can also use the function 'cconv' to compute the same thing

    cc = cconv (x, conj (flip (y)), length (x));
    cc = circshift (cc, 1);
    corr_circ2 = real ([cc(2:end); cc]);

I would be in favor of documenting ways to compute a circular
cross-correlation using 'fft' and/or 'cconv'.

--
mike

Reply | Threaded
Open this post in threaded view
|

Re: Fwd: Potential serious problem with xcorr

JuanPi
Hi Mike,

Thanks for the answers. So there are two possible actions here:
1. Document the circular xcorr implementation
2. Offer an xcorrcirc (or other name, I rather have it starting with
xcorr to make it TAB findable) function with the signal package

In both cases I would add the note to xcorr and mention the alternative.

What would be your choice?

Reply | Threaded
Open this post in threaded view
|

Re: Fwd: Potential serious problem with xcorr

Mike Miller-4
On Wed, Apr 03, 2019 at 04:53:58 +0200, JuanPi wrote:

> Hi Mike,
>
> Thanks for the answers. So there are two possible actions here:
> 1. Document the circular xcorr implementation
> 2. Offer an xcorrcirc (or other name, I rather have it starting with
> xcorr to make it TAB findable) function with the signal package
>
> In both cases I would add the note to xcorr and mention the alternative.
>
> What would be your choice?

I recommend option 1, document in the doc strings of 'xcorr' and/or
'cconv' how to compute a circular cross-correlation using existing
functions.

--
mike