next up previous index
Next: Praktische Umsetzung als Excel-Funktion Up: RechnenmitSpektren Previous: Addition zweier Spektren   Index


Multiplikation zweier Spektren

Zwei Spektren durchzumultiplizieren bedeutet ihre Funktionswerte Punkt für Punkt zu multiplizieren. Das Ergebnis ist wieder ein Spektrum, welches sich jedoch in den Werten und dem Gültigkeitsbereich von den Ausgangsspektren unterscheiden kann. Gehen wir davon aus, dass uns die Spektralwerte als Stützstellen in bereits sortierten Listen vorliegen. Zwei durchzumultiplizierende Spektren werden mit zwei Variablen $i$ und $j$ durchlaufen. Abbildung [*] zeigt den prinzipiellen Ablauf und die Probleme beim Durchmultiplizieren:
Figure: Spektren
\includegraphics[width=0.60\textwidth]{HTMLBilder/FaltungPlot.pdf}
  1. $x_{i}$ und $x_{j}$ liegen nicht immer an gleichen Stellen. Somit lassen sich die $y$-Werte nicht direkt einander zuordnen. Man kann versuchen direkt verwertbare Werte zu interpolieren (hier Rot).
  2. Spektrum $i$ fängt erst an, wo Spektrum $j$ bereits Werte angenommen hat. Der Gültigkeitsbereich kann also erst ab $x_{i}$ beginnen. Ebenso verhält es sich damit, dass die Stelle $x_{j+4}$ über den Bereich der $i$-Werte ragt.
Durch Interpolation auf die roten Werte erhält jede Stützstelle eine Partnerstelle an gleicher $x$-Position um die Multiplikation durchzuführen. Das Ergebnis ist im Beispiel die graue Kurve. Sie erstreckt sich über den Bereich, in dem sich beide Spektren überlappen.

Ein Algorithmus muß beide Spektren abhängig vom lokalen Umfeld durchlaufen um die richtigen Stützstellen $i$ , $i+1$, $j$ und $j+1$ zur Interpolation parat zu haben und zu entscheiden, wie wir zum nächsten Schritt weiterrücken, indem wir entweder $i$ oder $j$ oder Beide hochzählen. Dabei dürfen $i$ und $j$ Werte zwischen dem Anfang 0 und ihrem Maximalwert minus 1 annehmen.

$\displaystyle 0 \leq i \leq i_{\text{max}}-1 \qquad 0 \leq j \leq j_{\text{max}}-1 $

Die letzte Ergebnisstützstelle müssen wir separat berechnen, da die letzte Stützstelle $i_{\text{max}}$ oder $j_{\text{max}}$ keine nächste Stelle mehr hat.

Abbildung [*] zeigt die Situation, bei der die aktuellen Zählerpositionen zu weit auseinanderliegen.

Figure: Kein Überlapp in lokalem Umfeld wenn $x_{i+1}\leq x_{j}$ oder $x_{j+1}\leq x_{i}$.
\includegraphics[width=0.90\textwidth]{HTMLBilder/Heranruecken.pdf}
Durch sukzessives Erhöhung von $i$ bzw. $j$ können wir den niedrgeren Bereich heranrücken. Liegt $x_{i+1}\leq x_{j}$, dann erhöhen wir $i$ um Eins. Liegt $x_{j+1}\leq x_{i}$, dann erhöhen wir $j$. Das geschieht so lange, bis wir einen Überlapp der beiden lokalen Bereiche erreicht haben. Wir erfassen alle Szenarien separiert nach den Randbedingungen In Abbildung [*] unterscheiden wir weiter nach der relativen Lage von $x_{j+1}$. Sie bestimmt den Weg zum nachfolgenden Schritt.
Figure: $x_{j}<x_{i}<x_{j+1}$
\includegraphics[width=0.90\textwidth]{HTMLBilder/XigroesserXj.pdf}
Wir können jetzt den Wert an der Stelle $i$ multiplizieren mit dem interpolierten Wert zwischen den Stellen $j$ und $j+1$ des zweiten Spektrums. Die gestrichelte rote Linie mit dem dicken roten Punkt für den interpolierten Wert steht wie in Abbildung [*] für diesen Vorgang. Das Wertepaar $[ x , y ]$ des durchmultiplizierten Spektrums ist dann

$\displaystyle [ x_{i} , ( y_j + \frac{y_{j+1}-y_{j}}{x_{j+1}-x_{j}} (x_{i}-x_{j}) ) \cdot y_{j} ] $

