Programmazione generica: Log FFT o convoluzione ad alta precisione (python)

Question

Ho un problema leggermente insolito, ma sto cercando di evitare la ricodifica FFT.

In generale, Voglio sapere questo: se io ho un algoritmo che viene implementato per il tipo float, ma sarebbe lavorare dove un certo insieme di operazioni è definito (numeri ad esempio complessi, per i quali anche definiscono +, *, ...), qual è il modo migliore per utilizzare tale algoritmo su un altro tipo che supporta tali operazioni? In pratica ciò è complicato perché gli algoritmi numerici generalmente sono scritti per velocità, non generalità.

In particolare:
sto lavorando con i valori con una gamma dinamica molto elevata, e così vorrei memorizzarli nello spazio di log (per lo più per evitare underflow).

Quello che mi piace è il registro della FFT di alcune serie:

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

Anche questo si tradurrà in underflow significativa. Quello che mi piacerebbe è quello di memorizzare i valori di registro e utilizzare + al posto di * e logaddexp al posto di +, ecc

Il mio pensiero su come fare questo è stato quello di implementare una semplice classe LogFloat che definisce queste operazioni primitive (ma opera nello spazio log). Quindi potrei semplicemente eseguire il codice FFT lasciandolo usare i miei valori registrati.

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 

Poi, l'idea sarebbe quella di costruire un elenco di LogFloat s, e quindi utilizzarlo nella FFT:

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

Questo può sicuramente essere fatto se re-implementare FFT (e semplicemente utilizzare LogFloat ovunque io preferisca usare float, ma ho pensato di chiedere un consiglio.Questo è un problema abbastanza ricorrente: ho un algoritmo di magazzino che voglio operare nello spazio di log (e usa solo una manciata di operazioni come '+ ',' - ',' ','/', ecc.).

Questo mi ricorda di scrivere funzioni generiche con modelli, in modo che gli argomenti di ritorno, i parametri, ecc. Siano costruiti dallo stesso tipo. Per exmaple, se si può fare un FFT di float, si dovrebbe essere in grado di fare facilmente uno su valori complessi (semplicemente usando una classe che fornisce le operazioni necessarie per valori complessi).

Allo stato attuale, sembra che tutte le implementazioni FFT siano scritte per la velocità di bleeding-edge, e quindi non saranno molto generali. Quindi, fin d'ora, sembra che avrei dovuto reimplementare FFT per i tipi generici ...

La ragione per cui sto facendo questo è perché voglio molto circonvoluzioni ad alta precisione (e la N^2 runtime è estremamente lento). Qualsiasi consiglio sarebbe molto apprezzato.

* Nota, potrebbe essere necessario implementare le funzioni trigonometriche per LogFloat e ciò andrebbe bene.

EDIT: Questo funziona perché LogFloat è un anello commutativo (e non richiede implementazione di funzioni trigonometriche per LogFloat).Il modo più semplice per farlo è stato reimplementare FFT, ma @ J.F.Sebastian ha anche indicato un modo di usare la convoluzione generica Python, che evita di codificare FFT (che, ancora una volta, era abbastanza semplice usando un libro di testo DSP o the Wikipedia pseudocode).

Answer 1

Confesso di non aver seguito la matematica nella tua domanda. Tuttavia, sembra proprio che quello che ti stai chiedendo sia come gestire numeri estremamente piccoli e grandi (in valore assoluto) senza colpire underflow e overflow. A meno che non ti fraintenda, penso che questo sia simile al problema che ho riscontrato con le unità di denaro, senza perdere i centesimi sulle transazioni da miliardi di dollari a causa degli arrotondamenti. Se questo è il caso, la mia soluzione è stata la mia soluzione risolta con il modulo matematico decimale di Python . La documentazione è valida sia per Python 2 sia per Python 3. La versione breve è che la matematica decimale è un tipo a virgola mobile e a virgola fissa di precisione arbitraria. I moduli Python sono conformi agli standard IBM/IEEE per la matematica decimale. In Python 3.3 (che è attualmente in forma alfa, ma l'ho usato senza problemi), il modulo è stato riscritto in C per una velocità fino a 100x (nei miei test rapidi).

Answer 2

Si potrebbe scalare i campioni nel dominio del tempo da un certo numero di grandi dimensioni s per evitare underflow, e poi, se

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

poi

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

Ora utilizzando il teorema di convoluzione si può capire come scalare il risultato finale per rimuovere l'effetto del fattore di scala.