Verfasst am: 18.10.2010, 11:01
Titel: Vektorisierung der Berechnung eines Beugungsmusters
Hallo ihr,
ich möchte ein Beugungsbild ausrechnen, dazu muss ein Integral über die geometrische Projektion des beugenden Körpers gerechnet werden. Ich habe es jetzt näherungsweise einfach durch eine Summe ersetzt und schon teilweise vektorisiert. Diese Phi-Vektorisierung hat schon einen Faktor ~100 in der Rechenzeit gebracht.
Kann man die for-Schleife über die Theta auch noch vektorisieren oder irgendwie anders schneller machen?
for iTheta=1:nTheta
[xObs,yObs,zObs] = sph2cart(thetaVal(iTheta),phiVal+pi/2,1);
D(:,iTheta) = D(:,iTheta) + exp(-1i*k*scale*[x y 0]*[xObs; yObs; zObs])';
end end end
k, scale etc. sind Konstanten des Problems, in thetaVal und phiVal stehen die Raumrichtungen, unter denen die Beugungsfigur berechnet wird.
Im Moment rechne ich das für eine kreisförmige Blende/Kugel, später habe ich als Eingangsdaten nur die (x,y)-Koordinaten von Punkten, an denen ich den Körper "getroffen" habe, die also einen Beitrag leisten.
Würde mich über Anregungen oder andere Ideen freuen.
Verfasst am: 18.10.2010, 18:42
Titel: Re: Vektorisierung der Berechnung eines Beugungsmusters
Hallo bligg,
Zitat:
Code:
for iTheta=1:nTheta
[xObs,yObs,zObs] = sph2cart(thetaVal(iTheta),phiVal+pi/2,1);
D(:,iTheta) = D(:,iTheta) + exp(-1i*k*scale*[x y 0]*[xObs; yObs; zObs])';
end
Nun ist "r" bei Dir immer 1, und "elev" ist auch konstant, und damit auch der Z-Wert. Damit kannst Du die Berechnung auch deutlich vereinfachen. Ausserdem kannst Du SPH2CART auch gleich für einen ganzen Vektor ausrechnen lassen.
Es sieht so aus, als sei D pre-allociert, wohl mit ZEROS, oder?
Ich bin davon ausgegangen, dass sph2cart built-in sein würde, und deshalb gar nicht erst auf die Idee gekommen, mal in den Code zu schauen (obwohl natürlich klar ist, was dort passiert).
Insofern war deine Anregung mit den konstanten Werten, die ich jedes Mal auf's Neue ausrechne, natürlich absolut richtig.
Für's Protokoll und eventuell später drauf Stoßende, das Problem lässt sich mittels
for iTheta = 1:nTheta
for iPhi = 1:nPhi
[xObs,yObs,zObs] = sph2cart(thetaVal(iTheta),phiVal(iPhi)+pi/2,1);
cartVal(iPhi,iTheta,:) = [xObs yObs zObs];
end end
Allerdings wird ein großer Teil der Rechenzeit (plausibelerweise) in der Berechnung der e-Funktion verbraucht. Viel hat das also leider nicht mehr gebracht.
Dafür hat man in dieser Form mit bsxfun ziemlich genau lineare Laufzeitabhängigkeit von nRay (klar), aber auch der Anzahl der Beobachtungsstellen (phiVal, thetaVal).
Leider ist die Laufzeit mit ~60s / 10^5 Strahlen immer noch einige Größenordnungen zu teuer für mich, so dass ich mich dennoch auf die Suche nach einem anderen Ansatz zur Berechnung machen werde.
wenn Du eine vollständige lauffähige Version posten würdest, könnten wir damit experimentieren.
Es sieht so aus, als habest Du den Aufruf von SPH2CART nicht beschleunigt, oder? In den Schleifen über iTheta und iPhi könnte man also noch viel Zeit sparen. Mit BSXFUN lassen sich sogar beide Schleifen ganz weglassen.
Aber die Schleife über nRay ist wohl aufwändiger.
Zunächst mal kannst Du die Berechnung von "x^2+y^2<=1" aus der Schleife heraus ziehen und so unnötige Iterationen sparen.
Dann multiplizierst Du "cartVal" mit "[x y 0]". Das Multiplizieren mit 0 kann man sich aber sparen.
Wahrscheinlich ist
schneller als der BSXFUN Aufruf. (Kommt das gleiche Ergebnis heraus??)
Ich würde noch "-1i*k*scale" als Konstante definieren.
Ich vermute, das Einführen der 3. Dimension für BSXFUN ist nicht nötig, aber das würde ich besser noch mal auspobieren - deswegen meinte ich, es wäre gut eine lauffähige Version zu posten.
ich will das Beugungsmuster eigentlich als Modul eines größeren Programms benutzen (wie gesagt, lochförmige Blende ist nur zur Validierung), deshalb funktioniert das Abfrage rausziehen zur Beschleunigung nicht.
Ja, die sph2cart habe ich nicht beschleunigt, weil die Schleifen tatsächlich nur einen kleinen Teil ausmachen, und nur einmal bei der Initialisierung anfallen. Insbesondere sind sie von der Anzahl der Testrays unabhängig.
Das Ausrechnen der Konstante ist natürlich prinzipiell eine komfortable Möglichkeit, Performance zu gewinnen, ich suche aber Rechenzeitersparnis in der Größenordnung 100 oder mehr.
Ich denke, da muss man diesen Ansatz einfach als gescheitert bezeichnen.
Es geht darum, dass ich später auf eine beliebig komplizierte Geometrie (die Kollisionsabfrage und Konsorten stehen) diese Strahlen haue, um die Brechung zu berechnen. Ich interessiere mich insbesondere nur für die über (hier) Theta gemittelten Werte.
(Kleiner Einschub: Ist es im englischsprachigen Raum üblich, die Kugelkoordinaten mit gegenüber der mir im deutschsprachigen Raum bekannten Winkelnamenkonvention verdrehten Namen zu bezeichnen?)
Etwas klarer: mir reicht die über den Azimut gemittelte Beugungsfigur in Abhängigkeit der Elevation/Polarwinkel.
Prinzipiell muss man dafür ein Integral auswerten, das wie eine Fourier-Transformation aussieht (Fourier-Optik), allerdings ist mir bei diesem Ansatz nicht klar, wie ich dort meine Winkel wieder rausfinde. Es gibt ein Beispiel von TMW: http://www.mathworks.com/help/techdoc/math/brentm1-1.html#brfb2uq
Ab dem Bild fehlt mir die Idee, wie ich an die Winkel komme.
Noch ein anderer Ansatz wäre, das Oberflächenintegral über die Projektion des Körpers in ein äquivalentes Konturintegral umzuwandeln (Satz von Stokes), das man dann wahrscheinlich noch ein Stück weit analytisch lösen kann.
Dort stellt sich aber das Problem, dass ich die Kontur irgendwie bestimmen müsste. Das ist ohne weitere Einschränkungen nur für konvexe Hüllen möglich.
Ich bekomme aus meiner Kollisionsabfrage sozusagen entweder an zufälligen Stellen oder in einem regelmäßigen Grid eine "getroffen" oder "nicht getroffen"-Ansage. Ich bekomme also eine Punktwolke, so oder so, um die ich eine Hülle bestimmen müsste, um das Konturintegral ausrechnen zu können.
Das nur zur Einordnung allgemein, ich hänge noch den etwas unaufgeräumten, weil so erstmal nur zum Testen geschriebenen, Code an. Es ist aber meines Erachtens verschenkte Liebesmüh', dort nach den letzten Optimierungsmöglichkeiten zu suchen; dieser stumpfe Zugang, das Oberflächenintegral direkt auszurechnen, scheint mir dann in Hinblick auf computational aspects ungeeignet.
for iTheta = 1:nTheta
for iPhi = 1:nPhi
[xObs,yObs,zObs] = sph2cart(thetaVal(iTheta),phiVal(iPhi)+pi/2,1);
cartVal(iPhi,iTheta,:) = [xObs yObs zObs];
end end
Die Aufteilung und das einzelne Auslesen von x und y auch nur, weil das die Form ist, in der ich später die Punkte eingeben würde, da bin ich mir der Unsauberkeiten bewusst!
Liebe Grüße
bligg
Einstellungen und Berechtigungen
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS