cross correlationpythontime series
I have two voltage time series from two sides of a battery with 20 seconds resolution. The signals are non-periodic and voltage2 is usually half the voltage. There are reasons to believe that there may be a slight delay in the signals and I wonder how I can estimate it.
Example of the data:
This is the very problem cross-spectral analysis is good for. Next you have an example of code using consumer prices (in differences) and price of oil, and estimating the coherency (roughly, a squared correlation coefficient broken by frequency band) and phase (lag in radians, again by frequency band).
Crudo <- dget(file="Crudo.dge")
IPC <- dget(file="ipc2001.dge")[,1]
dIPC <- diff(IPC)
datos <- ts.union(dIPC,
Crudo)
datos <- window(datos,
start=c(1979,1),
end=c(2002,1))
sp <- spectrum(datos,
main="Petróleo e IPC",
spans=rep(3,5))
par(mfrow=c(2,1))
plot(sp,plot.type="coh")
plot(sp,plot.type="phase")
These are the graphs produced by the last instructions. You can probably adapt this to your setup.
I believe that calculating the coherence would be of use to you.
The coherence at a frequency, $C_{XY}(f)$, between two series, $x(n)$ and $y(n)$, is defined as:
$C_{XY}(f) = \frac{|P_{XY}(f)|^2}{P_{X}(f)P_{Y}(f)}$
Where $P_X(f)$ and $P_Y(f)$ are the power spectrum of the series $x(n)$ and $y(n)$ respectively. $P_{XY}(f)$ is called the cross-spectrum.
Let's say that we calculate the spectrum of $x(n) = x_n$ by first Fourier transforming the windowed data:
$X(f) = \sum_{n=0}^{N-1}h_n x_n e^{-i2\pi fn}$
where $h_n = h(n)$ is the $n^{\text{th}}$ value of the window we are using. This implies that the power spectrum of $x(n)$ is:
$P_{X}(f) = |X(f)|^2$
The cross spectrum is then:
$P_{XY}(f) = |X(f) \cdot Y^{*}(f)|^{1/2}$, where $*$ denotes complex conjugate.
Thus, the coherence is:
$C_{XY}(f) = \frac{|X(f) \cdot Y^{*}(f)|}{|X(f)|^2 |Y(f)|^2}$.
After this, you could calculate the average coherence I suppose - I don't think averaging is the correct approach honestly, I will do some more checking. This is hopefully a start though.
Using a multitaper method for estimating your power spectra is a good idea (Spectrum Estimation and Harmonic Analysis, DJ Thomson, 1982 (paper) or Spectral Analysis for Physical Applications, Pervcival and Walden, 1993 (book)) as you are able to get a much reduced variance for your estimate.
If you decide to go that route, I would suggest checking out: http://www.spectraworks.com/Help/mtmtheory.html for more info.
Best Answer
To estimate the time delay between two signals you can use the cross-correlation (
np.correlate) between them and find the argmax of the cross-correlation function$$\tau_{\text{delay}} = \text{argmax }((f * g)(t)),$$
this will estimate the time offset where the signals are best aligned.
Another possible way is to use peak-detection (
scipy.signal.find_peak) and find matching peaks (e.g. with max or min in each signal or more sophisticated methods) and calculate the offset.