Calculating the PDF of a waveform from its samples A while ago I was trying different ways to draw digital waveforms, and one of the things I tried was, instead of the standard silhouette of the amplitude envelope, to display it more like an oscilloscope.  This is what a sine and square wave look like on a scope:



The naïve way to do this is:


Divide up the audio file into one chunk per horizontal pixel in the output image
Calculate the histogram of sample amplitudes for each chunk
Plot the histogram by brightness as a column of pixels


It produces something like this:


This works fine if there are a lot of samples per chunk and the signal's frequency is unrelated to the sampling frequency, but not otherwise.   If the signal frequency is an exact submultiple of the sampling frequency, for instance, the samples will always occur at exactly the same amplitudes in each cycle and the histogram will just be a few points, even though the actual reconstructed signal exists between these points.  This sine pulse should be as smooth as the above left, but it isn't because it's exactly 1 kHz and the samples always occur around the same points:



I tried upsampling to increase the number of points, but it doesn't solve the issue, just helps smooth things out in some cases.

So what I'd really like is a way to calculate the true PDF (probability vs amplitude) of the continuous reconstructed signal from its digital samples (amplitude vs time).  I don't know what algorithm to use for this.  In general, the PDF of a function is the derivative of its inverse function.

PDF of sin(x): $\frac{d}{dx} \arcsin x = \frac{1}{\sqrt{1-x^2}}$

But I don't know how to calculate this for waves where the inverse is a multi-valued function, or how to do it fast.  Break it up into branches and calculate the inverse of each, take the derivatives, and sum them all together?  But that's pretty complicated and there's probably a simpler way.

This "PDF of interpolated data" is also applicable to an attempt I made to do kernel density estimation of a GPS track.  It should have been ring shaped, but because it was only looking at the samples and not considering the interpolated points between the samples, the KDE looked more like a hump than a ring.  If the samples are all we know, then this is the best we can do.  But the samples are not all we know.  We also know that there's a path between the samples.  For GPS, there's no perfect Nyquist reconstruction like there is for bandlimited audio, but the basic idea still applies, with some guesswork in the interpolation function.
What I'd go with is essentially Jason R's "random resampler", which in turn is a presampled-signal based implementation of yoda's stochastic sampling.

I've used simple cubic interpolation to one random point between each two samples. For a primitive synth sound (decaying from a saturated non-bandlimited square-like signal +even harmonics to a sine) it looks like this:



Let's compare it to a higher-sampled version,



and the weird one with the same samplerate but no interpolation.



Notable artifact of this method is the overshoot in the square-like domain, but this is actually what the PDF of the sinc-filtered signal (as I said, my signal is not bandlimited) would also look like and represents the perceived loudness much better than the peaks, if this were an audio signal.

Code (Haskell):

cubInterpolate vll vl v vr vrr vrrr x
    = v*lSpline x + vr*rSpline x
      + ((vr-vl) - (vrr-vll)/4)*ldSpline x
      + ((vrr-v) - (vrrr-vl)/4)*rdSpline x
     where lSpline x = rSpline (1-x)
           rSpline x = x*x * (3-2*x)
           ldSpline x = x * (1 + x*(x-2))
           rdSpline x = -ldSpline (1-x)

                   --  rand list   IN samples  OUT samples
stochasticAntiAlias :: [Double] -&gt; [Double] -&gt; [Double]
stochasticAntiAlias rs (lsll:lsl:lsc:lsr:lsrr:[]) = []
stochasticAntiAlias (r:rLst) (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)
    = ( cubInterpolate lsll lsl lsc lsr lsrr lsrrr r )
          : stochasticAntiAlias rLst (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)


rand list is a list of random variables in range [0,1].