通用编程：日志FFT或高精度卷积（python）

Question

我有一个稍微不寻常的问题，但我试图避免重新编码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过，但我想我会请教这是一个相当反复出现的问题：我有一只股票的算法，我想在日志空间操作（和它仅使用像“操作了一把+ '，' - '，''，'/'等）。

这让我想起使用模板编写泛型函数，以便返回参数，参数等由相同类型构造而成。例如，如果你可以做一个浮点数的FFT，你应该能够很容易地对复数值做一个（通过简单地使用一个为复数值提供必要操作的类）。

当前标准的，它看起来像所有的FFT实现是带血的速度写的，所以不会很一般。因此，截至目前，它看起来像我不得不重新实现FFT泛型类型...

的原因，我这样做是因为我要非常高精度的卷积（和N^2运行非常缓慢）。 任何意见将不胜感激。

*注意，我可能需要实现三角函数的LogFloat，这将是罚款。

编辑： 这并不工作，因为LogFloat是一个交换环（和它不要求执行的三角函数的LogFloat）。最简单的方法是重新实现FFT，但@ J.F.Sebastian还指出了一种使用通用卷积的方法，该方法避免了对FFT进行编码（使用DSP教科书或使用DSP教科书或the Wikipedia pseudocode也很容易）。

Answer 1

我承认我并没有完全跟上你的问题的数学。然而，这听起来像你真正想知道的是如何处理极小和大（绝对值）的数字，而不会造成下溢和溢出。除非我误解了你，否则我认为这与我使用货币单位的问题类似，而不是因四舍五入而损失数十亿美元交易的便士。如果是这样，我的解决方案是Python的内置十进制数学模块。该文档适用于Python 2和Python 3。简短版本是十进制数学是一种任意精度的浮点和定点类型。 Python模块符合IBM/IEEE的十进制数学标准。在Python 3.3中（目前使用的是alpha版本，但我一直在使用它，完全没有问题），该模块已被C语言重写，速度提高了100倍（在我的快速测试中）。

Answer 2

您可以通过大量的规模的时域样本，以避免下溢，然后，如果

˚F（F（T））= X（j * W）

然后

˚F（SF（S * T））< - >X（W/S）

现在使用卷积定理可以解决如何扩展您的最终结果删除比例因子的影响。