从样本中计算波形的 PDF

信息处理 算法 插值 声音的
2021-12-29 23:00:55

前段时间我在尝试不同的方法来绘制数字波形,我尝试的其中一件事是,而不是幅度包络的标准轮廓,将其显示得更像示波器。这是示波器上正弦波和方波的样子:

在此处输入图像描述

天真的方法是:

  1. 将音频文件划分为输出图像中每个水平像素的一个块
  2. 计算每个块的样本幅度直方图
  3. 按亮度将直方图绘制为一列像素

它产生这样的东西: 在此处输入图像描述

如果每个块有很多样本并且信号的频率与采样频率无关,则此方法可以正常工作,但其他情况则不然。例如,如果信号频率是采样频率的精确约数,则样本在每个周期中总是以完全相同的幅度出现,并且直方图将只是几个点,即使这些点之间存在实际的重构信号。这个正弦脉冲应该和左上图一样平滑,但这不是因为它正好是 1 kHz,并且样本总是出现在相同的点附近:

在此处输入图像描述

我尝试上采样以增加点数,但这并不能解决问题,只是在某些情况下有助于解决问题。

所以我真正想要的是一种从其数字样本(幅度与时间)中计算连续重建信号的真实PDF (概率与幅度)的方法。我不知道为此使用什么算法。一般来说,函数的 PDF 是其反函数的导数。

sin(x) 的 PDF:ddxarcsinx=11x2

但我不知道如何计算反函数为多值函数的波,或者如何快速计算。把它分解成分支并计算每个分支的倒数,取导数,然后把它们加在一起?但这很复杂,可能有更简单的方法。

这个“插值数据的 PDF”也适用于我对 GPS 轨迹进行核密度估计的尝试。它应该是环形的,但是因为它只查看样本而不考虑样本之间的插值点,所以 KDE 看起来更像是一个驼峰而不是一个环。如果样本是我们所知道的,那么这是我们能做的最好的。但样本并不是我们所知道的全部。我们也知道样本之间有一条路径。对于 GPS,没有像带限音频那样完美的奈奎斯特重建,但基本思想仍然适用,插值函数有一些猜测。

4个回答

插值几倍于原始速率(例如 8x 过采样)。这允许您假设分段线性信号。与波形的无限分辨率、连续 sin(x)/x 插值相比,该信号的误差非常小。

假设每一过采样值都有一条从一个值到下一个值的连续线。使用之间的所有值。这为您提供了一个从 y1 到 y2 的薄水平切片,可以累积成任意分辨率的 PDF。每个矩形概率切片必须缩放到 1/nsamples 区域。

使用样本之间的线而不是样本本身可以防止“尖峰”PDF,即使在采样周期和波形之间存在基本关系的情况下也是如此。

我要使用的基本上是 Jason R 的“随机重采样器”,它又是 yoda 随机采样的基于预采样信号的实现。

我对每两个样本之间的一个随机点使用了简单的三次插值。对于原始合成声音(从饱和的非带限方形信号 + 偶次谐波衰减到正弦),它看起来像这样:

随机重采样合成器 PDF

让我们将其与更高采样的版本进行比较,

在此处输入图像描述

以及具有相同采样率但没有插值的怪异。

在此处输入图像描述

这种方法的显着伪影是平方域中的过冲,但这实际上是 sinc 滤波信号的 PDF (正如我所说,我的信号没有带限)看起来也更好地代表了感知响度比峰值,如果这是一个音频信号。

代码(哈斯克尔):

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是 [0,1] 范围内的随机变量列表。

虽然您的方法在理论上是正确的(并且需要对非单调函数稍作修改),但计算泛型函数的逆非常困难。正如您所说,您将不得不处理分支点和分支切割,这是可行的,但您真的不想这样做。

正如您已经提到的,常规采样对同一组点进行采样,因此在不采样的区域中极易受到不良估计的影响(即使满足奈奎斯特标准)。在这种情况下,较长时间的采样也无济于事。

一般来说,在处理概率密度函数和直方图时,从随机抽样的角度考虑比常规抽样要好得多(请参阅链接的答案以获取介绍)。通过随机抽样,您可以确保每个点都具有相同的“命中”概率,并且是估计 pdf 的更好方法。

这是一个例子:考虑函数现在,如果我以采样频率 Hz(奈奎斯特频率, Hz)对其进行采样,则得到的概率密度是左侧的图(-2 和 2 之间的 401 个等间距箱)。如果我采样 10 秒或 100 秒都没关系。它仍然保持不变。另一方面,以每秒样本(均匀分布)的速率随机采样 30 秒(我在这里没有使用 Hz,因为这意味着不同的含义)给出了右侧的图(相同的分箱)。f(x)=sin(20πx)+sin(100πx)fs=1000fN=1001000

您可以很容易地看到,虽然它很嘈杂,但它比右边的 PDF 更接近实际 PDF,右边的 PDF 在几个间隔中显示零,在其他几个间隔中显示大错误。通过更长的观察时间,您可以降低右侧的方差,最终在大型观察的范围内收敛到精确的 PDF(黑色虚线)。

在此处输入图像描述

核密度估计

估计波形 PDF 的一种方法是使用核密度估计器

这需要你所有的样本和一个核函数(例如高斯),并将它们与卷积(狄拉克增量)给出 PDF 的估计值,x(n)K(x)δ(xx(n))P^

P^(x)=n=0NK(xx(n))


更新:有趣的附加信息。

假设您有那么 --- 正如您所说 --- 您还可以知道它的 DFTx(n)n=0,1,...,N1X(k)

X(k)=n=0N1x(n)eȷ2πnk/N

所以的系数:X(k)eȷ2πnk/N

x(n)=1Nk=0N1X(k)eȷ2πnk/N

因此,猜测您可能会将每个傅立叶分量的所有 PDF 卷积在一起:

|X(k)|11x2

然而,这并没有考虑到中的加法有贡献(或没有)的方式X(k)x(n)

不过,需要更多的思考!