対数減衰を実現するときの具体的手法について

Question

特定の（単純なsinやcos、矩形波なども含めて、とにかく基本周波数が既知な）波形を対数減衰させたい場合
対数減衰率などのパラメータは既知だとして

対数減衰の概念は対数減衰率をδ
ある地点での振幅をAn、その1周期先の振幅をAn+1（下付き添字）と表す場合

δ=ln(An/An+1)=ln(An+1/An+2)= … =ln(An+m-1/An+m)

なので

An/An+1=An+1/An+2= … =An+m-1/An+m

という解釈でいいのですよね？

ということは、δが既知で振幅を出すときは

An+1 = An/(e^δ)
An+2 = An+1/(e^δ)
・
・
・

で、e^δが一定
なので前の振幅にpow( e, δ);の戻り値をかければ次の振幅が出る

ただし、コストの高いpowを何度も呼び出すのは実用的ではないので

あらかじめ
double d = pow( e, δ);

とかやっといて、dを何度もかける

とかで、よほど厳密さを求めない限り（浮動小数を使う時点でそんなに厳密さを求めてもないような気もしますが）実用的には十分なのでしょうか？

ただ、何度も掛けると誤差が累積すると思うのですが
それが問題になるようなら
例えば
1000回に一回ぐらいpowで計算し直すとか、適当に調整すれば良さ気ですかね？

maru_yoshi_ · Accepted Answer

前回の私の以下の回答ですが、

> 誤差の比率（精度と言ってもよいかもしれない）が乗算回数に比例するというのが正しいでしょう。
これも用語の使い方がおかしいですね。
「精度が乗算回数に比例する」という言い方だと、計算を重ねると精度がよくなっていくように感じられてします。
誤差の比率は「相対誤差」という言い方が誤解がなくで済みそうです。

さて、本題。
> この場合xは相当小さいから
> そんな程度の回数でこの結果ではないですよね。
その通りです。
実際 0.1に対して0.1fは0.10000000149011612ですから相対誤差は1.49011612×10の-8乗。
Nがかなり大きくないとxの２乗以降の項は影響を与えません。

> 丸め誤差が発生する（可能性がある）のもn回でしたねw
これは忘れていました。
前回の計算結果は元々の数値が持っている誤差が計算によってどれだけ拡大するかの計算誤差だけを求めていました。
実際には計算結果も丸められるので、誤差がそれ以上に拡大するんですね。

> 2nx以内には収まるんではないかと考えて良さ気、でしょうかね。
たぶんそうなんでしょうね。
実際には手動で計算した値とプログラムで計算した値を比べて評価することが必要だと思いますよ。
特に今回の場合、N回計算後の値は計算機で簡単に求めることができますからね。

さて、元の問題である減衰率を計算する場合の誤差ですが、以下のプログラムでpowと掛け算の繰り返しの差を測定してみました。

#include <stdio.h>
#include <tchar.h>
#include <iostream>

// 乗算繰り返しによる累乗計算
template <typename T>
T my_power(T f, int n)
{
Td = 1.0;
for (int i = n; i != 0; --i)d *= f;
return d;
}

// 誤差測定
template <typename T>
void func(T delta)
{
cout << "type:"<< typeid(T).name() <<"\tdata=" << delta << endl;
cout << "count:\tpowered\t\tmy_powered\trelative_error" << endl;
for (int i = 4; i <= 65536; i <<= 1)
{
T powered= pow(delta, i);
T my_powered = my_power(delta, i);
T error = powered - my_powered;// 誤算
T relative_error = error / powered;// 相対誤差
cout << i << ":\t" << powered << "\t" << my_powered << "\t" << relative_error << endl;
if (powered < delta/32768.0)break;// 演算打ち切り
}
cout << endl;
}

int _tmain(int argc, _TCHAR* argv[])
{
// deltaが1に近い場合のみ検証する
func(0.99f);
func(0.999f);
func(0.99);
func(0.999);
return 0;
}

結果
type:float data=0.99
count: powered my_powered relative_error
4: 0.960596 0.960596 0
（途中省略）
256: 0.0763153 0.0763152 2.14784e-006
512: 0.00582403 0.00582401 3.35812e-006
1024: 3.39193e-005 3.39191e-005 8.15129e-006
2048: 1.15052e-009 1.1505e-009 1.77555e-005

type:float data=0.999
count: powered my_powered relative_error
4: 0.996006 0.996006 0
（途中省略）
4096: 0.0166065 0.0166059 3.24153e-005
8192: 0.000275775 0.000275757 6.58536e-005
16384: 7.60519e-008 7.60416e-008 0.000135565

type:double data=0.99
count: powered my_powered relative_error
4: 0.960596 0.960596 0
（途中省略）
256: 0.076315 0.076315 -2.72773e-015
512: 0.00582398 0.00582398 -4.46788e-015
1024: 3.39187e-005 3.39187e-005 -8.19096e-015
2048: 1.15048e-009 1.15048e-009 -1.67165e-014

type:double data=0.999
count: powered my_powered relative_error
4: 0.996006 0.996006 0
（途中省略）
4096: 0.016605 0.016605 3.55197e-014
8192: 0.000275727 0.000275727 7.21551e-014
16384: 7.60255e-008 7.60255e-008 1.4658e-013