und würde dem grauen Punkt in Abbildung [*] entsprechen. Die gestrichelte rote Linie mit dem kleinen roten Punkt zeigt jeweils die Vorschau auf den nächsten Schritt.
Figure: $x_{i}=x_{j}< \{x_{i+1} \; ; \; x_{j+1}\}$
\includegraphics[width=0.90\textwidth]{HTMLBilder/XigleichXj.pdf}
Im Szenario von Abbildung [*] können wir den Wert an der Stelle $i$ multiplizieren mit dem Wert an Stelle $j$ des zweiten Spektrums. Die gestrichelte rote Linie mit dem dicken roten Punkt steht für diesen Vorgang. Das Wertepaar $[ x , y ]$ des durchmultiplizierten Spektrums ist dann

$\displaystyle [ x_{i} , y_i \cdot y_{j} ] $

und würde dem grauen Punkt in Abbildung [*] entsprechen. Die gestrichelte rote Linie mit dem kleinen roten Punkt zeigt die Vorschau auf den nächsten Schritt.
Figure: $x_{i}<x_{j}<x_{i+1}$
\includegraphics[width=0.90\textwidth]{HTMLBilder/XikleinerXj.pdf}
Im Szenario von Abbildung [*] können wir den Wert an der Stelle $i$ multiplizieren mit dem interpolierten Wert zwischen den Stellen $j$ und $j+1$ des zweiten Spektrums. Die gestrichelte rote Linie mit dem dicken roten Punkt für den interpolierten Wert steht wie in Abbildung [*] für diesen Vorgang. Das Wertepaar $[ x , y ]$ des durchmultiplizierten Spektrums ist dann

$\displaystyle [ x_{j} , ( y_j + \frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}} (x_{j}-x_{i}) ) \cdot y_{j} ] $

und würde dem grauen Punkt in Abbildung [*] entsprechen. Die gestrichelte rote Linie mit dem kleinen roten Punkt zeigt die Vorschau auf den nächsten Schritt.

Haben wir uns durch Erhöhung von $i$ und $j$ immer weiter durch die Ausgangsspektren gearbeitet, kommen wir irgendwann an den Endpunkt, an dem $i = i_{\text{max.}}$ oder $j= j_{\text{max.}}$ ist. Es gibt an dieser Stelle keine Vorschau über diesen Punkt hinaus und wir müssen den letzten Punkt von der rückwärtigen Seite interpolieren. Die rückwärtige Seite $x_{i-1}$ oder $x_{j-1}$ gibt es natürlich nur, wenn $i$ oder $j$ bereits einmal hochgezählt wurden, also größer Null sind. Das ist nicht immer der Fall, beispielsweise wenn ein Spektrum mit 2 Stützstellen alle Stützstellen des anderen Spektrums umschließt. In den Fällen, bei denen $i$ oder $j$ Null sind, können wir jedoch immer getrost um Eins hochzählen, denn zwei Stützstellen gibt es immer. Sonst ist es kein Bereich und es gibt kein Integral. Für den Fall dass $i$ den Endpunkt erreicht hat und $j>0$ ist wäre die letzte Stützstelle dann

$\displaystyle [ x_{i} ; ( y_{j-1} + \frac{y_{j}-y_{j-1}}{x_{j}-x_{j-1}} (x_{i}-x_{j-1}) ) \cdot y_{i} ] $

und für den Fall dass $j$ den Endpunkt erreicht hat und $i>0$ ist wäre die letzte Stützstelle

$\displaystyle [ x_{j} ; ( y_{i-1} + \frac{y_{i}-y_{i-1}}{x_{i}-x_{i-1}} (x_{j}-x_{i-1}) ) \cdot y_{j} ] $

und damit sind wir fertig. Die bisherige Vorgehensweise ist in Tabelle [*] zusammengefasst.

Table: Zusammenfassung
Bedingung $[ x ; y ]$ Inkrement  
$x_{i}\geq x_{j+1}$   $j++$  
$x_{j}\geq x_{i+1}$   $i++$  
$x_{i}=x_{j}$ $[ x_{i} ; y_i \cdot y_{j} ] $ $x_{i+1}=x_{j+1} \Rightarrow i++ ; j++$  
$x_{i+1}<x_{j+1} \Rightarrow i++$  
$x_{i+1}>x_{j+1} \Rightarrow j++$  
$x_{i}<x_{j}$ $[ x_{i} ; ( y_i + \frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}} (x_{j}-x_{i}) ) \cdot y_{j} ] $ $x_{i+1}=x_{j+1} \Rightarrow i++ ; j++$  
$x_{i+1}<x_{j+1} \Rightarrow i++$  
$x_{i+1}>x_{j+1} \Rightarrow j++$  
$x_{i}>x_{j}$ $[ x_{j} ; ( y_j + \frac{y_{j+1}-y_{j}}{x_{j+1}-x_{j}} (x_{i}-x_{j}) ) \cdot y_{i} ] $ $x_{i+1}=x_{j+1} \Rightarrow i++ ; j++$  
$x_{i+1}<x_{j+1} \Rightarrow i++$  
$x_{i+1}>x_{j+1} \Rightarrow j++$  
$i = i_{\text{max.}}$ oder $j= j_{\text{max.}}$ $[ x_{i} ; ( y_{j-1} + \frac{y_{j}-y_{j-1}}{x_{j}-x_{j-1}} (x_{i}-x_{j-1}) ) \cdot y_{i} ] $    


