5

少し変わった問題がありますが、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教科書またはウィキペディア擬似コードを使用すると非常に簡単でした)。

4

2 に答える 2

1

私はあなたの質問の数学に完全に追いついていないことを告白します。しかし、あなたが本当に疑問に思っているのは、アンダーフローやオーバーフローにぶつかることなく、非常に小さい数と大きい数(絶対値)をどのように処理するかということのように思えます。私があなたを誤解しない限り、これは私がお金の単位で作業している問題に似ていると思います。丸めによって数十億ドルの取引でペニーを失うことはありません。その場合、私の解決策はPythonの組み込みのdecimal-mathモジュールです。ドキュメントはPython2Python3の両方に適しています。短いバージョンでは、10進演算は任意精度の浮動小数点型と固定小数点型です。Pythonモジュールは、10進演算に関するIBM/IEEE標準に準拠しています。Python 3.3(現在はアルファ形式ですが、問題なく使用しています)では、モジュールがCで書き直され、最大100倍の速度が得られました(私のクイックテストで)。

于 2012-05-31T02:18:00.443 に答える
0

アンダーフローを回避するために、時間領域のサンプルを多数の s でスケーリングできます。

F ( f(t) ) = X (j*w)

それから

F (sf(s*t)) <-> X (w/s)

畳み込み定理を使用して、最終結果をスケーリングしてスケーリング係数の影響を取り除く方法を理解できます。

于 2012-06-12T05:15:32.843 に答える