この結果からすると、乗算の繰り返しで問題はなさそうですね。

maru_yoshi_ · Answer

>＞最初の数値には丸め誤差が発生しており、その数値を使って1000回乗算を行うと誤差は1000倍に拡大します。

>
> に対する部分に関しては解決しておきたいです。

私の言葉の使い方が間違ったかもしれません。
誤差の比率（精度と言ってもよいかもしれない）が乗算回数に比例するというのが正しいでしょう。

数値 A の真の値が a とした時、誤差Δ（デルタ）は
Δ = |a - A|（真の値との差の絶対値）
∴ A = a ± Δ（真の値は数値Aの誤差範囲内）
通常、誤差は正の値で取り扱うので
A = a + Δ
Δ = a * δと置くと
A = a + Δ = a * (1 + δ)
δは誤差が真の値に対する比率（精度）である。

Aを使って計算（乗算）を行うと
A*A = a*(1+δ)*a*(1+δ)
= a*a*(1+2*δ+δ*δ)
従ってA * A の計算結果はその真の値 a * a に対して a * a * (2 * δ + δ * δ)分の誤差を含む。
通常δ<<1と考えるので　2 * δ >> δ * δとなり δ * δの項は無視できる。
つまり、A*Aの誤差は a*a*2*δに近似でき、真の計算結果に対する誤差の比率は2*δになる。
つまり、誤差（の比率）が2倍になった。

同様にAのN乗は
(AのN乗)=(aのN乗)*((1+δ)のN乗)
=(aのN乗)*(1+N*δ+N*(N-1)*δ*δ/2 + ...)
≒(aのN乗)*(1+N*δ)
(AのN乗)-(aのN乗)≒(aのN乗)*N*δ
誤差（の比率）がN倍になった。

> すくなくとも、浮動小数の丸め誤差、ではない、現実の実数の場合は
>
> 10を8乗すると100000000
> 9を8乗すると43046721
>
> ( 100000000 - 43046721 ) > (10-9)*8
>
> 1を8乗すると1
> 0.9を8乗すると0.43046721
>
> ( 1 - 0.43046721 ) < (1-0.9)*8
>
> となり、一致しないと思うのですが
>
> 浮動小数の丸め誤差
> では、これが一致するような計算になってくる
> ということでしょうか？
10に対する9の差は1、精度は1/10で0.1、8乗数に対する精度は0.1*8＝0.8
10の8乗に対する9を8乗の差は(100000000 - 43046721)＝56953279、精度は0.56953279
0.56953279 < 0.8
が成り立つでしょ。つまり、
(100000000 - 43046721)/100000000 < ((10-9)/10)*8
ということです。

maru_yoshi_ · Answer

> で、e^δが一定
> なので前の振幅にpow( e, δ);の戻り値をかければ次の振幅が出る
一周期先の振幅はいいんですが、周期の途中はどうするんですか？
周期の途中も計算しないと、波形は連続にならないんですけど。

> よほど厳密さを求めない限り（浮動小数を使う時点でそんなに厳密さを求めてもないような気もしますが）実用的には十分なのでしょうか？
それはどの程度の精度を求めているのか不明なので、十分かどうかはわかりません。
δの具体的な値を使って、繰り返し回数、powの精度を検討した上でどの程度の誤差が発生するかを評価する必要があります。

> 何度も掛けると誤差が累積すると思うのですが
浮動小数点で掛け算を行った時に累積する誤差は乗算回数に比例します。
最初の数値には丸め誤差が発生しており、その数値を使って1000回乗算を行うと誤差は1000倍に拡大します。
倍精度の場合、丸め誤差は2の-52乗程度（10の-16乗程度）なので、1000回乗算を行った時の誤差は10の-13乗程度です。
もっとも今回は元の値はpow関数の結果なので既にそこにある程度の誤差があることを考慮に入れておく必要があります。
残念ながらpow関数の精度は私も知りません。

浮動小数点では誤差が累積するというのは誤差の性質を理解していないための誤解だと思います。
誤差が拡大するケースは桁落ちや情報落ちがほとんどであり、それらは加減算で発生します。
乗除算では丸め誤差がほとんどなのでそれほど問題は発生しないはずです。

> 1000回に一回ぐらいpowで計算し直すとか、適当に調整すれば良さ気ですかね？ 
事前に d をN回かけた結果と、pow(e, N * δ)の値を計算してみて誤差がどの程度になるか確かめてみたらいいんじゃないですかね。
Nは1024（2の10乗）程度から指数関数的に増やしていって、無視出来ない誤差になったら再計算というのはどうですか。
それでも、この方法でわかる差は乗算繰り返しとpow関数の差であり、真の結果との誤差ではないことに注意！

対数減衰を実現するときの具体的手法について

前回の私の以下の回答ですが、

>＞最初の数値には丸め誤差が発生しており、その数値を使って1000回乗算を行うと誤差は1000倍に拡大します。

この回答への補足

> で、e^δが一定

この回答への補足

似たような質問が見つかりました

関連するカテゴリからQ&Aを探す

デイリーランキングこのカテゴリの人気デイリーQ&Aランキング

マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング