Komplexe Funktionen invers Fast Fourier Transformieren

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
Fredi22
User
Beiträge: 10
Registriert: Montag 7. November 2016, 16:28

Mit großer Freude bin ich auf dieses Forum hier gestoßen :).

Ich möchte zwei ziemlich umfangreiche komplexe Funktionen Invers Fast Fourier Transformieren und habe dazu folgendes Mini Skript gefunden:

Python: Recursive

Code: Alles auswählen

from cmath import exp, pi
 
def fft(x):
    N = len(x)
    if N <= 1: return x
    even = fft(x[0::2])
    odd =  fft(x[1::2])
    T= [exp(-2j*pi*k/N)*odd[k] for k in range(N//2)]
    return [even[k] + T[k] for k in range(N//2)] + \
           [even[k] - T[k] for k in range(N//2)]
 
print( ' '.join("%5.3f" % abs(f) 
                for f in fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])) )

Hierzu habe ich mehrere Fragen:

1. Wie mache ich aus der FFT eine inverse FFT?
2. Wo würde ich meine komplexe Funktion eingeben? Müsste ich einfach nach dem Skript etwas eingeben wie fft(f(x)) oder wie würde das ganze dann aussehen?

Vielen Dank schonmal im vorraus!
Zuletzt geändert von Anonymous am Montag 7. November 2016, 16:44, insgesamt 1-mal geändert.
Grund: Quelltext in Python-Codebox-Tags gesetzt.
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

Schau Dir diese Seite mal an:
http://www.magben.de/?h1=mathematik_fue ... urier_code
dort habe ich in einem Beispiel Python Code zu
  • Fourier-Transformation,
  • Fourier-Reihe und
  • FFT
Die transformierte Funktion ist in dem Beispiel reell, der Code sollte aber auch mit einer komplexen Funktion funktionieren.
Bild
a fool with a tool is still a fool, www.magben.de, YouTube
Fredi22
User
Beiträge: 10
Registriert: Montag 7. November 2016, 16:28

Danke für die Antwort. Die verlinkte Seite ist jedoch leider wenig hilfreich.

Ich verwende nun die numpy funktion numpy.fft.ifft().
Auf der Seite https://docs.scipy.org/doc/numpy/refere ... y.fft.ifft wird die funktion kurz erklärt.

fft.ifft(a,n,axis norm). In der () sind also die parameter.

Bei der FFT bzw.IFFT geht es darum das Integral über -undendlich bis unendlich von (e^(-i*k*w) * f(w)) mit simpson weigtings numerisch anzunähern, wobei das ergebnis die Inverse fast fourier transformierte der Funktion f(w) ist, nennen wir sie F(k). Hierfür wird der definitionsbereich von w mithilfe von n-vielen w´s diskretisiert, welches entsprechend n diskrete werte für f(w) liefert.

In der verlinkten Beschreibung wird gesagt, dass der parameter "a" der Input und ein array ist. Ich gehe also mal davon aus es handelt sich um die diskreten Werte von f(w). Das array hat dann natürlich die n viele elemente. Die funktion fft.ifft funktioniert mit nur einem eintrag in dem klammern, nämlich dem input array a. numpy.fft.ifft(a) funktioniert also und gibt mir auch n komplexe diskrete Werte aus. Allerdings muss ich ja noch die schrittweite von w, nennen wir es dw eingeben, da es die ristriktion gibt bei FFT, dass

dw*dk = 2pi/n ergibt.

D.h. mann muss normalerweise immer irgendwo die Schrittweite dk bzw. dw angeben. Ich frage mich ob der parameter axis damit zusammenhängt. Ich verstehe allerdings nicht was axis und norm bedeuten sollen.
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

Wenn Du eine Funktion f(x) im Orts-Raum hast, die Du mit der FT in den k-Raum zu F(k) transformieren willst, dann brauchst Du dazu numpy.fft.fft.

Die Parameter von numpy.fft.fft:
a: das zu transformierende Numpy-Array. a muss entweder periodisch sein oder bei 0 anfangen und bei 0 aufhören.
n: das Array a wird auf die Länge n vergrößert und mit Nullen aufgefüllt. Das macht dann Sinn, wenn das Array sehr lang (>1e5) ist und Dir numpy.fft.fft zu lang dauert, dann kannst Du fft beschleunigen, wenn n eine Zweierpotenz ist, also wenn das Array 750000 Elemente hat, dann übergibst Du n=2**20.
axis: Diesen Parameter brauchst Du nur bei mehrdimensionalen Arrays, wenn Du z.B. die Frequenzanalyse von mehreren Signalen gleichzeitig machen möchtest.
norm: Habe ich noch nie benutzt.

f sei nun ein Array mit den Funktionswerten und x enthält die äquidistanten Ortskoordinaten dazu, dann ist

dx=x[1]-x[0]
k = 2 * numpy.pi * numpy.fft.fftfreq(x.size, dx)
F = numpy.fft.fft(f) * dx
a fool with a tool is still a fool, www.magben.de, YouTube
Fredi22
User
Beiträge: 10
Registriert: Montag 7. November 2016, 16:28

Also ich Sache ist, es handelt sich bei meinem problem um eine Call Options preis der Fourier transformiert ist und nun von mir zum Call Preis wieder zurück transformiert werden soll. der Call Preis ist im K Space und der Fourier Transformierte Preis im angular frequency space w.

ich habe also f(w) dass ich INVERS zurück transformieren muss zu C(k).

Wie auch immer ich bekomme keine anständige FFT bzw. IFFT hin. Ich poste mal ein simples Beispiel, was nicht mit meinen funktionen zu tun hat.

Code: Alles auswählen

import math
import numpy as np
import matplotlib.pylab as plt 
from math import *
import time
from scipy.fftpack import fft, ifft

N=32

l=[]
l2 =[]
for i in range(0,N):
    x=2*i
    y= np.sin(2*x)
    l.append(y)
    l2.append(x)

dx = l2[1]-l2[0]    
dw = 2*np.pi / (dx*N)

a=np.fft.fft(l) * dx 

b=np.fft.ifft(a) * dw


print (l)
print (a)
print (b)
    
Output
l=

[0.0, -0.7568024953079282, 0.98935824662338179, -0.53657291800043494, -0.2879033166650653, 0.91294525072762767, -0.90557836200662378, 0.27090578830786904, 0.55142668124169059, -0.99177885344311578, 0.74511316047934883, 0.017701925105413577, -0.76825466132366682, 0.98662759204048534, -0.52155100208691185, -0.30481062110221668, 0.92002603819679074, -0.8979276806892913, 0.25382336276203626, 0.56610763689818033, -0.99388865392337522, 0.73319032007329221, 0.03539830273366068, -0.77946606961580467, 0.98358774543434491, -0.50636564110975879, -0.32162240316253093, 0.92681850541778499, -0.88999560436683334, 0.23666139336428604, 0.58061118421231428, -0.99568698688917945]
a=fft(l)
>> >[-1.49580427+0.j -1.50491260-0.08606522j -1.53318920-0.17491292j
..., -1.58372464+0.26973505j -1.53318920+0.17491292j
-1.50491260+0.08606522j]
b=ifft(a)
[ -2.17991781e-17 +0.00000000e+00j -1.48597822e-01 -5.43504982e-17j
1.94260037e-01 -7.90220207e-17j ..., 4.64683559e-02 +4.02058386e-17j
1.14002739e-01 +5.44979453e-18j -1.95502683e-01 +4.61758064e-17j]

man sieht denke ich ziemlich deutlich das nur murks raus kommt, da B ja wieder l seien müsste. Ich habe jetzt sin(2x) verwendet weil es periodisch ist. Meine Call preis funktionen sind nicht periodisch, ich bin mir aber sehr sicher dass man FT IFT auch auf nicht periodische funktionen anwenden kann.
Fredi22
User
Beiträge: 10
Registriert: Montag 7. November 2016, 16:28

Die scipy version scheint zu laufen:

Code: Alles auswählen

import math
import numpy as np
import matplotlib.pylab as plt
from math import *
import time
from scipy.fftpack import fft, ifft

N=32

l=[]
l2 =[]
for i in range(0,N):
    x=-N/2*2+2*i
    l.append(x)
    
al=np.array(l)    
y= np.sin(2*al)
yfft=fft(y)
bl=ifft(yfft)  

print (y)
print (yfft)    
print (bl)
Output
original funktion
[-0.92002604 0.30481062 0.521551 ..., 0.98662759 -0.521551 -0.30481062]

fft
>>> [-0.92002604+0.j -0.92002604-0.10981706j -0.92002604-0.2231845j
..., -0.92002604+0.34417516j -0.92002604+0.2231845j
-0.92002604+0.10981706j]

ifft(fft)=wieder original funktion
[-0.92002604 +0.00000000e+00j 0.30481062 -7.47756040e-17j
0.52155100 +1.04083409e-17j ..., 0.98662759 -8.43884260e-17j
-0.52155100 -2.42861287e-17j -0.30481062 +1.99675694e-16j]
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

Probier das mal, dann siehst Du wie "periodisch" Deine zu transformierende Funktion ist;

Code: Alles auswählen

plt.figure()
plt.plot(l2, l)
plt.show()
Die Sinus-Funktion ist nicht besonders gut geeignet um damit FT und FFT zu üben. Die Fourier-Transformierte der Sinus-Funktion ist die Delta-Funktion und die Delta-Funktion lässt sich nicht mal eben so mit Python implementieren. Such Dir zum Üben ein Funktionspaar ohne Delta-Funktion aus, nimm lieber was mit der Exponential-Funktion:
https://de.wikipedia.org/wiki/Fourier-T ... ions-Paare

FT und FFT ist auch eher ein fortgeschrittenes Numpy-Thema und wenn ich Deinen Code sehe, da fehlen erstmal die Grundlagen von Numpy.
So erzeugst Du periodische und ausreichend abgetastete Sinus-Werte:

Code: Alles auswählen

x = np.linspace(0, 1, 1000, endpoint=False)
y = np.sin(2*np.pi * x)
plt.figure()
plt.plot(x,y)
plt.show()
a fool with a tool is still a fool, www.magben.de, YouTube
Fredi22
User
Beiträge: 10
Registriert: Montag 7. November 2016, 16:28

Wie auch immer: Das wichtige für mich ist zu verstehen was die fft module genau machen. Ich brauche ein Hülle welche INTRGERAL [e**(-ikw) f(w) dw] approximiert. Unabhängig davon was f(w) ist. Ich kann dies mit IFT machen oder mit IFFT approximieren. Nach meinem Verständnis ist IFFT die numerische Lösung des angegeben Integral mit Simpson weighting und genau da brauche ich. Scipy scheint unter fft was anderes zu verstehen.

Code: Alles auswählen

import math
import numpy as np
import matplotlib.pylab as plt 
from math import *
import time
from scipy.fftpack import fft, ifft
start_time = time.clock()

# Definition od the option variables

r = -0.0032
S = 64.97

# ln (S)
s= np.log(S)

# Maturity
T = 81/365

vw =[]
# Number of Inputs. Equals the number of samples in the angular frequency space. Must be a power of 2 for FFT.
N = 4096

# Setting the angular frequency interval (spacing size)
dw = 0.25

# dw implies the following intervall dk in the log-strike space 
dk = 2*np.pi / (N*dw)

logK = []
K = []

# Set dumping parameter
a = 1

# Set Model parameter ----------------------------------------------------------------------------------------------------------- 

# sigma
vola = 0.2725523671

# -------------------------------------------------------------------------------------------------------------------------------

# phi, characteristic function (CF) of the chosen Lévy process. TO APPLY FFT FOR OTHER lÉVY PROCESSES JUST CHANGE THE CF (phi) below
def phi(w):
    return np.exp(1j*(s+(r-0.5*vola**2)*T)*w-(vola**2*T*w**2*0.5))

    # Fills the input array with N discrete values of the CF
for i in range(0, N):
 # discrete angular frequencies values w     
 wi = i*dw  
 vw.append(wi) 
 logstrikes = -(N/2) * dk + i*dk  
 logK.append(logstrikes)
 strike = np.exp(logK[i])
 K.append(strike)

w = np.array(vw)
 
g = (np.exp(-r*T)*phi(w-(a+1)*1j))/((a**2)+a-w**2+(1j*(2*a+1)*w))
fg = ifft(g)

alogK = np.array(logK)


C= (np.exp(-a*alogK)/(2*np.pi)) * (fg)


plt.plot(K[2048:2800], (C[2048:2800]))

#plt.plot(K[:2900], logK[:2900])
# calculates the computation time   
time = time.clock() - start_time 



Wäre die definition die gewünschte, müssten reale werte herauskommen bei der ifft. Dies wären die Call Preise, da "g" definitiv(!) die FT der Callpreise sind, im sinne der Definition des obigen integrals. Also g(w)= INTEGRAL (e(ikw) C(k) dk). Wenn ich verstehen würde was scipy ifft genau macht könnte ich da sproblem auch selber lösen.
Fredi22
User
Beiträge: 10
Registriert: Montag 7. November 2016, 16:28

MagBen hat geschrieben:Probier das mal, dann siehst Du wie "periodisch" Deine zu transformierende Funktion ist;

Code: Alles auswählen

plt.figure()
plt.plot(l2, l)
plt.show()
Die Sinus-Funktion ist nicht besonders gut geeignet um damit FT und FFT zu üben. Die Fourier-Transformierte der Sinus-Funktion ist die Delta-Funktion und die Delta-Funktion lässt sich nicht mal eben so mit Python implementieren. Such Dir zum Üben ein Funktionspaar ohne Delta-Funktion aus, nimm lieber was mit der Exponential-Funktion:
https://de.wikipedia.org/wiki/Fourier-T ... ions-Paare

FT und FFT ist auch eher ein fortgeschrittenes Numpy-Thema und wenn ich Deinen Code sehe, da fehlen erstmal die Grundlagen von Numpy.
So erzeugst Du periodische und ausreichend abgetastete Sinus-Werte:

Code: Alles auswählen

x = np.linspace(0, 1, 1000, endpoint=False)
y = np.sin(2*np.pi * x)
plt.figure()
plt.plot(x,y)
plt.show()
Sry aber dann ist numpy nicht der richtige ansatz. Scipy kann problems los sin(x) fftn und dann zurrück transformieren. Ich werde von nun an nur noch scipy benutzen. Dein Code unterscheided sich nur von meinem durch das engere grid bei der diskretisierung. Mein Problem ist aber unabhängig vom grid, dass wurde bereits in diversen papers getestet und sollte für jedes grid funktionieren, da es nicht um Signale geht sondern lediglich um eine mathematische Ausführung des Integrals. Diese ausführung MUSS aber durch FFT geschehen. Zumindest muss ich es am ende so nennen können. Die Ausführung des von mir angegebeben Integrals mit simpson weightings ist ausreichend um es am ende FFT nennen zu können. Trotzdem danke für die Antwort.
BlackJack

@Fredi22: Wenn Dir die Numpy-Grundlagen fehlen, dann hilft auch Scipy nicht weiter, denn das baut ja auf Numpy auf.
Fredi22
User
Beiträge: 10
Registriert: Montag 7. November 2016, 16:28

Ich hab das problem gelöst... die scipy fft funktioniert und ich bekomm die richtige Werte raus.
Abgesehen davon fehlen mir keine numpy basics, zumindest nicht um sowas Einfaches wie mein Beispiel zu lösen. Ich hab lediglich das array über ne schleife gefüllt und er über linspace, daraus hat er dann geschlossen dass mir basics fehlen ;-). p.s. man kann sinus sehr wohl über python transformieren, dafür gibs eine million Beispiele im Netz und ich habs eben selber noch gemacht. Python wär ja ne tolle Sprache wenn sowas einfaches nicht ginge :).
Antworten