# Inverse of Cross Correlation?
**URL:** https://discourse.julialang.org/t/inverse-of-cross-correlation/64357
**Category:** Signal and Image Processing
**Tags:** question, dsp
**Created:** [July 9, 2021, 2:08pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357 "2021-07-09T14:08:57Z")
**Posts on this page:** 18
**Page:** 1
### Author:  [@marianoarnaiz](https://discourse.julialang.org/u/marianoarnaiz)
#### Post date: [July 9, 2021, 2:08pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/1 "2021-07-09T14:08:57Z")
Hi all.
I have a signal that is a cross correlation between to signals:
```julia
X=crosscor(A,B)
```
But in my case I get X and I can compute B… any idea how to make the inverse of the cross correlation so that I get:
```julia
A=icrosscor(X,B)
```
---
### Author:  [@JesperMartinsson](https://discourse.julialang.org/u/JesperMartinsson)
#### Post date: [July 9, 2021, 3:22pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/2 "2021-07-09T15:22:15Z")
It depends really if there is any noise components involved? Also, initial conditions might be worth considering. If signals are noise free and you can live with periodically extending them, then a frequency domain approach might be straight forward: using the DFT (i.e. using the FFT).
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [July 9, 2021, 3:33pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/3 "2021-07-09T15:33:29Z")
The relevant theory to look into is [Deconvolution - Wikipedia](https://en.wikipedia.org/wiki/Deconvolution). I don’t know what is available out of the box in Julia packages but ignore any deconvolution you find in deep learning packages; that’s something completely different.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [July 9, 2021, 3:36pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/4 "2021-07-09T15:36:44Z")
> [@GunnarFarneback](#):
>
> I don’t know what is available out of the box in Julia packages
The first thing to try would probably be a Wiener filter (= Tikhonov regularization), which is implemented e.g. in the [Deconvolution.jl package](https://github.com/JuliaDSP/Deconvolution.jl).
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [July 9, 2021, 4:33pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/5 "2021-07-09T16:33:28Z")
@marianoarnaiz, in general there is no inverse because the spectra of the 2 signals are multiplied and any spectral notches (zeros) are unrecoverable, thus the need for deconvolution. For trivial cases spectral division followed by ifft should work.
---
### Author:  [@cgeoga](https://discourse.julialang.org/u/cgeoga)
#### Post date: [July 9, 2021, 5:37pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/6 "2021-07-09T17:37:04Z")
The responses you have here seem appropriate based on your use of the word signal, but as a sanity check about an alternative reading, are you saying that your `X` is the empirical estimate for a cross correlation between vectors in the sense of \mathbb{E} \mathbf{x} \mathbf{y}^T (assuming zero means, and I recognize that this is covariance)? Because in that case your problem is much easier than the inverse problem currently being discussed, although of course if you have the correlation then you’ll also need to provide the marginal variances if you want to fully “invert” the problem, and if the means aren’t zero then you’ll need to know what was used for the empirical estimate for similar reasons.
EDIT: I didn’t realize `crosscor` was a built in function in `StatsBase`. Looking at the docs, I gather that this is not what you’re trying to do. Oops.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [July 14, 2021, 2:09pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/7 "2021-07-14T14:09:08Z")
A simple example to illustrate this problem using two time signals with different bandwidths, with and without noise added.
```julia
using DSP, FFTW
sweep(t,t0,t1,f0,f1) = (t0 < t < t1) ? sin(2π*(t-t0)*(f0 + 0.5*(f1 - f0)/(t1 - t0)*(t - t0))) : 0.0
Δt = 0.004 # sample period (s)
n = 500
t = 0:Δt:(n-1)*Δt
nl = [0.05, 0.0] # 5% noise level and noise-free
ϵ = 0.1; # white noise for deconvolution
s1 = [sweep.(t,0.5,1.5,8.,64.) .+ nl*(2*rand(n) .- 1) for nl in nl] # add random noise
s2 = [sweep.(t,0.,0.5,8.,36.) .+ nl*(2*rand(n) .- 1) for nl in nl] # add random noise
xc = [xcorr(s1,s2)[n:end] for (s1,s2) in zip(s1,s2)] # xcorr's length = 2*n - 1, cut to length n
f = LinRange(-0.5/Δt, 0.5/Δt, n+1)[n÷2+1:end]
S1 = fft.(s1); S1a = [fftshift(abs.(S1))[n÷2:end] for S1 in S1]
S2 = fft.(s2); S2a = [fftshift(abs.(S2))[n÷2:end] for S2 in S2]
Xc = fft.(xc); Xca = [fftshift(abs.(Xc))[n÷2:end] for Xc in Xc]
# deconvolve Xc for S2
S2_decon = [@. conj(Xc)*S1 / (S1*conj(S1) + ϵ^2) for (S1,Xc) in zip(S1,Xc)]
S2a_decon = [fftshift(abs.(S2_decon))[n÷2:end] for S2_decon in S2_decon]
s2_decon = [real(ifft(S2_decon)) for S2_decon in S2_decon]
# deconvolve Xc for S1
S1_decon = [@. Xc*S2 / (S2*conj(S2) + ϵ^2) for (S2,Xc) in zip(S2,Xc)]
S1a_decon = [fftshift(abs.(S1_decon))[n÷2:end] for S1_decon in S1_decon]
s1_decon = [real(ifft(S1_decon)) for S1_decon in S1_decon]
```
Input signals with and without noise, their cross-correlations and respective spectra:
Deconvolution to obtain narrower band S2 (left with noise, right without):
Deconvolution to obtain wider band S1 (left with noise, right without):
---
### Author:  [@marianoarnaiz](https://discourse.julialang.org/u/marianoarnaiz)
#### Post date: [July 14, 2021, 3:53pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/8 "2021-07-14T15:53:40Z")
This is exactly what I needed! Thanks!!!
---
### Author:  [@marianoarnaiz](https://discourse.julialang.org/u/marianoarnaiz)
#### Post date: [July 14, 2021, 5:28pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/9 "2021-07-14T17:28:07Z")
Hi @rafael.guerra would you mind sending me the plots lines too?
Thanks a lot!
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [July 14, 2021, 5:55pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/10 "2021-07-14T17:55:22Z")
@marianoarnaiz, sure, it is posted below using Plots; gr(). Need to run each plot block one at a time.
If some good soul could translate it to Makie, it would be awesome.
```julia
using Plots; gr()
default(fontfamily="Computer Modern",guidefontsize=5,tickfontsize=5,legendfontsize=5,
legend=:outertop, fg_color_legend=nothing, dpi=300)
plot(plot(t,s1, xlabel="Time(s)", label=["s1n(t)" "s1(nt)"]),
plot(t,s2, xlabel="Time(s)", label=["s2n(t)" "s2(t)"]),
plot(t,xc, xlabel="Time(s)", label=["s1n(t)*s2n(t)" "s1(t)*s2(t)"], xlims=(0,1)),
plot(f,S1a, xlabel="Frequency (Hz)",label=["|S1n(f)|" "|S1(f)|"]),
plot(f,S2a, xlabel="Frequency (Hz)",label=["|S2n(f)|" "|S2(f)|"]),
plot(f,Xca, xlabel="Frequency (Hz)",label=["|S1n(f)×S2n(f)|" "|S1(f)×S2(f)|"]),
lw=[1 0.3], lc=[:blues :black], layout=(2,3))
plot(plot(t,[s2[1], s2_decon[1]], xlabel="Time (s)", label=["s2n(t)" "s2n_decon(t)"], xlims=(0,1)),
plot(t,[s2[2], s2_decon[2]], xlabel="Time (s)", label=["s2(t)" "s2_decon(t)"], xlims=(0,1)),
plot(f,[S2a[1], S2a_decon[1]], xlabel="Frequency (Hz)", label=["|S2n(f)|" "|S2n_decon(f)|"]),
plot(f,[S2a[2], S2a_decon[2]], xlabel="Frequency (Hz)", label=["|S2(f)|" "|S2_decon(f)|"]),
lw=[1 0.3], lc=[:blues :black], layout=(2,2))
plot(plot(t,[s1[1], s1_decon[1]], xlabel="Time (s)", label=["s1n(t)" "s1n_decon(t)"]),
plot(t,[s1[2], s1_decon[2]], xlabel="Time (s)", label=["s1(t)" "s1_decon(t)"]),
plot(f,[S1a[1], S1a_decon[1]], xlabel="Frequency (Hz)", label=["|S1n(f)|" "|S1n_decon(f)|"]),
plot(f,[S1a[2], S1a_decon[2]], xlabel="Frequency (Hz)", label=["|S1(f)|" "|S1_decon(f)|"]),
lw=[1 0.3], lc=[:blues :black], layout=(2,2))
```
---
### Author:  [@marianoarnaiz](https://discourse.julialang.org/u/marianoarnaiz)
#### Post date: [July 14, 2021, 6:18pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/11 "2021-07-14T18:18:22Z")
Thanks man! This is very awesome.
It is all very similar to what I am currently trying to figure out.
My problem is about a mechanical wave with several components. I can measure one of them and use it as a filter to extract its part from the original one. It looks A LOT like the last plot. I need to study your example to see if I can apply it to my case.
Here is a picture of an example of what I am trying to do (which is NOT working hahaha)
---
### Author:  [@marianoarnaiz](https://discourse.julialang.org/u/marianoarnaiz)
#### Post date: [July 14, 2021, 6:49pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/12 "2021-07-14T18:49:59Z")
> [@rafael.guerra](#):
>
> `S1_decon = [@. Xc*S2 / (S2*conj(S2) + ϵ^2) for (S2,Xc) in zip(S2,Xc)]`
Dear Rafael. Could you explain this line:
S1\_decon = [@. Xc_S2 / (S2_conj(S2) + ϵ^2) for (S2,Xc) in zip(S2,Xc)]
The división by (S2\*conj(S2) + ϵ^2) would not just be ϵ^2 always?
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [July 14, 2021, 6:55pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/13 "2021-07-14T18:55:36Z")
The white noise term `ϵ` is small and avoids division by zero over the frequencies where S2(f) is close to zero (spectral notches). Everywhere else, the term S2.\*conj(S2) dominates.
---
### Author:  [@marianoarnaiz](https://discourse.julialang.org/u/marianoarnaiz)
#### Post date: [July 14, 2021, 7:01pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/14 "2021-07-14T19:01:42Z")
But S2\*conj(S2) wouldn’t always be 1?
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [July 14, 2021, 7:03pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/15 "2021-07-14T19:03:21Z")
@marianoarnaiz, the picture is not helpful for general audience not familiar with specialized topics. It is not my specific area of expertise but it seems that what you are trying to do is covered in this [thesis](https://etda.libraries.psu.edu/files/final_submissions/7912), which handles Rayleigh waves obtained by cross-correlation of ambient seismic noise.
---
### Author:  [@marianoarnaiz](https://discourse.julialang.org/u/marianoarnaiz)
#### Post date: [July 14, 2021, 7:08pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/16 "2021-07-14T19:08:06Z")
> [@rafael.guerra](#):
>
> thesis
this is exactly what I am doing. Jus that the explanation to the problem is not clear anywhere. Your script gave me a few ideas to try to solve my problem.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [July 14, 2021, 7:14pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/17 "2021-07-14T19:14:08Z")
The product of a complex number and its conjugate is the absolute square or square of the magnitude. The different spectra magnitudes were shown in the plots, so this is just their square.
---
### Author:  [@marianoarnaiz](https://discourse.julialang.org/u/marianoarnaiz)
#### Post date: [July 14, 2021, 7:28pm UTC](https://discourse.julialang.org/t/inverse-of-cross-correlation/64357/18 "2021-07-14T19:28:33Z")
Thanks @rafael.guerra. All was very helpfull!