# Using ListDeconvolve to find velocity distrbution

Consider a Distribution of Particles correspondent with the NormalDistribution:

PDF[NormalDistribution[], x]

The Particles also have a velocity distribution:

Exp[-x] UnitStep[x]

Which looks like this

So all particles move to the right. After a certain timestep, the particles will therefore have another spatial distribution. The new spatial distribution can be obtained by the convolution of the old distribution (which was the first picture) with the associated velocity distribution (associated with the timestep). For a certain timestep we get:

Convolve[PDF[NormalDistribution[], x], Exp[-x] UnitStep[x]]

For testing purposes I want to use the initial and the final spatial distribution, to find the velocity distribution.

What I get, is the following:

cv[y_]=Convolve[PDF[NormalDistribution[], x], Exp[-x] UnitStep[x], x, y];
list1 = Table[PDF[NormalDistribution[], x], {x, -10, 10, 0.01}];
list2 = Table[cv[y], {y, -10, 10, 0.01}];
list3 = ListDeconvolve[list1, list2];
ListPlot[LowpassFilter[list3, 1000], PlotRange -> All]

I kept the "lowpassfilter" because of lazyness, but set cutoff to 1000, to make it idle. Decreasing the cutoff frequency of the low pass also does not help, because signal energy is distributed all over the place. That's why the values are so small. How can I reconstruct the velocity distribution? By the way: IF I would take the Exponential function as my initial spatial distribution, and the Normaldistribution as my velocity distribution instead, then ListDeconvolve would give a good result, this is why I did not totally discarded ListDeconvolve.

EDIT:

And this is the data I am really face with:

The blue curve is the initial spatial dist, and the orange is that of the next timestep. As you can see, it is a little wider. If I repeat the steps from above, the velocity distribution claims to be:

ListPlot[LowpassFilter[list3, 0.001], PlotRange -> All]

This is how it looks like after applying a lowpass filter with very low cutoff. Without lowpass it's a mess.

• shouldn't you be looking at the smoothhistogram of list3 instead of plotting the list itself? I expect something symmetric from how your data spreads in the timestep you showed. – tsuresuregusa Feb 28 '16 at 21:54
• Save yourself huge time: Replace "cv[y_] :=" with "cv[y_] =". The convolution is exactly computable, so need not be repeatedly (numerically) evaluated (since you are feeding decimals through list2's Table[]). This eliminates the numerical ringing in list2 over the first 300-ish list entries at the $10^{-12}$ level. – Eric Towers Feb 28 '16 at 21:59
• SmoothHistogram does not give the desired result. I slipped up before: there shouldn't be negative numbers in the distribution. – Anton Alice Feb 28 '16 at 22:07
• your last plot is not a distribution of velocities, unless your mean velocity is 400. You need to have negative velocities in the distribution since your particles are moving both to left and right. – tsuresuregusa Feb 29 '16 at 0:00
• the last plot has nothing to do with the first few plots. and I did not say, that it is the actual velocity dist. It is the ASSOCIATEd velocity distribution. for the actual velocity in meters/s have to convert the values on the x-axis according to the time-step (which is measured in seconds). the mean velocity is however zero. – Anton Alice Feb 29 '16 at 0:34

The ringing you are seeing is caused by the sharp jump in the derivative of the velocity distribution as it is reflected at $10$ to $-10$ in list2. (ListDeconvolve[] reflects the list at the boundaries. If you instead pad with zeroes, you'll get different ringing, but it won't really improve things.) If you will generate your data as

cv[y_] = Convolve[PDF[NormalDistribution[], x], Exp[-x] UnitStep[x],
x, y];
list1 = Table[PDF[NormalDistribution[], x], {x, -100, 100, 0.01}];
list2 = Table[cv[y], {y, -100, 100, 0.01}];
list3 = ListDeconvolve[list1, list2];
Print[ListPlot[list1, PlotRange -> {-10^-21, 10^-21}]];
Print[ListPlot[list2, PlotRange -> {-10^-11, 10^-11}]];
Print[ListPlot[list3, PlotRange -> All]];
Print[ListPlot[LowpassFilter[list3, 1000], PlotRange -> All]];

You'll find three things:

• You have a clear peak at $10\,000$, as expected.
• You have some ringing -- again, this is a boundary effect.
• Your low-pass filter isn't doing very much.

If you put list1 back to the interval $[-10,10]$ (and swap the order of the lists in ListDeconvolve[]) your peak is there and crazy ringing is back.

If instead you put list2 back to the interval $[-10,10]$ (and leave the arguments to ListDeconvolve[] as you provided them) your peak is there and the crazy ringing has been crushed a bit.

This tells you something you already know: the normal distribution tapers to horizontal faster than the exponential distribution does.

This is happening because your position and velocity distribution functions have infinite bandwidth. The ringing you are seeing is normally "offset" by all the Fourier bins that have been chopped off.

• Never forget: the derivative is multiplied by the frequency. So a bin $\sin kt$ is $k \sin kt$ in the derivatives. High frequency truncated bins contain rather more energy than one would first expect. In the above, you see the energy in the ringing that comes and goes based on how early or late the truncation is. Other than regularization using a serious nonlinear deconvolution technique, I have never found any simple method that worked. This text is not a bad start. (I've used the 1996 edition.) – Eric Towers Mar 1 '16 at 4:30