Skip to content
This repository has been archived by the owner on Sep 25, 2023. It is now read-only.

Commit

Permalink
correlation lags function (#459)
Browse files Browse the repository at this point in the history
I implemented a `cusignal.correlation_lags`  which works as expected, with notable increase in speed as expected compared to the [Scipy function]( https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlation_lags.html).

I am still working on the pytest and will update that in [`test_convolution.py`](https://github.com/rapidsai/cusignal/blob/branch-22.04/python/cusignal/test/test_convolution.py) once completed soon.

Authors:
  - Esosa E (https://github.com/sosae0)

Approvers:
  - Adam Thompson (https://github.com/awthomp)

URL: #459
  • Loading branch information
sosaeio authored Mar 10, 2022
1 parent ff66ec4 commit 40a5e76
Show file tree
Hide file tree
Showing 2 changed files with 87 additions and 1 deletion.
2 changes: 1 addition & 1 deletion python/cusignal/convolution/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,4 +19,4 @@
convolve1d2o,
convolve1d3o,
)
from cusignal.convolution.correlate import correlate, correlate2d
from cusignal.convolution.correlate import correlate, correlate2d, correlation_lags
86 changes: 86 additions & 0 deletions python/cusignal/convolution/correlate.py
Original file line number Diff line number Diff line change
Expand Up @@ -261,3 +261,89 @@ def correlate2d(
out = out[::-1, ::-1]

return out


def correlation_lags(in1_len, in2_len, mode="full"):
r"""
Calculates the lag / displacement indices array for 1D cross-correlation.
Parameters
----------
in1_size : int
First input size.
in2_size : int
Second input size.
mode : str {'full', 'valid', 'same'}, optional
A string indicating the size of the output.
See the documentation `correlate` for more information.
See Also
--------
correlate : Compute the N-dimensional cross-correlation.
Returns
-------
lags : array
Returns an array containing cross-correlation lag/displacement indices.
Indices can be indexed with the np.argmax of the correlation to return
the lag/displacement.
Notes
-----
Cross-correlation for continuous functions :math:`f` and :math:`g` is
defined as:
.. math::
\left ( f\star g \right )\left ( \tau \right )
\triangleq \int_{t_0}^{t_0 +T}
\overline{f\left ( t \right )}g\left ( t+\tau \right )dt
Where :math:`\tau` is defined as the displacement, also known as the lag.
Cross correlation for discrete functions :math:`f` and :math:`g` is
defined as:
.. math::
\left ( f\star g \right )\left [ n \right ]
\triangleq \sum_{-\infty}^{\infty}
\overline{f\left [ m \right ]}g\left [ m+n \right ]
Where :math:`n` is the lag.
Examples
--------
Cross-correlation of a signal with its time-delayed self.
>>> import cusignal
>>> import cupy as cp
>>> from cupy.random import default_rng
>>> rng = default_rng()
>>> x = rng.standard_normal(1000)
>>> y = cp.concatenate([rng.standard_normal(100), x])
>>> correlation = cusignal.correlate(x, y, mode="full")
>>> lags = cusignal.correlation_lags(x.size, y.size, mode="full")
>>> lag = lags[cp.argmax(correlation)]
"""

# calculate lag ranges in different modes of operation
if mode == "full":
# the output is the full discrete linear convolution
# of the inputs. (Default)
lags = cp.arange(-in2_len + 1, in1_len)
elif mode == "same":
# the output is the same size as `in1`, centered
# with respect to the 'full' output.
# calculate the full output
lags = cp.arange(-in2_len + 1, in1_len)
# determine the midpoint in the full output
mid = lags.size // 2
# determine lag_bound to be used with respect
# to the midpoint
lag_bound = in1_len // 2
# calculate lag ranges for even and odd scenarios
if in1_len % 2 == 0:
lags = lags[(mid - lag_bound): (mid + lag_bound)]
else:
lags = lags[(mid - lag_bound): (mid + lag_bound) + 1]
elif mode == "valid":
# the output consists only of those elements that do not
# rely on the zero-padding. In 'valid' mode, either `in1` or `in2`
# must be at least as large as the other in every dimension.

# the lag_bound will be either negative or positive
# this let's us infer how to present the lag range
lag_bound = in1_len - in2_len
if lag_bound >= 0:
lags = cp.arange(lag_bound + 1)
else:
lags = cp.arange(lag_bound, 1)
return lags

0 comments on commit 40a5e76

Please sign in to comment.