In die Programmiersprache Javascript umgesetzt wäre die in Listing [*] umgesetzte Funktion eine Möglichkeit.
Listing: Multiplikation zweier Spektren in Javascript 
    1 
    2 function durchhangelnmul(spek1, spek2){
    3     var s=0,i=0,j=0, arr=[];
    4     while (i < spek1.length-1 && j< spek2.length-1) {
    5         if(spek1[i]["x"]>=spek2[j+1]["x"]){j++; continue;};
    6         if(spek2[j]["x"]>=spek1[i+1]["x"]){i++; continue;};
    7 
    8         if(spek1[i]["x"]==spek2[j]["x"]){
    9             arr.push({x:spek1[i]["x"], y: spek1[i]["y"]*spek2[j]["y"]});
   10         };
   11         if(spek1[i]["x"]<spek2[j]["x"]){
   12             arr.push({x:spek2[j]["x"], y: (spek1[i]["y"] + (spek1[i+1]["y"]- spek1[i]["y"]) / (spek1[i+1]["x"] - spek1[i]["x"]) * (spek2[j]["x"] - spek1[i]["x"])) * spek2[j]["y"] });
   13         };
   14 
   15         if(spek1[i]["x"]>spek2[j]["x"]){
   16             arr.push({x:spek1[i]["x"], y: (spek2[j]["y"] + (spek2[j+1]["y"]- spek2[j]["y"]) / (spek2[j+1]["x"] - spek2[j]["x"]) * (spek1[i]["x"] - spek2[j]["x"])) * spek1[i]["y"] });
   17         };
   18         if(spek1[i+1]["x"]==spek2[j+1]["x"]){i++; j++; continue;};
   19         if(spek1[i+1]["x"]<spek2[j+1]["x"]){i++; continue;};
   20         if(spek1[i+1]["x"]>spek2[j+1]["x"]){j++; continue;};
   21     }
   22     /*************Letztes Element************/
   23     if(i==spek1.length-1){
   24         arr.push({x:spek1[i]["x"], y: (spek2[j-1]["y"] + (spek2[j]["y"]- spek2[j-1]["y"]) / (spek2[j]["x"] - spek2[j-1]["x"]) * (spek1[i]["x"] - spek2[j-1]["x"])) * spek1[i]["y"] });
   25         //console.log("i:arr[arr.length-1]",arr[arr.length-1]);
   26     }else{
   27         arr.push({x:spek2[j]["x"], y: (spek1[i-1]["y"] + (spek1[i]["y"]- spek1[i-1]["y"]) / (spek1[i]["x"] - spek1[i-1]["x"]) * (spek2[j]["x"] - spek1[i-1]["x"])) * spek2[j]["y"] });
   28     }
   29     return arr;
   30 }
