少し変わった問題がありますが、FFTの再コーディングを避けようとしています。
一般的に、私はこれを知りたいです:タイプに対して実装されたアルゴリズムがあるが、特定の操作のセットが定義されている場所(たとえば、、、、 ...float
も定義する複素数)で機能する+
場合、*
それらの操作をサポートする別のタイプでそのアルゴリズムを使用する最良の方法は?一般的に数値アルゴリズムは一般性ではなく速度のために書かれているため、実際にはこれは注意が必要です。
具体的
には、ダイナミックレンジが非常に高い値を処理しているため、ログスペースに保存したいと思います(主にアンダーフローを回避するため)。
私が欲しいのは、いくつかのシリーズのFFTのログです。
x = [1,2,3,4,5]
fft_x = [ log( x_val ) for x_val in fft(x) ]
これでも、かなりのアンダーフローが発生します。私が欲しいのは、ログ値を保存し、の代わりに、などの代わりに使用することです+
。*
logaddexp
+
これを行う方法についての私の考えは、これらのプリミティブ操作を定義する(ただし、ログスペースで動作する)単純なLogFloatクラスを実装することでした。次に、ログに記録された値を使用させることで、FFTコードを実行することができます。
class LogFloat:
def __init__(self, sign, log_val):
assert(float(sign) in (-1, 1))
self.sign = int(sign)
self.log_val = log_val
@staticmethod
def from_float(fval):
return LogFloat(sign(fval), log(abs(fval)))
def __imul__(self, lf):
self.sign *= lf.sign
self.log_val += lf.log_val
return self
def __idiv__(self, lf):
self.sign *= lf.sign
self.log_val -= lf.log_val
return self
def __iadd__(self, lf):
if self.sign == lf.sign:
self.log_val = logaddexp(self.log_val, lf.log_val)
else:
# subtract the smaller magnitude from the larger
if self.log_val > lf.log_val:
self.log_val = log_sub(self.log_val, lf.log_val)
else:
self.log_val = log_sub(lf.log_val, self.log_val)
self.sign *= -1
return self
def __isub__(self, lf):
self.__iadd__(LogFloat(-1 * lf.sign, lf.log_val))
return self
def __pow__(self, lf):
# note: there may be a way to do this without exponentiating
# if the exponent is 0, always return 1
# print self, '**', lf
if lf.log_val == -float('inf'):
return LogFloat.from_float(1.0)
lf_value = lf.sign * math.exp(lf.log_val)
if self.sign == -1:
# note: in this case, lf_value must be an integer
return LogFloat(self.sign**int(lf_value), self.log_val * lf_value)
return LogFloat(self.sign, self.log_val * lf_value)
def __mul__(self, lf):
temp = LogFloat(self.sign, self.log_val)
temp *= lf
return temp
def __div__(self, lf):
temp = LogFloat(self.sign, self.log_val)
temp /= lf
return temp
def __add__(self, lf):
temp = LogFloat(self.sign, self.log_val)
temp += lf
return temp
def __sub__(self, lf):
temp = LogFloat(self.sign, self.log_val)
temp -= lf
return temp
def __str__(self):
result = str(self.sign * math.exp(self.log_val)) + '('
if self.sign == -1:
result += '-'
result += 'e^' + str(self.log_val) + ')'
return result
def __neg__(self):
return LogFloat(-self.sign, self.log_val)
def __radd__(self, val):
# for sum
if val == 0:
return self
return self + val
次に、のリストを作成し、LogFloat
それをFFTで使用するというアイデアがあります。
x_log_float = [ LogFloat.from_float(x_val) for x_val in x ]
fft_x_log_float = fft(x_log_float)
これは、FFTを再実装すれば間違いなく実行できます(以前LogFloat
はどこでも使用できfloat
ますが、アドバイスを求めたいと思いました。これはかなり繰り返し発生する問題です。ログスペースで操作したいストックアルゴリズムがあります。 (そして、'+'、'-'、' '、'/'などの少数の操作のみを使用します)。
これは、テンプレートを使用してジェネリック関数を作成することを思い出させます。そのため、戻り引数、パラメーターなどは同じ型から構築されます。たとえば、floatのFFTを実行できる場合は、複素数値に対して簡単にFFTを実行できるはずです(複素数値に必要な操作を提供するクラスを使用するだけです)。
現在のところ、すべてのFFT実装は最先端の速度で作成されているように見えるため、あまり一般的ではありません。したがって、現時点では、ジェネリック型のFFTを再実装する必要があるようです...
私がこれを行っている理由は、非常に高精度の畳み込みが必要なためです(そして、N ^ 2ランタイムは非常に遅いです)。アドバイスをいただければ幸いです。
*注:の三角関数を実装する必要があるかもしれませんが、それでLogFloat
問題ありません。
編集:
これはLogFloat
可換環であるため機能します(そして、の三角関数の実装を必要としませんLogFloat
)。これを行う最も簡単な方法はFFTを再実装することでしたが、@ JFSebastianは、FFTのコーディングを回避する一般的な畳み込みの使用方法も指摘しましたPython
(これもDSP教科書またはウィキペディア擬似コードを使用すると非常に簡単でした)。