Programación genérica: Log FFT O convolución de alta precisión (python)

Question

Tengo un problema levemente inusual, pero estoy tratando de evitar volver a codificar FFT.

En general, quiero saber esto: Si tengo un algoritmo que se implementa para el tipo float, pero funcionaría siempre que sea un cierto conjunto de operaciones se define (por ejemplo, números complejos, por lo que también definen +, *, ...), ¿cuál es la mejor manera de usar ese algoritmo en otro tipo que admita esas operaciones? En la práctica, esto es complicado porque, en general, los algoritmos numéricos están escritos para la velocidad, no para la generalidad.

Específicamente:
estoy trabajando con los valores con un rango dinámico muy alto, y por eso me gustaría almacenarlos en el espacio de registro (en su mayoría para evitar desbordamiento).

Lo que me gustaría es el logaritmo de la FFT de algunas series:

x = [1,2,3,4,5] 
fft_x = [ log(x_val) for x_val in fft(x) ] 

Incluso Esto dará lugar a desbordamiento significativo. Lo que me gustaría es almacenar los valores de registro y utilizar en su lugar de +*logaddexp y en lugar de +, etc.

Mi idea de cómo hacer esto fue implementar una clase LogFloat simple que define estas operaciones primitivas (pero opera en el espacio de registro). Entonces podría simplemente ejecutar el código FFT dejándolo usar mis valores registrados.

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 

Entonces, la idea sería construir una lista de LogFloat s y, a continuación, utilizarlo en la FFT:

x_log_float = [ LogFloat.from_float(x_val) for x_val in x ] 
fft_x_log_float = fft(x_log_float) 

Esto definitivamente se puede hacer si volver a implementar FFT (y simplemente use LogFloat donde quiera que use float, pero pensé que pediría consejo. Este es un problema bastante recurrente: tengo un algoritmo de stock que quiero operar en el espacio de registro (y solo usa un puñado de operaciones como '+ ',' - ',' ','/', etc.).

Esto me recuerda escribir funciones genéricas con plantillas, de modo que los argumentos de retorno, los parámetros, etc. se construyan del mismo tipo. Por ejemplo, si puedes hacer una FFT de flotadores, deberías poder hacer fácilmente uno en valores complejos (simplemente usando una clase que proporcione las operaciones necesarias para valores complejos).

Tal como está en la actualidad, parece que todas las implementaciones de FFT están escritas para obtener velocidad de borde de sangrado, por lo que no serán muy generales. Así que a partir de ahora, parece que tendría que volver a implementar FFT para los tipos genéricos ...

La razón por la que estoy haciendo esto es porque quiero mismas circunvoluciones de alta precisión (y el N^2 tiempo de ejecución es extremadamente lento). Cualquier consejo sería muy apreciado.

* Nota, podría necesitar implementar funciones trigonométricas para LogFloat, y eso estaría bien.

EDIT: Esto funciona porque LogFloat es un anillo conmutativo (y que no requiere la implementación de las funciones trigonométricas para LogFloat).La forma más sencilla de hacerlo fue volver a implementar FFT, pero @J.Sebastian también señaló una forma de utilizar la convolución genérica Python, que evita codificar la FFT (que, de nuevo, fue bastante fácil al usar un libro de texto DSP o the Wikipedia pseudocode).

Answer 1

Confieso que no me mantuve al día con las matemáticas en su pregunta. Sin embargo, parece que lo que realmente se está preguntando es cómo lidiar con números extremadamente pequeños y grandes (en valor absoluto) sin golpear subdesbordamientos y desbordamientos. A menos que te malinterprete, creo que esto es similar al problema que tengo trabajando con unidades de dinero, sin perder centavos en transacciones de miles de millones de dólares debido al redondeo. Si ese es el caso, mi solución ha sido el módulo decimal-matemático incorporado de Python. La documentación es válida tanto para Python 2 como para Python 3. La versión corta es que las matemáticas decimales son de precisión arbitraria y flotante y de punto fijo. Los módulos de Python cumplen con los estándares de IBM/IEEE para matemática decimal. En Python 3.3 (que actualmente está en formato alfa, pero lo he estado usando sin problemas), el módulo ha sido reescrito en C para una velocidad de hasta 100x (en mis pruebas rápidas).

Answer 2

Se podía escalar sus muestras en el dominio de tiempo por un número grande s para evitar flujo inferior, y luego, si

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

continuación

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

Ahora, utilizando el teorema de convolución se puede encontrar la manera de escalar el resultado final para eliminar el efecto del factor de escala.