In die Programmiersprache VBA umgesetzt wäre die in Listing [*] umgesetzte Funktion eine Möglichkeit. Übergeben werden zwei Arrays U und V und der Rückgabewert ist das durchmultiplizierte Array.
Listing: Funkton Falten in VBA 
    1 
    2 Public Function Durchmultiplizieren(ByRef U As Variant, ByRef V As Variant) As Variant()
    3 Dim F() As Variant: ReDim F(1, 0)
    4 Dim i As Long, j As Long: i = 0: j = 0
    5 Do While i < UBound(V, 2) And j < UBound(U, 2)
    6     If V(0, i) >= U(0, j + 1) Then j = j + 1: GoTo ContinueLoop
    7     If U(0, j) >= V(0, i + 1) Then i = i + 1: GoTo ContinueLoop
    8     If U(0, j) = V(0, i) Then
    9         If Not IsEmpty(F(1, UBound(F, 2))) Then ReDim Preserve F(1, UBound(F, 2) + 1)
   10         F(0, UBound(F, 2)) = V(0, i)
   11         F(1, UBound(F, 2)) = U(1, j) * V(1, i)
   12         'Debug.Print i; j, F(0, UBound(F, 2)); F(1, UBound(F, 2)); "=  "; U(1, j); "*"; V(1, i)
   13         GoTo Weiter
   14     End If
   15     If U(0, j) < V(0, i) Then
   16         If Not IsEmpty(F(1, UBound(F, 2))) Then ReDim Preserve F(1, UBound(F, 2) + 1)
   17         F(0, UBound(F, 2)) = V(0, i)
   18         F(1, UBound(F, 2)) = (U(1, j) + (U(1, j + 1) - U(1, j)) / (U(0, j + 1) - U(0, j)) * (V(0, i) - U(0, j))) * V(1, i)
   19         'Debug.Print i; j, F(0, UBound(F, 2)); F(1, UBound(F, 2)); "= ("; U(1, j); "+("; U(1, j + 1); "-"; U(1, j); ")/("; U(0, j + 1); "-"; U(0, j); ")*("; V(0, i); "-"; U(0, j); "))*"; V(1, i)
   20         GoTo Weiter
   21     End If
   22     If U(0, j) > V(0, i) Then
   23         If Not IsEmpty(F(1, UBound(F, 2))) Then ReDim Preserve F(1, UBound(F, 2) + 1)
   24         F(0, UBound(F, 2)) = U(0, j)
   25         F(1, UBound(F, 2)) = (V(1, i) + (V(1, i + 1) - V(1, i)) / (V(0, i + 1) - V(0, i)) * (U(0, j) - V(0, i))) * U(1, j)
   26         'Debug.Print i; j, F(0, UBound(F, 2)); F(1, UBound(F, 2)); "= ("; V(1, i); "+("; V(1, i + 1); "-"; V(1, i); ")/("; V(0, i + 1); "-"; V(0, i); ")*("; U(0, j); "-"; V(0, i); "))*"; U(1, j)
   27         GoTo Weiter
   28     End If
   29 Weiter:
   30     If V(0, i + 1) = U(0, j + 1) Then i = i + 1: j = j + 1: GoTo ContinueLoop
   31     If V(0, i + 1) < U(0, j + 1) Then i = i + 1: GoTo ContinueLoop
   32     If V(0, i + 1) > U(0, j + 1) Then j = j + 1: GoTo ContinueLoop
   33 ContinueLoop:
   34 Loop
   35 '***************letztes Element******************
   36 If j = UBound(U, 2) Then
   37     i = IIf(i = 0, 1, i)
   38     'Debug.Print ; ; , V(0, i); (U(1, j - 1) + (U(1, j) - U(1, j - 1)) / (U(0, j) - U(0, j - 1)) * (V(0, i) - U(0, j - 1))) * V(1, i);
   39     'Debug.Print "= ("; U(1, j - 1); "+("; U(1, j); "-"; U(1, j - 1); ")/("; U(0, j); "-"; U(0, j - 1); ")*("; V(0, i); "-"; U(0, j - 1); "))*"; V(1, i)
   40     If Not IsEmpty(F(1, UBound(F, 2))) Then ReDim Preserve F(1, UBound(F, 2) + 1)
   41     F(0, UBound(F, 2)) = U(0, j)
   42     F(1, UBound(F, 2)) = (V(1, i - 1) + (V(1, i) - V(1, i - 1)) / (V(0, i) - V(0, i - 1)) * (U(0, j) - V(0, i - 1))) * U(1, j)
   43     'Debug.Print "X"
   44 ElseIf i = UBound(V, 2) Then
   45     j = IIf(j = 0, 1, j)
   46     'Debug.Print ; ; , V(0, i); (U(1, j - 1) + (U(1, j) - U(1, j - 1)) / (U(0, j) - U(0, j - 1)) * (V(0, i) - U(0, j - 1))) * V(1, i);
   47     'Debug.Print "= ("; U(1, j - 1); "+("; U(1, j); "-"; U(1, j - 1); ")/("; U(0, j); "-"; U(0, j - 1); ")*("; V(0, i); "-"; U(0, j - 1); "))*"; V(1, i)
   48     If Not IsEmpty(F(1, UBound(F, 2))) Then ReDim Preserve F(1, UBound(F, 2) + 1)
   49     F(0, UBound(F, 2)) = V(0, i)
   50     F(1, UBound(F, 2)) = (U(1, j - 1) + (U(1, j) - U(1, j - 1)) / (U(0, j) - U(0, j - 1)) * (V(0, i) - U(0, j - 1))) * V(1, i)
   51 End If
   52 Durchmultiplizieren = F
   53 End Function

Subsections
next up previous index
Next: Praktische Umsetzung als Excel-Funktion Up: RechnenmitSpektren Previous: Addition zweier Spektren   Index