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

質問

特定の（単純な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で計算し直すとか、適当に調整すれば良さ気ですかね？

回答No.3

maru_yoshi_ · ベストアンサー

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

さて、本題。
> この場合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)
{
T d = 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

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

お礼

どうもどうも、ご丁寧にありがとうございます。

＞誤差の比率は「相対誤差」という言い方が誤解がなくで済みそうです。

ローカルな用語については、そもそも込み入ったところの話なので、どう定義しても通じれば問題はないと思います。

私もその後funcの中身を書き換えて0.99fの500乗くらいまでは試していましたが

powと乗算の繰り返しの（絶対値の）誤差チェックにシンプルに特化したmaru_yoshi_さん作成のfuncの誤差チェック法でも

func(0.999f);
func(0.99);
func(0.999);

などの2048乗や16384乗とかでも大差はないようで、安心しました。

よーく考えてみましたが
「おそらく～～～とすれば、2nx以内には収まるんではないか」より正確だと考えられるのは

浮動小数の、仮数部のみに注目した場合、仮数部のみで表現できる最小精度（といいますか単位と言いますか）を

yとしたとき（yは符号部も指数部も無視しているので、一意の実数や絶対値ではない）

最終的な「仮数部の誤差」が2ny以内（私は詳しいことは知らないのですが、丸めの方法が四捨五入的なことやってくれてればny以内、かな）に収まるんではないか、という感じですかね。

この論法なら、相当0に近づいたとき

従って、指数部と仮数部がともに限界すれすれまで行ったとき「2nx」は遥かに超えて比率が跳ね上がる可能性が高い、ということも説明できるのではないかと思います。
（ただし0に近づいたときの相対誤差が跳ね上がったとしても、絶対値としては微小な差にすぎないので、結局は大した問題ではないでしょうが）

もしかしたら何かまだ見落としがあるかもしれないけど、概ね2ny以内には収まるということであれば
使用する浮動小数の仮数部に使用されるビット数と乗算回数から

おそらく、とり得る最大相対誤差は「2の-b乗」付近として十分予想可能という事になるはずです。

まぁ、まだ見落としがあったとしても実際には2nyよりだいぶ小さくなることも多々あるでしょうし、今回の状況では実用上の問題はやはりまずないだろう、ということで
解決とさせていただきます。

重ねて、どうもありがとうございました♪

回答No.2

maru_yoshi_

>＞最初の数値には丸め誤差が発生しており、その数値を使って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
ということです。

お礼

なるほど、そう言う事でしたか。
あとちょこっとだけ確認させてください
（最も影響が大きいところを見れば、それで近似できるであろう、という見解ではないかと思うのですが）

質問で減衰率をδとしているので
紛らわしくないように

Δ = a * x

と、xで表記する場合

3乗では
A^3=a^3( 1 + 3x + 3x*x + x*x*x )

となり

以下係数は

パスカルの三角形

であらわすことが出来る

このとき、n回の乗算を行うなら

a^n の内側について

単なる1のところが基と同じで
次にx<1から
nx が最も誤差の数値に対して大きく影響する（かもしれない）

ただし

乗算回数がかなり増えた場合

パスカルの三角形の関係上

たった10回で

1,10,45,120,210,252,210…

までになってしまうので
今度は係数が余りに大きくなると
その後の項の影響が大きくなってしまうかもしれない。

こういったことを踏まえて下記結果になった。

ということで、やはり結果的には、その時の様子見で決めるべき
ということですかね。

なお、先ほどのfunc関数の最後に

sprintf_s( c, 64, "%.40f\r\n", powf( f, (float)n ) );
OutputDebugStringA( c+left+n-1 );

を追加して再度実行してみると、このようになりました。

10000000149011612
10000000474974513
10000000474974513

10000000149011612
10000001715582130
10000000827403710

10000000149011612
10000003477346285
10000002153857305

10000000149011612
100000083
100000046

同じ場合もあるけど、概ねpowfで再計算する方が誤差が少なくなる
と考えてよさそうですね。

補足

いやまてよ

この場合xは相当小さいから
そんな程度の回数でこの結果ではないですよね。

と思ったら

丸め誤差が発生する（可能性がある）のもn回でしたねw

しかもそれは絶対値は一律ではないし、+方向か-方向かも一律では…ない（…のかな？）
ただし浮動小数の表現できる精度はfloatもdoubleも決まっていることから

仮に一定で正または負に一律に伸びるとして計算しても
最終的な誤差(の比率)は
2nx以内には収まるんではないかと考えて良さ気、でしょうかね。

回答No.1

maru_yoshi_

> で、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関数の差であり、真の結果との誤差ではないことに注意！

お礼

まず　pow( e, δ); →　pow( e, -δ);　か　1.0 / pow( e, δ);　とかでしたねｗ　
すみません。

ご回答ありがとうございます♪

大よその考え方としてはあってる、っぽいということでしょうかね。

＞一周期先の振幅はいいんですが、周期の途中はどうするんですか？
＞周期の途中も計算しないと、波形は連続にならないんですけど。

質問後にイメージしてたら、ちょうどそう思ってたとこでした

解決法としては、「既知のパラメータ」のなかに、振幅変化に対する分解能を示せるパラメータを含めばいいのではないかと、暫定しました。

今回の要件はWAVE（音声）ファイルに対してなのですが

例えば
サンプリング周波数48000Hz
でsin波、400Hzの波形を扱う場合
48000/400 = 240(サンプル)
が１周期に相当しますので
const double d = pow( e, -δ/240.0 );
としたのを、1サンプルごとにかけていけば、最大限の分解能が得られます。

ただ、そこまでの連続性を求めなくても
sin波、cos波、矩形波や三角波　…などなど
周期的に0(付近)になる箇所が分かってるものであれば

今までの経験上
そこを狙って、サンプル数分の計算結果を一括で乗算して振幅変化をさせる程度の分解能でも、十分それっぽく聞こえるはずです。

なお、この「連続性」
というのは量子化ビット数にも影響を受けます

そもそも量子化した時点で「完全な連続」でないのは確かなのですが、今回、整理してみると

このフィルタを通す段階では「ノイズが混じったように聞こえなければ問題ない」という状況でした。
つまり、よほどの急激な変化とならなければ大丈夫です。

また
仮に16bit符号付き整数でのデータであれば、電圧は
65536段階での調整となるので、それより細かい差になってくると、元々入力側が表現できてないところからフィルタをかけることになる、ので問題なし、となります。

（実際には、キャッシュミスと型変換のコストのバランスを考慮して、波形もパラメータも多大な計算時を行う時はfloatで扱う可能性が高いのですが（なのでpowf））

＞乗除算では丸め誤差がほとんどなのでそれほど問題は発生しないはずです。

ややっ？
そうなんですか？

浮動小数は内部的に2進でやってる、という事ぐらいしか知らなかったので
丸めとかの正確な内部処理をどうやってるかは知らなかったのですが

＞浮動小数点で掛け算を行った時に累積する誤差は乗算回数に比例します。
＞最初の数値には丸め誤差が発生しており、その数値を使って1000回乗算を行うと誤差は1000倍に拡大します。

すくなくとも、浮動小数の丸め誤差、ではない、現実の実数の場合は

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

となり、一致しないと思うのですが

浮動小数の丸め誤差
では、これが一致するような計算になってくる
ということでしょうか？

＞事前に d をN回かけた結果と～

全くその通りですね。

これを使いたい理由は、説明には非常に長い文字数が必要になってしまいますが
冷静に分析してみたら、ある程度ラフでもそこそこ大丈夫な理由がある
と判明しました。

ただ、上記疑問

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

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

補足

試してみました。（ほとんどCのC++です）
手抜き関数ですみません

void func( const int n, const float f, const int left ){

char c[ 64 ];
sprintf_s( c, 64, "%.40f\r\n", f );
OutputDebugStringA( c + left );

float d = f;
for ( int i=n-1;i--; ) d *= f;

sprintf_s( c, 64, "%.40f\r\n", d );
OutputDebugStringA( c + left + n - 1 );

}

後半の0は省いて書くと、この結果は以下のようになりました

func( 4, 0.1f, 2 );
の場合

10000000149011612
10000000474974513

func( 8, 0.1f, 2 );

10000000149011612
10000001715582130

func( 16, 0.1f, 2 );

10000000149011612
10000003477346285

func( 32, 0.1f, 2 );

10000000149011612
100000083

という事で、これはあくまで
10進表記で表示させているのであれですが

とはいえ、最初の誤差に対する単純な乗算回数倍よりかはだんだん差は開いていく（？）ものの

乗算回数倍というのは一つの目安にはなりそうですね。

って、こういうことであっているのでしょうか？