Adding delay calculators (delay-utils) to utils.jl

[stefano-tronci]

Hi there!

I wrote few simple methods similar to Matlab's finddelay and alignsignals. They are very simple, reflecting that I am Julia newbie. I make massive use of these methods and I though someone else might like them, so here my Pull Request. I have also other methods working with N-dimensional arrays, but a bug in xcorr prevents them to be fast and efficient. Hope you find the methods an interesting suggestion for your project.

[jfsantos]

These are interesting, could you add some tests as well?

[stefano-tronci]

Yep! I have been working on few tests, I just have to polish them up and make them a little more complete. I will commit soon.

[stefano-tronci]

Hi there, I pushed the test I was working on. It is the first time I work with tests, so let me know whether there is something else I should do.

[rob-luke]

The convention in the rest of the package is to indent the content of functions.

[rob-luke]

Nice tests and comments. Would it be possible to add Inline documentation as discussed in issue #116?

[stefano-tronci]

Yes, I never did it but looks straightforward. I will start working on it as soon as I am back from holidays early next week

[stefano-tronci]

Sorry guys, I underestimated how much of a noob I am. At first I imagine that I was supposed to add documentation directly in src/util.jl following this. Not sure is that the case tho, as after skimming through the rest of the source code I couldn't see similar lines. I am afraid I have to ask what is the correct way to add inline documentation.

[jfsantos]

We currently don't have inline docs for the library (this is the reason #116 is still open!). We currently use the documentation from the docs folder to automatically generate the documentation pages at readthedocs.

[stefano-tronci]

Hi there cool people!

It has been a while since I had feedback on this PR, so I was wondering whether I completed it correctly on my side. Maybe I did not mention it, but it is my first PR so I might have missed some step...

[galenlynch]

Hi @CrocoDuckoDucks do you have any interest in rebasing this and updating it for Julia 1.0?

[stefano-tronci]

Hi @galenlynch, sorry for the very late reply, I missed all of my notifications for a couple of weeks...

I guess I could have a go. I'll see what I can do and revert back.

[stefano-tronci]

So, I rebased on master and changed some stuff to make it work. I also changed the inline documentation using that of other functions as an example. I got a feeling that simpler tests are preferred, so I slimmed down the testing too.

[galenlynch]

Awesome, thank you! I'll review this soon.

[galenlynch]

These seem like great additions. I have made some modifications to your code to limit redundancy and to avoid changing the size of arrays.

I think your documentation should be expanded to include examples to make your sign conventions clear. In addition, some critical details, like filling in spaces after shifting with zeros, are absent from your documentation. Finally, you may want to contrast shiftsignal with circshift.

You should add these functions to the documentation in /docs/src/util.md.

I would personally prefer it if you didn't use Unicode in your code unless it really clarifies your notation, because many editing environments make it difficult to write Unicode characters.

[galenlynch]

1. Why not just use argmax instead, which returns the index of the maximum value?
2. In the case where the maximum cross-correlation is not unique, this will return the index of the first maximum. Since you are reporting lags relative to the center index, I don't think this is the behavior that you want. You would rather report the lag with the smallest absolute value, correct?

[galenlynch]

This isn't clear to me.

[galenlynch]

Your implementation seems to fill in the spaces with zeros, but that's not clear from the documentation. You might want to explicitly contrast this function with circshift in base.

[stefano-tronci]

@galenlynch Good stuff, thanks. I will address the various changes in the next few days.

[stefano-tronci]

So, I changed a few things.

1. Fixed the typos.
2. Since we prefer to mutate the first argument in mutating functions, I changed all the function signatures so that mutating and not mutating versions are used in the same way.
3. Using argmax and abs to locate the absolute peak of the cross-correlation.
4. Improved documentation and added examples and references.
5. Removed unicode characters.
6. Removed unnecessary parametric types.
7. Redundancy removed as per suggestions.

It looks like tests are now failing for Julia 0.6. Is this a problem?

[galenlynch]

@CrocoDuckoDucks Thanks for making all of those changes! Yes, we do still support Julia 0.6, but fixing those errors should be simple: argmax was added in 0.7, and replaced the function indmax. To write code that works in both 0.6 and 0.7+, you have to use Compat.jl's argmax. As Compat is already imported in the Util submodule, then if you just replace argmax with Compat.argmax it should work as expected in Julia 0.6-1.x.

[galenlynch]

Ok, this PR looks great to me now. Once the tests pass on 0.6, I will be more than happy to merge it.

Edit: finddelay seems to have one problem after the fixes: if the maximum absolute correlation is not a unique maximum, then the delay returned will still be the maximum closest to the start of the cross-correlation. Instead, I believe your intended behaviour is to report the lag closest to zero. I think the following code should achieve what you want:

function finddelay(x::AbstractVector{<: Real}, y::AbstractVector{<: Real})
    s = xcorr(y, x)
    if all(v -> v == first(s), s)
        return 0
    end

    # Find indices of maximum absolute cross-correlations
    max_corr = mapreduce(abs, max, s, init = typemin(eltype(s)))
    max_ind = findall(x -> x == max_corr, s)

    # Delay is position of peak cross-correlation relative to center.
    # If the maximum absolute cross-correlation is not unique, use the position
    # closest to the center.
    d = minimum(cld(length(s), 2) .- max_ind)
end

[galenlynch]

@CrocoDuckoDucks you should be able to fix the errors on 0.6 by using Compat.mapreduce. Sorry for suggesting something that didn't work!

Edit: ... and now that I look at the code I suggested I see that it had another error in it. To find the absolute maximum, we need to compare the aboslute value of each cross-correlation to the maximum:
max_ind = findall(x -> abs(x) == max_corr, s)

[stefano-tronci]

@galenlynch No problem, I went for a bit of trial and error. Sorry if that is making a bit of noise with notifications and the like. I am not entirely sure the behaviour of finddelay is what I want. I will do some experiments locally and report back.

[galenlynch]

No problem! Now that I think about it, I don't really understand why you are taking the absolute value of the correlation in the first place, it seems like only positive correlations would be desired. Also, there's another bug in my suggested implementation (finding the minimum lag, instead of the minimum absolute value of lag).

Here's my suggestion with all the bug fixes:

function finddelay(x::AbstractVector{<: Real}, y::AbstractVector{<: Real})
    s = xcorr(y, x)
    max_corr = maximum(s)
    max_idxs = Compat.findall(x -> x == max_corr, s)

    center_idx = cld(length(s), 2)
    # Delay is position of peak cross-correlation relative to center.
    # If the maximum cross-correlation is not unique, use the position
    # closest to the center.
    d_ind = Compat.argmin(abs.(center_idx .- max_idxs))
    d = center_idx - max_idxs[d_ind]
end

[stefano-tronci]

No problem! Now that I think about it, I don't really understand why you are taking the absolute value of the correlation in the first place, it seems like only positive lags would be desired.

Often one wants to find the delay of some kind of filter, which could have negative gain. If we pick the absolute peak of the cross-correlation we find the delay also in that case. For example:

julia> finddelay([0, 0, -1, -2, -3], [1, 2, 3])
2

So, the delay was correctly calculated also in case of phase inversion. Instead if we peak the maximum:

julia> finddelay([0, 0, -1, -2, -3], [1, 2, 3])
-1

[galenlynch]

I see, that makes sense.

[stefano-tronci]

I committed a revision I am happy with. Let me know what you think.

[galenlynch]

It seems you've changed finddelay to find the lag with maximum correlation, instead of maximum absolute correlation. From our last conversation, it sounded like you want it to be absolute correlation?

[galenlynch]

Thanks for sticking with it!