Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New xcorr convenience features #293

Open
ssfrr opened this issue Apr 29, 2019 · 11 comments
Open

New xcorr convenience features #293

ssfrr opened this issue Apr 29, 2019 · 11 comments

Comments

@ssfrr
Copy link
Contributor

ssfrr commented Apr 29, 2019

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
Copy link
Contributor Author

ssfrr commented Apr 29, 2019

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
Copy link
Contributor

HDictus commented May 1, 2019

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
Copy link
Member

galenlynch commented May 1, 2019 via email

@ssfrr
Copy link
Contributor Author

ssfrr commented May 1, 2019

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
Copy link
Member

galenlynch commented May 1, 2019

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
Copy link
Member

galenlynch commented May 1, 2019

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

@ssfrr
Copy link
Contributor Author

ssfrr commented May 2, 2019

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
Copy link

cwindolf commented May 22, 2021

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
Copy link

cbrnr commented Jul 28, 2021

+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
Copy link
Member

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
Copy link

cbrnr commented Jul 30, 2021

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).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

6 participants