5

我有一个稍微不寻常的问题，但我试图避免重新编码 FFT。

一般来说，我想知道这一点：如果我有一个为 type 实现的算法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

然后，想法是构造一个LogFloats 列表，然后在 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，但我想我会征求意见。这是一个相当反复出现的问题：我有一个我想在日志空间中操作的股票算法（而且它只使用了一些操作，如“+”、“-”、“ ”、“/”等）。

这让我想起了用模板编写泛型函数，这样返回参数、参数等都是由同一类型构造的。例如，如果您可以对浮点数进行 FFT，那么您应该能够轻松地对复杂值执行一次（只需使用为复杂值提供必要操作的类）。

就目前而言，看起来所有 FFT 实现都是为前沿速度而编写的，因此不会很普遍。所以到目前为止，看起来我必须为泛型类型重新实现 FFT ......

我这样做的原因是因为我想要非常高精度的卷积（并且N^2运行时间非常慢）。任何建议将不胜感激。

*注意，我可能需要为 实现三角函数LogFloat，这很好。

编辑： 这确实有效，因为LogFloat它是一个交换环（并且它不需要实现三角函数LogFloat）。最简单的方法是重新实现 FFT，但@JFSebastian 还指出了一种使用Python通用卷积的方法，它可以避免对 FFT 进行编码（同样，使用 DSP 教科书或Wikipedia 伪代码非常容易）。

4

2 回答 2

1

我承认我没有完全跟上你问题中的数学。但是，听起来您真正想知道的是如何处理极小和极大（绝对值）的数字而不会遇到 underflows 和 overflows。除非我误解了你，否则我认为这与我在处理货币单位时遇到的问题相似，不会因为四舍五入而在十亿美元的交易中损失几分钱。如果是这样，我的解决方案是 Python 的内置十进制数学模块。该文档适用于Python 2和Python 3. 简短的版本是十进制数学是一种任意精度的浮点和定点类型。Python 模块符合十进制数学的 IBM/IEEE 标准。在 Python 3.3（目前是 alpha 形式，但我一直在使用它没有任何问题）中，该模块已用 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 回答