Search code examples
rfftconvolutiondeconvolution

Convolution & Deconvolution of a vector in R

I have a convolution a*b = c, which I can compute with

a <- c(0,0,0,100,0,0,0)
b <- c(0,0,1, 2 ,1,0,0)/4

c = zapsmall(convolve(a, b[3:5], type = "f"))
#c = c(0, 25, 50, 25,  0)

Now I want to run a deconvolution assuming I only know b and c to get a. My initial idea was to use the code

a_reconstructed = deconv(c,b[3:5])
> a_reconstructed 
$q
[1]   0 100   0

$r
[1] 0 0

However, when I try this out with my actual data set, it does not seem to work. My data set is:

a = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0.00222518702742444, 0.00619744432122357, 0.0111974024392919, 
0.0169298527049305, 0.0231934447961803, 0.0298203447440921, 0.0366496261404753, 
0.0435075340291462, 0.0501839012880515, 0.0563909814495668, 0.0616617691947156, 
0.0648170048964798, 0.0648170046456146, 0.0648170046456146, 0.0648170046456146, 
0.0648170048964798, 0.0616617691947156, 0.0563909814495668, 0.0501839012880515, 
0.0435075340291462, 0.0366496261404753, 0.0298203447440921, 0.0231934447961803, 
0.0169298527049305, 0.0111974024392919, 0.00619744432122357, 
0.00222518702742444, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

b = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 0.0232558139534884, 
0.0232558139534884, 0.0232558139534884, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0)

c = zapsmall(convolve(a, b[80:122], type = "f"))
> dput(c)
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 5.1749e-05, 0.000195875, 0.00045628, 0.000849997, 
0.00138938, 0.002082876, 0.002935193, 0.003946996, 0.005114064, 
0.006425482, 0.007859476, 0.009366849, 0.010874221, 0.012381593, 
0.013888965, 0.015396337, 0.016830332, 0.01814175, 0.019308818, 
0.020320621, 0.021172938, 0.021866434, 0.022405817, 0.022799534, 
0.023059939, 0.023204065, 0.023255814, 0.023255814, 0.023255814, 
0.023255814, 0.023255814, 0.023255814, 0.023255814, 0.023255814, 
0.023255814, 0.023255814, 0.023255814, 0.023255814, 0.023255814, 
0.023255814, 0.023255814, 0.023255814, 0.023255814, 0.023204065, 
0.023059939, 0.022799534, 0.022405817, 0.021866434, 0.021172938, 
0.020320621, 0.019308818, 0.01814175, 0.016830332, 0.015396337, 
0.013888965, 0.012381593, 0.010874221, 0.009366849, 0.007859476, 
0.006425482, 0.005114064, 0.003946996, 0.002935193, 0.002082876, 
0.00138938, 0.000849997, 0.00045628, 0.000195875, 5.1749e-05, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0)

a_reconstructed = deconv(c,b[80:122])

> dput(a_reconstructed$q)
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0.002225207, 0.006197418, 0.011197415, 0.016929831, 
0.023193469, 0.029820328, 0.036649631, 0.043507529, 0.0501839239999999, 
0.0563909740000001, 0.061661742, 0.064817039, 0.064816996, 0.0648169960000001, 
0.0648169959999999, 0.064816996, 0.061661785, 0.056390974, 0.0501839239999997, 
0.0435075290000002, 0.036649631, 0.0298203279999999, 0.0231934689999999, 
0.0169298310000002, 0.0111974149999999, 0.00619741800000028, 
0.00222520699999964, 1.13637940203637e-16, 1.38405183581353e-16, 
-5.01172454231425e-17, 2.62241400469931e-18, -2.03965533698836e-17, 
4.95344867554326e-18, -9.70293181738747e-17, 7.28448334638698e-17, 
4.02103480720561e-17, -2.30189673745828e-17, -1.04896560187973e-17, 
-1.82986221661241e-16, 2.68943125148608e-16, -1.74244841645577e-16, 
1.64920702962201e-16, -1.74827600313301e-18, -1.2179656155159e-16, 
3.41787958612477e-16, -2.70982780485596e-17, -2.18825879725465e-16, 
2.29898294411973e-16, -3.84620720689233e-17, -2.88465540516926e-17, 
8.74138001566417e-18, 1.54139667609549e-16, -2.11250017045223e-16, 
-4.29999998611063e-08, 4.29999999272494e-08, 1.85608644558132e-16, 
-8.24603514811005e-17, -2.08918973482187e-16, -4.29999999444408e-08, 
4.29999999488115e-08, -1.21213820668263e-16, -4.52220690574938e-16, 
2.77393116938224e-16, 8.62482650368385e-17, -1.42775873589185e-16, 
-7.25534363456311e-17, 2.03674136580597e-16, -7.22620659039672e-17, 
1.69291366293459e-16, -2.94793802041224e-16, 1.15396468903667e-16, 
1.02309465968249e-16)

> dput(a_reconstructed$r)
c(3.89630606243601e-19, 1.31393359644378e-18, -1.63829696401998e-18, 
-5.61885178543501e-20, -2.29098403315009e-18, 4.76327819835599e-19, 
8.53031661056526e-19, 8.52647435662004e-19, 1.61948977202266e-18, 
-3.57759068408242e-18, 2.31285675649752e-18, 2.96531841891069e-19, 
1.98260965322424e-18, 5.56346377161829e-19, -1.29645340894417e-18, 
5.97947502875264e-18, 2.06744240171553e-18, 1.67523099196064e-19, 
4.97558981899935e-18, 5.36951597416393e-18, 3.51852316410607e-18, 
2.86339963448168e-18, 7.47172270907357e-18, -1.75040516602959e-19, 
-9.99999993440607e-10, 3.78214748648309e-18, 8.77083178137432e-18, 
8.40165156141487e-18, 3.92122239684741e-18, -9.99999997263831e-10, 
2.13617403272263e-18, -9.01974613145322e-19, -1.38893610330793e-17, 
-5.77150302229999e-18, -4.17818084768659e-18, -5.45246783355247e-18, 
-8.09949340619467e-18, -3.41702793469095e-18, -5.09771382989791e-18, 
-5.20268358577031e-19, -8.74999682297592e-18, -5.24133381038652e-18
)

Here is a visualiuation of the two vectors a (=f) and b (=g) and the resulting convolution ab (=fg). How can I get the deconvolution working for my dataset?

enter image description here

Solution

  • Reversing convolve:

    bb <- b[80:122]
n1 <- length(bb) - 1
n <- length(a)
aa <- zapsmall(Re(fft(fft(c(numeric(n1), c, numeric(n1)))/Conj(fft(c(bb, numeric(n - 1)))), TRUE))[-1:-n1]/(n + n1), 6)
all.equal(a, aa, tolerance = 1e-6)
#> [1] TRUE

    Note that this loses some precision, hence the tolerance = 1e-6.