Hi all.
I have a signal that is a cross correlation between to signals:
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:
A=icrosscor(X,B)
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).
The relevant theory to look into is Deconvolution - Wikipedia. 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.
2 Likes
The first thing to try would probably be a Wiener filter (= Tikhonov regularization), which is implemented e.g. in the Deconvolution.jl package.
2 Likes
@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.
3 Likes
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.
1 Like
A simple example to illustrate this problem using two time signals with different bandwidths, with and without noise added.
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):
3 Likes
This is exactly what I needed! Thanks!!!
1 Like
Hi @rafael.guerra would you mind sending me the plots lines too?
Thanks a lot!
@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.
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))
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)
Dear Rafael. Could you explain this line:
S1_decon = [@. XcS2 / (S2conj(S2) + ϵ^2) for (S2,Xc) in zip(S2,Xc)]
The división by (S2*conj(S2) + ϵ^2) would not just be ϵ^2 always?
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.
But S2*conj(S2) wouldn’t always be 1?
@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, which handles Rayleigh waves obtained by cross-correlation of ambient seismic noise.
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.
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.
1 Like
Thanks @rafael.guerra. All was very helpfull!