4

I've got a distribution g that is the result of a signal convolved with itself f * f = g. I'd like to extract f using a deconvolution, but scipy.signal.deconvolve requires knowledge of the filter, which in this case is simply f. How might I go about doing this?

Ned Bingham
  • 2,749
  • 18
  • 23
  • 2
    I suspect it's not a well-defined problem: for a given `g`, there are probably multiple possible signals `f` that give rise to it. (Obviously, if `f * f = g` then `(-f) * (-f) = g` too, but I suspect that that may not be the only ambiguity.) – Mark Dickinson Feb 12 '20 at 19:44
  • 2
    Some interesting-looking theory here: https://en.wikipedia.org/wiki/Convolution_power – Mark Dickinson Feb 12 '20 at 19:48
  • Hmmm, the convolution theorem states that the fourier transform of the convolution is the point-wise product of the fourier transforms of the signals. So theoretically, I could take the fourier transform of g and then do a point-wise sqrt... Maybe? https://en.wikipedia.org/wiki/Convolution_theorem – Ned Bingham Feb 12 '20 at 19:53
  • 3
    Sounds plausible: so assuming a DFT, you'd be taking the complex square root of each coefficient, and you've got two choices for that square root for each coefficient; you'd want to align the choices so that the resulting transform is still the transform of something real. But it still seems like you'd get multiple possible `f`. Are there conditions on `f`? Does "distribution" mean "probability distribution", so that `f` has to be everywhere nonnegative? If so, I don't know how you'd go about choosing the square roots so that you get the nonnegative `f`. – Mark Dickinson Feb 12 '20 at 19:59
  • Yeah, f is a probability distribution. Non-negative and sum to 1. – Ned Bingham Feb 12 '20 at 20:01
  • 1
    Okay; I could believe that `f` is essentially uniquely determined by `g` in this case. No idea how to go about finding it, though, so I'll stop talking at this point, and upvote in the hope that someone else is attracted to answer. :-) – Mark Dickinson Feb 12 '20 at 20:01

1 Answers1

1

I was just futzing around with FFTW (version three which is the most recent and you can get is with sudo yum install fftw-devel) to try this same trick for probability zeros

look at these numbers that I copied into matlab to manipulate.

   7.111111111111112                   0   2.666666666666667                   0
  -3.190168624172971   5.606582639073907  -1.276809690904072  -2.195543579836102
  -2.606020723747883  -4.054749114715470  -1.052134133759009   1.926916438034795
   2.940900975857930  -0.483198603638057   1.720644222118320  -0.140412119317491
  -0.195012018994538   1.533362958067982  -0.821797402600774  -0.932932468036094
  -0.673564216400630  -0.170223824055756  -0.102899643706976   0.827135148011287
   0.046647507892872  -0.284192170782519   0.409049270804055  -0.347381343846294
   0.118497854845000  -0.014760807185239  -0.344899636379741   0.021398699256654
   0.027835307478515   0.044321199987332   0.200215406178045   0.110683790107336
  -0.011219166367803   0.021687271390731  -0.081234822937496  -0.133485065926827
  -0.011904309402313   0.001435596468739   0.006566970366959   0.109304320601307
  -0.004228749028931  -0.004431569979087   0.030795295176596  -0.071952062054847
   0.000442771735255  -0.003176290646744  -0.042718701210582   0.037176816672011
   0.001443162803394  -0.000877140261770   0.039572573554867  -0.011082679024575
   0.000820781607003   0.000292606798639  -0.029087458686252  -0.005029775921552
   0.000128767316409   0.000402056164143   0.016597296180403   0.012112098253013
  -0.000114815488799   0.000135337157772  -0.005597470979628  -0.012089134384515
  -0.000049480907422  -0.000025609798959  -0.001765593101901   0.007252463472756
   0.000018838838714                   0   0.004340373107703                   0

the right 2 columns are the real and imaginary components of the FFT of a beta distribution with alpha=12*(14/16) and beta=4*(14/6) I evaluated it on 17 grid points going from 0 to 1.0 on x axis, these were copied into an array that was 36 elements long (to give FFT a nice factorable with small factors number of points to work with) with the remaining elements filled with zeros. I did the forward FFT and those 2 columns are the output. The 2 left columns are the real and imaginary components of the right column times itself.

the first row is frequency 0 or a y axis shift above zero, if you sum the real components of the 2 through second to last row, and add half of the last row to that sum, THEN the sum is exactly HALF of row zero.

since row 0 should always be a positive real number you can target the signs of the rest of the rows to add up to half of it.

don't have an algorithm for that yet, but something similar to gram-schmidt (going with the largest magnitude real components first) or maybe a linear algebra approach might work. I know there are discrete optimization methods... maybe the discrete simplex method would work for this since all the signs are +/- 1 (after you normalize the last row by multiplying it by 0.5)

Edit: let S_i be the real values of the complex sqrts that we want to find signs for (i.e. it does not include the 0th and the last one which are strictly positive) we want to find signs, x_i, such that

sum(x_i*S_i)=goal

(where

goal =-0.5*(zeroth and last elements real component).

if we add the sum of the S_i (which without loss of generality we have preliminarily chosen to all be positive) then the equation becomes

sum(x_i*S_i+S_i) = goal + sum(S_i)

or

sum((x_i+1)*S_i) = goal + sum(S_i)

dividing both sides by 2 we arrive at

sum((x_i+1)/2 * S_i)] = 0.5*(goal+sum(S_i))

and the quantity

(x_i+1)/2

is a binary variable which means the problem is equivalent to the real values subset sum problem. This paper https://arxiv.org/pdf/2110.09901.pdf purports to having a good algorithm to solve the real valued subset sum problem, I only skimmed it and it does not yet make sense to me (but I am having a massive headache at the moment, my kids brought a virus home and I am sick and didn't sleep well last night)

Keith
  • 71
  • 1
  • 7
  • this is related https://stackoverflow.com/questions/66691350/algorithm-to-find-combination-of-signs-of-integers-in-a-set-such-that-the-set-su only problem is it's integers not real numbers – Keith Dec 19 '22 at 20:35