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:

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.

• Do you have an example of a multivalued function that you're interested in? You'll probably have to evaluate it along a branch cut that makes the most sense for your physical data. – Lorem Ipsum Sep 8 '11 at 21:10
• Are you more interested in ways to draw that kind of plot, or is the plot just motivation for the question about calculating the PDF? – datageist Sep 8 '11 at 21:17
• @yoda: Well, the function above for sine wave is found by taking only half a cycle, inverting and taking the derivative, because each half-cycle has the same PDF as the next. But to get the value for an entire arbitrary audio signal, you can't make that assumption. I think you'd need to divide it into "branch cuts", take the PDF of each in turn, and sum them all together? – endolith Sep 8 '11 at 21:28
• @datageist: Hmm. I'm interested in ways to draw that kind of plot, but that kind of plot is the PDF. A shortcut that produces the same or very similar result is ok. – endolith Sep 8 '11 at 21:32
• @endolith, Oh yeah, I understand. Just a question about emphasis really (i.e. which kinds of shortcuts are reasonable). – datageist Sep 8 '11 at 21:49

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.

• I've written a function for the linearly-interpolated histogram, but it's dodgy. Do you know of existing code for this? – endolith Oct 7 '11 at 2:33
• The linear interpolation makes a huge difference for most waveforms, even without the oversampling. 1 kHz sine looks mostly like 997 Hz sine now. Instead of just horizontal lines at the sample values, it's now horizontal bands of color between them. With oversampling, the bands are smoothed out, too. With FFT resampling and some overlap with adjacent chunks I should be able to make it hit the real intersample peaks. I need to make my interpolated histogram code faster, though... – endolith Oct 10 '11 at 18:00
• I completely rewrote my script for this, and I think I got the histogram and antialiasing right this time: gist.github.com/endolith/652d3ba1a68b629ed328 – endolith Jul 16 '15 at 2:28
• The latest version is at github.com/endolith/scopeplot – endolith Jul 26 '19 at 13:59

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.

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

• Looks awesome. +1 for the Haskell code. – datageist Sep 12 '11 at 15:09
• Yes, it should overshoot the sample values. I actually planned to have a peak value for each pixel column, too, possibly drawn differently, based on max intersample peaks and not just on max samples. Waveforms like flic.kr/p/7QAScX show why this is necessary. – endolith Sep 12 '11 at 20:37
• By "higher-sampled version" do you mean that it's upsampled, but still uniformly-sampled? And that's the blue dots? – endolith Sep 16 '11 at 16:05
• @endolith It's simply the original waveform calculated in a higher sample rate in the first place. Essentially like the blue points represent a sound sampled at 192 kHz, and the downmost yellow ones represent a naïvely-done downsample to 24 kHz. The upper yellow points are stochasticAntiAlias of this. But the higher-sampled version is indeed uniform rate in both cases. – leftaroundabout Sep 16 '11 at 23:21

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.

• "it is extremely hard to calculate the inverse of a generic function" Well, this isn't a function so much as a series of samples, so finding the inverse is just swapping the x and y coordinates of the samples and then resampling to fit the new coordinate system. I can't change the sampling anyway. We're talking about pre-existing data created using uniform sampling. – endolith Sep 10 '11 at 3:39

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) )$$

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!

• I thought of that, but density estimation is used for estimating an unknown probability density function. Because of the Nyquist sampling theorem, the entire waveform is known, exactly, and the exact probability density function should be knowable, too. I'm fine with estimation if it's a speed vs accuracy trade-off, but there must be a way to get the actual PDF out of it. Like, a reconstructed waveform can be made by putting a sinc function at each sample and summing them together. Can the PDF be created by using the PDF of a sinc function as the kernel? I don't think it works like that. – endolith Sep 8 '11 at 20:26
• Like, I don't think this solves the problem where the signal samples are a submultiple of the sampling frequency. It doesn't take into account the reconstructed waveform between samples, does it? It just blurs each point in the PDF to try to fill in gaps. I had a similar issue with trying to do kernel density estimation of a GPS trace, because it doesn't take into account the values between samples. – endolith Sep 8 '11 at 20:34

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.

• But what if the waveform's at 44.1/π kHz? :) This is good advice, though. Is there such a thing as random resampling? Or really, I guess what would work perfectly would be to resample non-uniformly, such that the new samples fit perfectly in bins in y dimension, instead of being evenly-spaced in the x dimension. Not sure if there's a way to do that – endolith Sep 11 '11 at 14:22
• You could easily implement a "random" resampler using a Farrow structure. It's a scheme that allows for arbitrary fractional-sample delay by interpolating using polynomials (often cubic). You could maintain an inter-sample phase accumulator, similar to that used in an NCO, that is incremented by pseudorandom fractions of a sampling interval for each output (resampled) sample. The value of the accumulator is used as an input to the Farrow interpolator, defining the amount of fractional delay for each output. – Jason R Sep 12 '11 at 0:05
• Hmm, to clarify, Farrow is just a processor/memory optimized version of regular old polynomial interpolation? – endolith Sep 12 '11 at 20:57
• Yes. It's just an efficient structure for implementing polynomial-based arbitrary fractional delay. – Jason R Sep 13 '11 at 2:31
• Cubic interpolation is only an approximation, though. I want to know true intersample peaks, and it doesn't seem to work well on extreme peaks: stackoverflow.com/questions/1851384/… Actually, it seems that an infinite series with a discontinuity like [..., -1, 1, -1, 1, 1, -1, 1, -1, ...] will produce an infinite intersample peak, though, so I'm not sure how much this would matter in practice. – endolith Oct 3 '11 at 1:46

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.

• Yeah, it really needs some type of interpolation as part of the algorithm. Smoothing the histogram alone won't do it, because the histogram is of discrete points, not the reconstructed waveform. The only way upsampling would work is if I do it to the point where there are many more samples than there are vertical pixels, but that's a heavy brute force method that takes a long time. – endolith Sep 8 '11 at 20:37
• or calculating the effect of interpolating on the output without actually interpolating – endolith Sep 8 '11 at 20:43