有限浮点精度和无限谐波信号生成问题
Limited float precision and infinitely harmonic signal generation problem
假设我们需要生成一个非常长的谐波信号,理想情况下是无限长。乍一看,解决方案似乎微不足道:
样本 1:
float t = 0;
while (runned)
{
float v = sinf(w * t);
t += dt;
}
不幸的是,这是一个行不通的解决方案。对于 t >> dt
,由于浮点精度有限,将获得不正确的值。幸运的是我们可以记住 sin(2*PI* n + x) = sin(x)
其中 n - 任意整数值,因此修改示例不难得到一个 "infinite" 类比
样本 2:
float t = 0;
float tau = 2 * M_PI / w;
while (runned)
{
float v = sinf(w * t);
t += dt;
if (t > tau) t -= tau;
}
对于一次物理模拟,我需要得到一个无限大的信号,即谐波信号之和,像这样:
样本 3:
float getSignal(float x)
{
float ret = 0;
for (int i = 0; i < modNum; i++)
ret += sin(w[i] * x);
return ret;
}
float t = 0;
while (runned)
{
float v = getSignal(t);
t += dt;
}
在这种形式下,代码无法正确处理大型 t
,原因与 Sample1 类似。问题是 - 如何获得 Sample3 算法的 "infinite" 实现?我假设解决方案应该看起来像 Sample2。一个非常重要的注意事项——一般来说,w[i] 是任意的而不是谐波,也就是说,所有频率都不是某个基频的倍数,所以我找不到常见的 tau
。不允许使用精度更高的类型(double、long double)。
感谢您的建议!
您可以选择一个任意的 tau
并在从 t
中减去它时存储每个 mod 的阶段提醒(正如@Damien 在评论中建议的那样)。
此外,将时间表示为 t = dt * it
其中 it
是一个整数可以提高数值稳定性(我认为)。
也许是这样的:
int ndt = 1000; // accumulate phase every 1000 steps for example
float tau = dt * ndt;
std::vector<float> phases(modNum, 0.0f);
int it = 0;
float t = 0.0f;
while (runned)
{
t = dt * it;
float v = 0.0f;
for (int i = 0; i < modNum; i++)
{
v += sinf(w[i] * t + phases[i]);
}
if (++it >= ndt)
{
it = 0;
for (int i = 0; i < modNum; ++i)
{
phases[i] = fmod(w[i] * tau + phases[i], 2 * M_PI);
}
}
}
假设我们需要生成一个非常长的谐波信号,理想情况下是无限长。乍一看,解决方案似乎微不足道:
样本 1:
float t = 0;
while (runned)
{
float v = sinf(w * t);
t += dt;
}
不幸的是,这是一个行不通的解决方案。对于 t >> dt
,由于浮点精度有限,将获得不正确的值。幸运的是我们可以记住 sin(2*PI* n + x) = sin(x)
其中 n - 任意整数值,因此修改示例不难得到一个 "infinite" 类比
样本 2:
float t = 0;
float tau = 2 * M_PI / w;
while (runned)
{
float v = sinf(w * t);
t += dt;
if (t > tau) t -= tau;
}
对于一次物理模拟,我需要得到一个无限大的信号,即谐波信号之和,像这样:
样本 3:
float getSignal(float x)
{
float ret = 0;
for (int i = 0; i < modNum; i++)
ret += sin(w[i] * x);
return ret;
}
float t = 0;
while (runned)
{
float v = getSignal(t);
t += dt;
}
在这种形式下,代码无法正确处理大型 t
,原因与 Sample1 类似。问题是 - 如何获得 Sample3 算法的 "infinite" 实现?我假设解决方案应该看起来像 Sample2。一个非常重要的注意事项——一般来说,w[i] 是任意的而不是谐波,也就是说,所有频率都不是某个基频的倍数,所以我找不到常见的 tau
。不允许使用精度更高的类型(double、long double)。
感谢您的建议!
您可以选择一个任意的 tau
并在从 t
中减去它时存储每个 mod 的阶段提醒(正如@Damien 在评论中建议的那样)。
此外,将时间表示为 t = dt * it
其中 it
是一个整数可以提高数值稳定性(我认为)。
也许是这样的:
int ndt = 1000; // accumulate phase every 1000 steps for example
float tau = dt * ndt;
std::vector<float> phases(modNum, 0.0f);
int it = 0;
float t = 0.0f;
while (runned)
{
t = dt * it;
float v = 0.0f;
for (int i = 0; i < modNum; i++)
{
v += sinf(w[i] * t + phases[i]);
}
if (++it >= ndt)
{
it = 0;
for (int i = 0; i < modNum; ++i)
{
phases[i] = fmod(w[i] * tau + phases[i], 2 * M_PI);
}
}
}