У меня немного необычная проблема, но я пытаюсь избежать перекодирования 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 s, а затем использовать его в FFT:
x_log_float = [ LogFloat.from_float(x_val) for x_val in x ]
fft_x_log_float = fft(x_log_float)
Это определенно можно сделать, если я повторно внедряю БПФ (и просто использую LogFloat везде, где бы я использовал float, но я подумал, что буду просить совета. Это довольно повторяющаяся проблема: У меня есть стандартный алгоритм, который я хочу использовать в пространстве журнала (и он использует только несколько операций, таких как '+', '-', '', '/' и т. Д. ).
Это напоминает мне о написании универсальных функций с шаблонами, так что возвращаемые аргументы, параметры и т. Д. Создаются из одного и того же типа. Например, если вы можете сделать БПФ с плавающей точкой, вы сможете легкосделайте одно для комплексных значений (просто используя класс, который обеспечивает необходимые операции для комплексных значений).
В его нынешнем виде, похоже, что все реализации FFT написаны для передовой скорости, и поэтому выиграл 'быть очень общимEral.Итак, на данный момент, похоже, мне пришлось бы переопределить FFT для универсальных типов ...
Причина, по которой я это делаю, заключается в том, что мне нужны очень высокоточные свертки (и N^ 2 время выполнения крайне медленное).Любой совет будет принята с благодарностью.
* Обратите внимание, мне может понадобиться реализовать тригонометрические функции для LogFloat, и это было бы хорошо.
РЕДАКТИРОВАТЬ: Это работает, потому что LogFloat является коммутативным кольцом (и не требует реализации тригонометрических функций для LogFloat).Простейшим способом сделать это было переопределение FFT, но @JFSebastian также указал способ использования универсальной свертки Python, который позволяет избежать кодирования FFT (что, опять же, было довольно легко, используя либо учебник по DSP, либо псевдокод Википедии ).