New `xcorr` convenience features

[ssfrr]

There are some commonly-used features that are useful when doing cross-correlations that we might want to add (maybe as keyword arguments to xcorr):

• demean::Bool - subtract the mean from each before cross-correlating, which would otherwise add an offset bias
• normalize::Bool - normalize by sqrt(sum(x->x^2, u)*sum(x->x^2, v)) to compensate for varying energies in signals
• maxlag::Int - specifying a maximum lag (often we care about lags shorter than our signals, and this allows less padding before the FFT, or fewer dot-products if doing it in the time-domain)
• unbiased::Bool - compensate for the varying numbers of zeros during the cross-correlation, which otherwise adds a scale bias

[ssfrr]

Specifying the lags is a little more complicated when you have signals of unequal size - for x of length N and y of length M the default range for xcorr(x, y) is -N+1 : M-1 - maybe we also allow maxlag::Tuple{Int, Int} to handle unequal cross-correlations.

[HDictus]

I think there is also merit in leaving the functions themselves straightforward and letting the user manage any pre and post processing they want to do. Unless these are best applied to the fourier transformed data, I think it might be cleaner to leave them out of the xcorr function itself (separation of concerns and all that).

It might be a good idea to add these features to DSP more genrally though, e.g. under a 'preprocessing' submodule, then they could be used for any signal processing, not just the functions into which the kwargs are built.

[galenlynch]

I don't know, having implemented all of these things for my own version of xcorr I think it would be nice to have the library handle some of them. Without accounting for the changing number of zeros, the result of xcorr is seldom useful for my research. However I am also sympathetic to having a simple function that's fast. Maybe we could use a different function, like "cross_correlation" that provides all the gloss, which but at the end of the day calls out to a fairly simple xcorr?

On Wed, May 1, 2019, 3:16 AM HDictus notifications@github.com wrote:

I think there is also merit in leaving the functions themselves straightforward and letting the user manage any pre and post processing they want to do. Unless these are best applied to the fourier transformed data, I think it might be cleaner to leave them out of the xcorr function itself (separation of concerns and all that).

It might be a good idea to add these features to DSP more genrally though, e.g. under a 'preprocessing' submodule, then they could be used for any signal processing, not just the functions into which the kwargs are built.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/JuliaDSP/DSP.jl/issues/293#issuecomment-488224982, or mute the thread https://github.com/notifications/unsubscribe-auth/ABUDZT23DVDBEQOZQS7FKBDPTE7VNANCNFSM4HJGJGOA .

[ssfrr]

I don't think this would affect the performance of xcorr, with the default options nothing it should be equivalent to what it's currently doing, and the cost of checking a few Bools should be negligible.

As for whether it makes sense to separate them, I think that might make sense in some cases but not others:

maxlag

This would affect how the cross-correlation is calculated (i.e. how much padding it needs to add before doing the FFT), so I don't think that makes sense as a separate function.

demean

this could work as a separate function, so the user would call xcorr(demean!(u), demean!(v)), though it's more verbose than xcorr(u, v, demean=true) (side note - LMK if you can think of a better name than demean which has not-great connotations as an english word)

normalize

there are a couple ways to handle this - LinearAlgebra defines normalize and normalize! so we could just have users normalize their arguments if they want. Rather than normalizing each input though, it's more efficient to apply the total normalization to whichever input is shorter, or the output if it's shorter than both the inputs (e.g. when maxlag is used). Additionally applying the normalizing to the output avoids needing to copy.

unbiased

This could potentially be a post-processing step, but would require access to all the arguments of xcorr so would be pretty verbose. It's also a step that doesn't really make sense outside the context of xcorr. Also if unbiased and normalize are both true then they can be combined into a single pass.

[galenlynch]

I call demean "centre" in my own code.

How would maxlags work? I only use it when doing direct (time-domain) cross correlation, especially with point processes. I can't think of how to do it with fft based convolution.

edit: I can't think of how to do it with FFT convolution in a way that's not just post processing, i.e. returning a smaller portion of the convolution.

[galenlynch]

Also normalize should probably imply demean and unbiased, because you'll be returning pearson correlation coefficients right?

[ssfrr]

Yeah, maybe centre would be better, though then there's the center vs. centre question.

The index twiddling to get maxlags working for FFT xcorr is a bit tricky, but the basic idea is that you only need nfft = max(sv - minlag, su + maxlag) (where sv and su are the sizes of your arguments). You still end up doing an FFT of data larger than what you eventually want so you have to truncate at the end, but you can use an FFT that's up to 50% smaller for small lags relative to your signal size. Also if the lags are limited and u and v are different sizes, you can also get some savings by truncating the longer signal, because part of it will never get used.

I've got a 1D implementation of maxlags (renamed to lagbounds) here. I get about a 45-55% speedup relative to master on a realistic load for my work (comparing 2 5s audio recordings with lags ranging from +/-1s).

using BenchmarkTools

u = rand(48000*5)
su = length(u)
v = rand(48000*5)
sv = length(v)
minlag = -48000
maxlag=48000
@assert xcorr(u,v; padmode=:none)[sv+minlag:sv+maxlag] ≈ xcorr_lin_fft(u,v; lagbounds=(minlag,maxlag))
display(@benchmark xcorr($u,$v; padmode=:none)[$(sv+minlag:sv+maxlag)])
display(@benchmark xcorr_lin_fft($u,$v; lagbounds=$(minlag,maxlag)))

I was thinking that normalization was orthogonal to centering and unbiasing. In general I'd think most people would want all three to be true, but if you know your data is centered from a previous step there's no reason for xcorr to re-center it. Also unbiasing is less relevant if you have a short signal you're searching for within a longer one.

[cwindolf]

Another enhancement would be to make xcorr n-dimensional. Right now it restricts to vectors, but I think there is nothing in the implementation which forces that. (Apart from the issue of padding -- for instance, scipy's correlate handles its equivalent of padmode=:longest by throwing an exception if neither array is larger than the other along all dimensions.) I'd be happy to implement that logic if that sounds interesting.

(Also wondering, maybe it would make sense to unify padmode in conv and xcorr as in #399 -- MATLAB, Octave, Scipy all offer padding modes called valid, same, and full for these functions, all of which are useful. This could be an alternative to the maxlag discussed above. However, I can see that there is some deprecation cycle happening to this argument in xcorr and I am not aware of the context, so maybe this is complicated.)

[cbrnr]

+1 for extending the padmode parameter to act like mode in np.correlate. Specifically, I'd find mode="same" very useful (so padmode=:same), which gives the middle samples.

[rob-luke]

Hello @cbrnr 😉 It seems this issue had stagnated and the initial author has less time for packages at the moment. Would you be able to open a PR for this?

[cbrnr]

So we meet in Julia land 😄 ! I'm not sure I have the time to tackle all enhancements mentioned in the top post, but I could at least try to implement :padmode=same. However, #403 seems to have implemented these options for conv – would this already solve most of the problem? I'm also not sure if the proposed implementation should try to get rid of the extra copy operations (e.g. for :padmode=same in xcorr I'd just return values from the middle of the array).