Entwickler-Ecke

Algorithmen, Optimierung und Assembler - Ableiten mittels FFT


Christian S. - Fr 23.11.07 16:50
Titel: Ableiten mittels FFT
Hallo!

Aktuell versuche ich eine Simulation zu schreiben, bei der zwei gekoppelte PDEs numerisch gelöst werden. Dabei möchte ich die Ortsableitung mittels FFT machen, da man dann die Diskretisierung nicht so fein machen muss.

Im Prinzip ist's ja simpel:
Array füllen ("a"),
Hintransformieren,
jedes Element a[j] mit i*j*dk (i = sqrt(-1)) multiplizieren,
wieder rücktransformieren.

Ich teste das mit einem Sinus. Ich pack einen Sinus ins Array rein und durchlaufe obige Prozedur. Der Realteil dessen was herauskommt ist auch richtig (bis auf Skalierung, aber das ist erstmal wurscht), nur hat's auch einen Imaginärteil. Und das sollte ja nicht so sein, denn eine reelle Funktion abgelitten ergibt ja immer noch eine relle Funktion :gruebel:

Habe ich da irgendwo einen prinzipiellen Fehler? Hab irgendwie gerade ein Brett vor'm Kopf ...

Grüße
Christian


Moderiert von user profile iconChristian S.: Topic aus Off Topic verschoben am Fr 23.11.2007 um 16:22


delfiphan - Fr 23.11.07 17:39
user defined imageuser defined imageuser defined image

Dies gilt für die Fourier-Transformation. So direkt kann man das aber wahrscheinlich nicht auf die DFT übertragen, da ein periodisches Signal angenommen wird und die Frequenz nicht entsprechend linear bis ins Unendliche steigen kann. Ich habe jedoch schon Beispiele gesehen/gerechnet, wo mittels DFT eine Ableitung berechnet wurde, jedoch habe ich das etwas komplizierter in Erinnerung. Es war iirc möglich, eine (vollbesetzte) "Ableitungsmatrix" herzuleiten bzw. zu berechnen, womit man diskretisierte Funktionen durch Multiplikation ableiten konnte.

Könnte auch sein, dass du etwas falsch implementiert hast? Z.B. das Array falsch adressiert (Ursprung der Frequenzachse am Anfang des Arrays, ab der Hälfte des Arrays kommen dann die negativen Frequenzen, ... kann aber auch anders aussehen).


Christian S. - Mo 26.11.07 15:10
Oh, Du hattest ja nochmal editiert. Gar nicht gesehen :lupe: Auf jeden Fall danke für Deine Mühe! :zustimm:

Der Fehler war folgender:
Die FFT sortiert einem ja die Stützstellen um (bildlich gesprochen), sodass man im Fourierraum nicht einfach mit i*j*dk multiplizieren kann. Man muss da ein bisschen umsortieren. Die k-Werte, mit denen man multiplizieren muss, bekommt man so:


Chrome-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
 
class method FFT.FourierIndex(var k : MyArray<Double>; dx : Double; k0 : Double); 
begin
  var dk := 2 * Math.PI / dx / k.Length;
  if k.Length = 0 then
    raise new Exception('Array length must not be zero!');
    
  var N_half := k.Length div 2;
  
  if N_half = 0 then
    k[0] := 0
  else
  begin
    for ik : Integer := k.FirstIdx to k.FirstIdx + N_half - 1 do
      k[ik] := dk  * (ik - k.FirstIdx) + k0;
      
    for ik : Integer := k.FirstIdx + N_half to k.LastIdx do
      k[ik] := dk  * (ik - k.Length - k.FirstIdx) + k0;  
  end;  
end;

Dabei ist MyArray<Double> ein Array of Double, nur dass man den Anfangsindex nicht auf 0 haben muss.


delfiphan - Mo 26.11.07 22:59
Ja sowas in der Art hatte ich auch vermutet (falsch adressiert/Klammerbemerkung am Schluss). Dann noch viel Spass mit deinen PDEs und hoffentlich explodieren sie dir nicht.


Entwickler-Ecke.de based on phpBB
Copyright 2002 - 2011 by Tino Teuber, Copyright 2011 - 2025 by Christian Stelzmann Alle Rechte vorbehalten.
Alle Beiträge stammen von dritten Personen und dürfen geltendes Recht nicht verletzen.
Entwickler-Ecke und die zugehörigen Webseiten distanzieren sich ausdrücklich von Fremdinhalten jeglicher Art!