Calculating the PDF of a waveform from its samples

Question

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:

1. Divide up the audio file into one chunk per horizontal pixel in the output image
2. Calculate the histogram of sample amplitudes for each chunk
3. 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.

Answer 1

Interpolate to several times the original rate (e.g. 8x oversampled). This allows you to assume a piecewise linear signal. This signal will have very little error compared to the infinite resolution, continuous sin(x)/x interpolation of the waveform.

Assume every pair of oversampled values has a continuous line from one value to the next. Use all the values between. This gives you one thin horizontal slice from y1 to y2 to be accumulated into an arbitrary resolution PDF. Each rectangular slice of probability must be scaled to a 1/nsamples area.

Using the line between samples rather than the sample themselves prevents a "spikey" PDF, even in the case when there is a fundamental relationship between the sampling period and the waveform.

Answer 2

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] -> [Double] -> [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].

Answer 3

While your approach is theoretically correct (and needs to be slightly modified for non-monotonic functions), it is extremely hard to calculate the inverse of a generic function. As you say you'll have to deal with branch points and branch cuts, which is doable, but you seriously wouldn't want to.

As you already mentioned, regular sampling samples the same set of points and as such is highly susceptible to poor estimates in regions where it doesn't sample (even if the Nyquist criterion is satisfied). In this case, sampling for a longer period doesn't help either.

In general, when dealing with probability density functions and histograms, its a much better idea to think in terms of stochastic sampling than regular sampling (see the linked answer for an introduction). By sampling stochastically, you can ensure that every point has an equal probability of being "hit" and is a much better way to estimate the pdf.

Here's an example: Consider the function $f(x)=\sin(20\pi x) + \sin(100\pi x)$. Now if I sample this with a sampling frequency $f_s=1000$ Hz (Nyquist frequency, $f_N=100$ Hz), the resulting probability density is the plot on the left (401 equispaced bins between -2 and 2). It doesn't matter if I sample for 10 seconds or 100. It still remains the same. On the other hand, sampling stochastically at a rate of $1000$ samples (uniform distribution) per second (I'm not using Hz here, because that implies a different meaning) for 30 seconds gives the plot on the right (same binning).

You can easily see that although it is noisy, it is a much better approximation to the actual PDF than the one on the right which shows zeros in several intervals and large errors in several others. By having a longer observation time, you can bring down the variance in the one on the right, eventually converging to the exact PDF (dashed black line) in the limit of large observations.

Answer 4

Kernel Density Estimation

One way to estimate the PDF of a waveform is to use a kernel density estimator.

This takes all your samples $x(n)$ and a kernel function (e.g. a Gaussian), $K(\mathrm{x})$ and convolves them with $\delta(\mathrm{x} - x(n))$ (the Dirac delta) to give an estimate of the PDF, $\hat{P}$ :

$$\hat{P}(\mathrm{x}) = \sum_{n=0}^N K( \mathrm{x} - x(n) ) $$

---

Update: Interesting additional information.

Suppose you have your signal $x(n)$ for $n= 0,1,...,N-1$, then --- as you say --- you can also know its DFT $X(k)$:

$$ X(k) = \sum_{n=0}^{N-1} x(n) e^{-\jmath 2\pi n k /N} $$

So that $X(k)$ is the coefficient of $e^{\jmath 2\pi n k /N}$:

$$x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X(k) e^{\jmath 2\pi n k /N}$$

So a guess at what you're aftermight be to convolve all the PDFs of each Fourier component together:

$$ |X(k)| \frac{1}{\sqrt{1-x^2}} $$

However that does not take account of the way the phase of $X(k)$ contributes (or not) to the addition in $x(n)$.

More thought required, though!

Answer 5

As you indicated in one of your comments, it would be attractive to be able to calculate the histogram of the reconstructed signal using only the samples and the PDF of the sinc function that interpolates bandlimited signals. Unfortunately, I don't think this is possible because the sinc's histogram doesn't have all of the information that the signal itself has; all information on the time-domain positions where each value is encountered is lost. This makes it impossible to model how the scaled and time-delayed versions of the sinc would sum together, which is what you would want in order to calculate the histogram of the "continuous" or up-sampled version of the signal without actually doing the up-sampling.

I think you're left with interpolation as the best option. You did indicate a couple of issues that prevented you from wanting to do this, which I think can be addressed:

• Computational expense: This is of course always a relative concern, depending upon the specific application that you want to use this for. Based on the link that you posted to the gallery of renderings that you've collected, I'm assuming you're wanting to do this for visualization of audio signals. Whether you're interested in this for a real-time or offline application, I would encourage you to prototype an efficient interpolator and see if it is truly too costly. Polyphase resampling is a good way to do this that is flexible (you can use any rational factor).

• Bias to periodic components related rationally to the sampling frequency: While you can't eliminate this completely, you should be able to mitigate it somewhat by interpolating by a "strange" factor: instead of up-sampling by 4, try 71/18 (just an example). This will be a somewhat more complicated structure, but it still can be implemented efficiently. This will give you a more uniform distribution of samples across the periods of components that have frequencies related to your sample rate. Alternatively, use a resampling scheme that allows you to select an arbitrary resampling ratio and then resample by an (approximately) irrational number, like $\pi$. This can be done efficiently using a Farrow interpolator, which uses polynomial interpolation.

Answer 6

You need to smooth the histogram (this will yield similar results as using a kernel method). Exactly how the smoothing has to be performed needs experimentation. Maybe it could also be done by interpolation. In addition to the smoothing I believe you will also get improved results if you upsample your waveform in such a way that the sampling frequency is 'significantly higher' than the highest frequency in your input. This should help in the 'tricky' case where a sine wave is related to the sampling frequency in such a way that only a few bins in the histogram gets populated. If taken to the extreme a sufficiently high sample rate should give you nice plots without smoothing. So upsampling combined with some kind of smoothing should yield better plots.

You give an example of a 1kHz tone, where the plot is not as you expect. Here is my proposal (Matlab/Octave code)

pixels_vertical = 100;
% This needs to be tuned to your configuration and acceptance
upsampling_factor = 16*(pixels_vertical/100); 
fs_original = 48000;
fsine = 1000; % in Hz
fs_up = upsampling_factor*fs_original;
duration = 1; % in seconds
x = sin(2*pi*fsine*[0:duration*fs_up]/fs_up);
period_in_samples = fs_up/fsine;
hist_points = linspace(-1,1,pixels_vertical);
istart = 1;
iend   = period_in_samples;
pixel_values = hist(x(istart:iend), hist_points);
% smooth pixel values
[b,a] = butter(2,0.2);
pixel_values_smooth = filtfilt(b,a,pixel_values);
figure;hold on;
plot(hist_points, pixel_values);
plot(hist_points, pixel_values_smooth,'r');

For your 1000Hz tone you get this

What you need to do is to tune the upsampling_factor expression to your preference.

Still not 100% sure exactly what your requirements are. But using the above principle of upsampling and smoothing you will get this for the 1kHz tone (made with Matlab). Note that in the raw histogram there are many bins with zero hits.