双精度 C# 的丢失答案

【问题标题】：Loss of double precision C#双精度 C# 的丢失
【发布时间】：2012-11-06 21:10:48
【问题描述】：

我在使用双精度时遇到精度损失，我似乎无法找到精度丢失的位置。我正在编写一个软件合成器，由于某种原因，振荡器（声波发生器）的输入频率在某处出现严重混叠。振荡器应该产生三角波形，并且正确的频率被传递给方法，但是当结果在扬声器中播放时，我可以清楚地听到声音捕捉到不同的频率。

我正在使用 NAudio (http://naudio.codeplex.com/)，这个方法对我想要生成的每个样本运行一次。

代码如下：

double counter = 0;
int samplecounter = 0;

public double GetNextTriangle(double frequency)
{
    double samplesperwave = (double)Parent.WaveFormat.SampleRate / frequency;
    double length = (double)samplecounter / samplesperwave;
    if (length.CompareTo(0.25) == -1)
    {
        counter = length * 4.0;
    }
    else if (length.CompareTo(0.75) == -1)
    {
        counter = 1.0 - (length - 0.25) * 4.0;
    }
    else
    {
        counter = -1.0 + (length - 0.75) * 4.0;
    }
    samplecounter++;
    samplecounter = samplecounter > samplesperwave ? 0 : samplecounter;
    return counter;
}

提前致谢！ //垫子

【问题讨论】：

请使用< 和> 运算符而不是CompareTo！！
@Dai 这样代码可以很容易被人类理解。
@DOOMDUDEMX 无论如何，如果你必须使用 CompareTo，你应该使用 if (a.CompareTo(b) < 0) 而不是 if (a.CompareTo(b) == -1)，因为如果接收者小于参数，该方法的协定允许它返回任何小于零的值.
@DOOMDUDEMX 如果您无法通过阅读代码来调试它，那么您只需绘制输出并将其与您的预期进行比较。问题不在于精度。当然。问题是您没有正确编码方程式。
除了下面已经提供的几个优秀答案之外，值得考虑将 WaveTable 用于软件合成器。使用这种方法，您可以预先计算波形，然后对其进行索引（如果您需要支持更改频率或滑音，则可以在表中的条目之间进行插值）

标签： c# casting floating-point double naudio

【解决方案1】：

这里的浮点运算条件良好。但如果 1/frequency 不是 1/sampleRate 的倍数，您可能会看到与采样相关的别名问题。

如果可以，试试这个 matlab 代码

sampleRate=5000;
frequency=700;
sampleperwave=sampleRate/frequency;
samplecounter=0:floor(sampleperwave);
samplecounter=repmat(samplecounter,1,5);
length=samplecounter/sampleperwave;
wave=-1+4*(length-0.75);
wave(length<0.75)=1-4*(length(length<0.75)-0.25);
wave(length<0.25)=4*length(length<0.25);
figure; stem(wave); hold on; plot(wave,'r')

您可以尝试将 samplecounter 声明为 double 并使用模数递增

samplecounter++;
samplecounter = samplecounter % samplesperwave;

或者返回matlab

samplecounter=0:length(samplecounter)-1;
samplecounter=rem(samplecounter,sampleperwave);
len=samplecounter/sampleperwave;
wave=-1+4*(len-0.75);
wave(len<0.75)=1-4*(len(len<0.75)-0.25);
wave(len<0.25)=4*len(len<0.25);
figure; stem(wave); hold on; plot(wave,'r')

【讨论】：

将 SampleCounter 更改为 double 效果很好！我确实使用 TrackBar-Control 的位置值实时绘制了这些值，但是使用我的方法很难发现这种轻微的不一致。这解释得很好！感谢您的回答！

【解决方案2】：

您的问题不是精度问题。问题是您的函数没有定义三角波。每次将samplecounter重置为0，函数的下一个返回值为0。下一次round length设置为0，执行第一个if分支，计数器设置为0。

当遇到这样的问题时，您应该绘制输出。如果你这样做了，你会立即看到你的函数不会产生三角波形。下面是一些代码：

static double GetNextTriangle(double frequency)
{
    double samplesperwave = Parent.WaveFormat.SampleRate / frequency;
    double t = samplecounter / samplesperwave;
    counter = 1.0 - 4.0 * Math.Abs(Math.Round(t - 0.25) - (t - 0.25));
    samplecounter++;
    return counter;
}

再次强调，您必须绘制和可视化代码的输出以深入了解其行为。

【讨论】：

该函数正在生成样本，这些样本共同构成声波。当波浪完成一个循环时，我将计数器设置回 0，以便该过程可以重新开始。该方法应该返回一个介于 -1 和 1 之间的值。如果我不将计数器设置为 0，则结果值将继续增长超出范围。
这不会改变任何事情。不连续性是您声音混乱的原因。
+1，当波周期（1/频率）不是采样周期的倍数时，您将遇到潜在混叠的采样问题。例如，频率=700Hz，采样率=5kHz。通过将 samplecounter 重置为零，您可以将混叠问题转换为每个周期的不连续问题（通过采样，它将在达到最终零之前停止）。将 samplecounter 表示为 double 并使用模(samplecounter+1,sampleperwave)。
@aka 这个评论是针对我的吗？我的答案中的代码没有别名。
@DavidHeffernan 很抱歉造成混淆，我的意思是 DOOMDUDEMX 存在采样问题，感谢您揭露与此别名问题相关的不连续性 +1

【解决方案3】：

为了更好地了解问题所在，让我们将数据可视化。这是每秒 41 个样本的 2 Hz 波的两个周期：

样本 20 周围有一个光点。看起来样本没有正确拟合三角波。

让我们回顾一下三角波的数学。您正在做的是在离散点对连续函数进行采样。首先，定义一个频率为f（或1/T）的连续三角波，单位为赫兹：

tri(t) =   {  4*f*t,            0     <= t < T/4   }
       =   { -4*f*(t - T/2),    T/4   <= t < 3*T/4 }
       =   {  4*f*(t - T),      3*T/4 <= t < T     }
       =   tri(t - T)  [it's periodic!]

您现在想要对这个连续函数进行采样。所以你定义了一个采样率，s（或1/U），以每秒采样数为单位。现在，nth 样本将只是 tri(n*U)：

tri[n] = tri(n*U)
       =   {  4*f*n*U,            0     <= n*U < T/4   }
       =   { -4*f*(n*U - T/2),    T/4   <= n*U < 3*T/4 }
       =   {  4*f*(n*U - T),      3*T/4 <= n*U < T     }

让我们通过定义归一化周期 P = T/U 和归一化频率 F = f/s = U/T 来稍微清理一下：

tri[n] =   {  4*F*n,            0     <= n < P/4   }
       =   { -4*F*(n - P/2),    P/4   <= n < 3*P/4 }
       =   {  4*F*(n - P),      3*P/4 <= n < P     }

现在我们得到这个：

不必担心提示处明显的“光点”；它们是意料之中的，你不应该试图避免它们。

代码如下：

public static double GetNextTriangle(int sample, double frequency, double sampleRate)
{
    double T = 1d / frequency;
    double U = 1d / sampleRate;

    double P = T / U;
    double F = U / T;

    double n = (double)sample;
    n %= P; // restrict n to the domain [0, P)

    if ((n >= 0) && (n < (P / 4d)))
    {
        return 4d * F * n;
    }
    else if ((n >= (P / 4d)) && (n < (3d * P / 4d)))
    {
        return -4d * F * (n - (P / 2d));
    }
    else // if ((n >= (3d * P / 4d)) && (n < P))
    {
        return 4d * F * (n - P);
    }
}

【讨论】：

【解决方案4】：

整个方法似乎是错误的。我会选择这样的：

public class RectWave
{ // triangular wave in range [0,1] with given frequency
    public double GetNextSample()
    {
        if (_ptr == _wave.Length)
            _ptr = 0;
        return _wave[_ptr++];
    }

    public RectWave(int sampleRate = 441000, int period = 20)
    { 
        _sampleRate = sampleRate;
        _period = period;
        _wave = new double[(int)(_sampleRate * _period / 1000.0)];  // one cycle
        _ptr = 0;
        BuildWave();
    }

    private void BuildWave()
    {
        double dt = 1000.0 / _sampleRate;   // in msec
        double slope = 1.0 / (_period / 2.0);   // 1.0 => peek 
        for (int i=0; i <= _wave.Length/2; ++i)
        {
            _wave[i] = _wave[_wave.Length - i - 1] = i * dt * slope; 
        }
    }

    private int _sampleRate;    // in Hz
    private int _period;        // in msec
    private double[] _wave;     // the data
    private int _ptr;           // current sample
}

【讨论】：