我有一個稍微不尋常的問題，但我試圖避免重新編碼FFT。通用編程：日誌FFT或高精度卷積（python）
在一般情況下，我想知道這一點：如果我有一個用於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也很容易）。
我不知道你的問題是什麼。如果fft是用python編寫的，那麼上面應該（模仿雄心勃勃的瘋狂）工作。如果它調用了一個c實現，那麼它將無法工作，因爲c代碼就是c代碼，它不會執行python的操作。 那麼問題是什麼？ – 2012-03-10 01:25:15
有[通用卷積]（http://www.math.ucla.edu/~jimc/mathnet_d/sage/reference/sage/rings/polynomial/convolution.html）。 'LogFloat'可能不是處理下溢的最佳方法。 – jfs 2012-03-10 05:22:41
DSP教科書和教程網站中的大量FFT提供了非常簡單的FFT源代碼示例，通常大約有1或2頁代碼。 （在我的DSP網站上有幾個FFT只有30到40行BASIC） – hotpaw2 2012-03-10 05:50:38