diff options
Diffstat (limited to 'buch/papers')
| -rw-r--r-- | buch/papers/erdbeben/Apperatur.jpg | bin | 0 -> 66346 bytes | |||
| -rw-r--r-- | buch/papers/erdbeben/Gausskurve2.pdf | bin | 0 -> 26978 bytes | |||
| -rw-r--r-- | buch/papers/erdbeben/Gausskurve2.tex | 39 | ||||
| -rw-r--r-- | buch/papers/erdbeben/Gausskurve3.pdf | bin | 0 -> 27445 bytes | |||
| -rw-r--r-- | buch/papers/erdbeben/Gausskurve3.tex | 47 | ||||
| -rw-r--r-- | buch/papers/erdbeben/teil1.tex | 289 | ||||
| -rw-r--r-- | buch/papers/mceliece/example_code/mceliece_simple.py | 327 |
7 files changed, 662 insertions, 40 deletions
diff --git a/buch/papers/erdbeben/Apperatur.jpg b/buch/papers/erdbeben/Apperatur.jpg Binary files differnew file mode 100644 index 0000000..d25381e --- /dev/null +++ b/buch/papers/erdbeben/Apperatur.jpg diff --git a/buch/papers/erdbeben/Gausskurve2.pdf b/buch/papers/erdbeben/Gausskurve2.pdf Binary files differnew file mode 100644 index 0000000..bee3bc0 --- /dev/null +++ b/buch/papers/erdbeben/Gausskurve2.pdf diff --git a/buch/papers/erdbeben/Gausskurve2.tex b/buch/papers/erdbeben/Gausskurve2.tex new file mode 100644 index 0000000..44319c3 --- /dev/null +++ b/buch/papers/erdbeben/Gausskurve2.tex @@ -0,0 +1,39 @@ +\documentclass{standalone} + +\usepackage{pgfplots} + +\pgfplotsset{compat = newest} + +\begin{document} + + +\begin{tikzpicture} + + +\begin{axis}[ + xmin = -1, xmax = 4, + ymin = -0.5, ymax = 2.50, + axis lines = center, + xlabel = $\sigma$, + ylabel = {$\mu$}, +] + +\addplot [ + domain=-2:5, + samples=200, + color=orange, +] +{(1/(2*pi*0.2^2)^0.2)*exp(-(x-2)^2/(2*0.2^2))}; + +\addplot [ + domain=-2:5, + samples=200, + color=blue, + ] + {1/(2*pi*0.5^2)^0.5)*exp(-(x-0.9)^2/(2*0.5^2))}; + +\end{axis} +\end{tikzpicture} + + +\end{document}
\ No newline at end of file diff --git a/buch/papers/erdbeben/Gausskurve3.pdf b/buch/papers/erdbeben/Gausskurve3.pdf Binary files differnew file mode 100644 index 0000000..e86a403 --- /dev/null +++ b/buch/papers/erdbeben/Gausskurve3.pdf diff --git a/buch/papers/erdbeben/Gausskurve3.tex b/buch/papers/erdbeben/Gausskurve3.tex new file mode 100644 index 0000000..85455ef --- /dev/null +++ b/buch/papers/erdbeben/Gausskurve3.tex @@ -0,0 +1,47 @@ +\documentclass{standalone} + +\usepackage{pgfplots} + +\pgfplotsset{compat = newest} + +\begin{document} + + +\begin{tikzpicture} + + +\begin{axis}[ + xmin = -1, xmax = 4, + ymin = -0.5, ymax = 2.50, + axis lines = center, + xlabel = $\sigma$, + ylabel = {$\mu$}, +] + +\addplot [ + domain=-2:5, + samples=200, + color=orange, +] +{(1/(2*pi*0.2^2)^0.2)*exp(-(x-2)^2/(2*0.2^2))}; + +\addplot [ + domain=-2:5, + samples=200, + color=blue, + ] + {1/(2*pi*0.5^2)^0.5)*exp(-(x-0.9)^2/(2*0.5^2))}; + +\addplot [ + domain=-2:5, + samples=200, + color=red, + ] + {((1/(2*pi*0.5^2)^0.5)*exp(-(x-0.9)^2/(2*0.5^2))*(1/(2*pi*0.2^2)^0.2)*exp(-(x-2)^2/(2*0.2^2)))/0.1}; + + +\end{axis} +\end{tikzpicture} + + +\end{document}
\ No newline at end of file diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index a89f303..0d21f84 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -3,16 +3,240 @@ % % (c) 2020 Prof Dr Andreas Müller, Hochschule Rapperswil % -\section{Teil 1 -\label{erdbeben:section:teil1}} -\rhead{Problemstellung} -Sed ut perspiciatis unde omnis iste natus error sit voluptatem -accusantium doloremque laudantium, totam rem aperiam, eaque ipsa -quae ab illo inventore veritatis et quasi architecto beatae vitae -dicta sunt explicabo. -Nemo enim ipsam voluptatem quia voluptas sit aspernatur aut odit -aut fugit, sed quia consequuntur magni dolores eos qui ratione -voluptatem sequi nesciunt +% +% teil2.tex -- Beispiel-File für teil2 +% +% (c) 2020 Prof Dr Andreas Müller, Hochschule Rapperswil +% + + + +\section{Kalman Filter} +\subsection{Geschichte} +Das Kalman Filter wurde 1960 von Rudolf Emil Kalman entdeckt und direkt von der NASA für die Appollo Mission benutzt. Der Filter kommt mit wenig Rechenleistung aus und war somit dafür geeignet die Rakete bei der Navigation zu unterstützen. Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den nächsten Zustand errechnen. Typische Anwendungen des Kalman-Filters sind die Glättung von verrauschten Daten und die Schätzung von Parametern und kommt heutzutage in jedem Satellit, Navigationssystem, Smartphones und Videospielen vor. + +\subsection{Wahrscheinlichkeit} +Das Kalman Filter versucht nichts anderes, als ein geeigneter Wert zwischen zwei Normalverteilungen zu schätzen. Die eine Kurve zeigt die errechnete Vorhersage des Zustands, bzw. deren Normal- Gauss-Verteilung. Die andere Kurve zeigt die verrauschte Messung des nächsten Zustand, bzw. deren Normal-Verteilung. Wie man in am Beispiel dieser zwei Gauss-Verteilungen sehen kann, ist sowohl der geschätzte Zustand als auch der gemessene Zustand nicht am selben Punkt. + + + +\begin{figure} + \begin{center} + \includegraphics[width=5cm]{papers/erdbeben/Gausskurve2.pdf} + \caption{System} + \end{center} +\end{figure} + + + +Um eine genauere Schätzung des Zustandes zu machen, wird nun ein Wert zwischen den beiden Verteilungen gesucht. An diesem Punkt wird nun eine Eigenschaft ausgenutzt. Durch das Multiplizieren zweier Normalverteilungen entsteht eine neue Normalverteilung. + +Wir haben eine Normalverteilung der Vorhersage: +\begin{equation} +{y_1}(x;{\mu_1},{\sigma_1})=\frac{1}{\sqrt{2\pi\sigma_1^2}}\quad e^{-\frac{(x-{\mu_1})^2}{2{\sigma_1}^2}} +\end{equation} +und für die Messung: + +\begin{equation} +{y_2}(x;{\mu_2},{\sigma_2})=\frac{1}{\sqrt{2\pi\sigma_2^2}}\quad e^{-\frac{(x-{\mu_2})^2}{2{\sigma_2}^2}}. +\end{equation} + +Diesen werden nun Multipliziert und durch deren Fläche geteilt um sie wieder zu Normieren: +\begin{equation} +{y_f}(x;{\mu_f},{\sigma_f})=\frac{ \frac{1}{\sqrt{2\pi\sigma_1^2}}e^{-\frac{(x-{\mu_1})^2}{2{\sigma_1}^2}} \cdot \frac{1}{\sqrt{2\pi\sigma_2^2}}e^{-\frac{(x-{\mu_2})^2}{2{\sigma_2}^2}}}{\int {y_1}*{y_2}\,} +\end{equation} + +Dadurch gleicht sich die neue Kurve den anderen an. Interessant daran ist, dass die fusionierte Kurve sich der genauere Normal-Verteilung anpasst. ist ${\sigma_2}$ klein und ${\sigma_1}$ gross, so wird sich die fusionierte Kurve näher an ${y_2}(x;{\mu_2},{\sigma_2})$ begeben. Sie ist also Gewichtet und die best mögliche Schätzung. + + +\begin{figure} + \begin{center} + \includegraphics[width=5cm]{papers/erdbeben/Gausskurve3.pdf} + \caption{System} + \end{center} +\end{figure} + + +Was in 2 Dimensionen erklärt wurde, funktioniert auch in mehreren Dimensionen. Dieses Prinzip mach sich der Kalman Filter zu nutze, und wird von uns für die Erdbeben Berechnung genutzt. + +\subsection{Anwendungsgrenzen} +Nicht lineare Systeme %Noch nicht Fertig + + +\section{Aufbau} +Um ein Erdbeben kenntlich zumachen werden in der Regel Seismographen mit vielen Sensoren verwendet. +Ein Seismograph besteht im Grunde aus einer federgelagerten Masse. Wirkt eine Bodenerregung auf das Gerät ein, bleibt die gekoppelte Masse in der regel stehen und das Gehäuse schwingt mit.Relativbewegung des Bodens kann damit als Längenänderung im Zeitverlauf gemessen werden. In modernen Seismographen wird die Bodenbewegung in alle Richtungen gemessen, sowohl Horizontal als auch Vertikal. +Wir konstruieren uns eine einfachere Version eines Seismographen, welcher rein mechanisch funktioniert. Zudem kann er nur in eine Dimension Messwerte aufnehmen. Würde das System ausgebaut werden, um alle Horizontalbewegungen aufzunehmen, würde der Verwendung des Kalman-Filters zu kompliziert werden. Für zwei Dimensionen (x,y) würde der Pythagoras für das System benötigt werden. Da sich der Pythagoras bekanntlich nicht linear verhält, kann kein lineares Kalman-Filter implementiert werden. Da das Kalman-Filter besonders effektiv und einfach für lineare Abläufe geeignet ist, würde eine Zweidimensionale Betrachtung den Rahmen dieser Arbeit sprengen. Für ein nicht-lineares System werden Extended Kalman-Filter benötigt, bei denen die System-Matrix (A) durch die Jacobi-Matrix des System ersetzt wird. + +\begin{figure} + \begin{center} + \includegraphics[width=5cm]{papers/erdbeben/Apperatur} + \caption{System} + \end{center} +\end{figure} + + +\subsection{Optionen} +Wollte man einen 2D Seismographen aufbauen, ohne den Pythagroas zu verwenden, kann dies mit der Annahme, das die Feder sehr lang sind erfolgen. Da sich bei langen Federn die Auslenkungen verkleiner...!!Noch nicht fertig! + +\section{Systemgleichung} +Da das Kalman-Filter zum Schätzen des nächsten Zustand verwendet wird, wird eine Gleichung, welche das System beschreibt. Das Kalman-Filter benötigt eine Beschreibung der Systemdynamik. Im Fall unseres Seismographen, kann die Differentialgleichung zweiter Ordnung einer gedämpften Schwingung am harmonischen Oszillator verwendet werden. Diese lautet: +\begin{equation} +m\ddot x + 2k \dot x + Dx = f +\end{equation} +mit den Konstanten $m$ = Masse, $k$ = Dämpfungskonstante und $D$ = Federkonstante. +Um diese nun in die Systemmatrix umzuwandeln, wird aus der Differentialgleichung zweiter Ordnung durch eine Substitution eine DGL erster Ordnung: + + +\begin{equation} +{x_1}=x, \qquad +{x_2}=\dot x, \qquad +{x_3}=\ddot x\qquad \mid \quad \text {Substitution} +\end{equation} + + +\begin{equation} +m{x_3}+ 2k{x_2} + D{x_1} = f\qquad \mid \quad \text {DGL 1. Ordnung} +\end{equation} + +\begin{equation} +{x_3}=-\frac{D}{m} {x_1} -\frac{2k}{m} {x_2} + \frac{f} {m} \qquad \mid \quad \text {nach} \quad{x_3} +\end{equation} +auch als Matrix-Vektor-Gleichung schreiben. +Hierbei beschreibt die Matrix $A$ die gesamte Systemdynamik in der Form, wie sie ein Kalman-Filter benötigt. + +Um die lineare Differentialgleichung in das Kalman-Filter zu implementieren, muss dieses als Vektor-Gleichung umgewandelt werden. Dafür wird die Gleichung in die Zustände aufgeteilt. Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und äussere Beschleunigung des ganzen System. Dabei muss unterschieden werden. um welche Beschleunigung es sich handelt. Das System beinhaltet sowohl eine Beschleunigung der Masse bzw. Feder (innere Beschleunigung), als auch eine Beschleunigung der ganzen Apparatur (äussere Beschleunigung). In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbeben Anregung gleich kommt. + + +\begin{equation} +\frac{d}{dt} \left(\begin{array}{c} {x_1} \\ {x_2} \end{array}\right) = \left( + \begin{array}{ccc} +0 & 1& 0 \\ +- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ +\end{array}\right) \left(\begin{array}{c} {x_1} \\ {x_2} \\ {x_3} \end{array}\right). +\end{equation} + +Durch die Rücksubstituion ergibt sich: +\begin{equation} +\frac{d}{dt} \left(\begin{array}{c} x(t) \\ v(t) \end{array}\right) = \left( + \begin{array}{ccc} +0 & 1& 0 \\ +- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ +\end{array}\right) \left(\begin{array}{c} x(t)\\ v(t)\\ f(t) \end{array}\right). +\end{equation} + + +Da die Kraft unbekannt ist, wird die letzte Zeile später mit Nullen bestückt, denn genau diese Werte wollen wir. + +\section{Kalman Filter} +Um den Kalman Filter zu starten, müssen gewisse Bedingungen definiert werden. In diesem Abschnitt werden die einzelnen Parameter/Matrizen erläutert und Erklärt, wofür sie nützlich sind. + + +\subsection{Anfangsbedingungen} +\subsubsection*{Anfangszustand $x$} +Das Filter benötigt eine Anfangsbedingung. In unserem Fall ist es die Ruhelage, die Masse bewegt sich nicht. Zudem erföhrt die Apparatur keine äussere Kraft. + +\begin{equation} +{x_0 }= \left( \begin{array}{c} 0\\ 0\\ 0\end{array}\right) +\end{equation} + +\subsubsection*{Anfangsfehler / Kovarianzmatrix $P$} +Da auch der Anfangszustand fehlerhaft sein kann, wird für den Filter einen Anfangsfehler eingeführt. Auf der Diagonalen werden die Varianzen eingesetzt, in den restlichen Felder stehen die Kovarianzen. +In unserem Fall ist der Anfangszustand gut bekannt. Wir gehen davon aus, dass das System in Ruhe und in Abwesenheit eines Erdbeben startet, somit kann die Matrix mit Nullen bestückt werden. Somit ergibt sich für die Kovarianzmatrix + +\begin{equation} +{P_0 }= +\left( +\begin{array}{ccc} +0 & 0 &0 \\ +0 &0 & 0 \\ +0 & 0 &0 \\ +\end{array} +\right). +\end{equation} +Diese Matrix beschreibt die Unsicherheit des geschätzten Zustandes und wird sowohl für die Vorhersage als auch die Korrektur benötigt. Sie wird nach jeder Schätzung aktualisiert.. Für einen gut bekannten Zustandsvektor können kleine Werte eingesetzt werden, für ungenaue Anfangsbedingungen sollten grosse Werte (1 Million) verwendet werden. Grosse Werte ermöglichen dem Filter sich schnell einzupendeln. + + +\subsubsection*{Dynamikmatrix $A$} +Die Dynamikmatrix bildet den Kern des Filters. Diese wurde weiter oben Bereits beschrieben. Dabei wollen wird die äussere Kraft des Systems ermitteln. +Da nichts über die äussere Kraft bekannt ist, müssen wir annehmen das deren Ableitung 0 ist. +Die System Vektor-Gleichung lautet daher: + + +\begin{equation} +A = \left( + \begin{array}{ccc} +0 & 1& 0 \\ +- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ +0 & 0& 0\\ +\end{array}\right) +\end{equation} + +\subsubsection*{Prozessrauschkovarianzmatrix $Q$} +Die Prozessrauschmatrix teilt dem Filter mit, wie sich der Systemzustand verändert. Kalman-Filter berücksichtigen Unsicherheiten wie Messfehler und -rauschen. Bei unserem Modell könnte das beispielsweise ein Windstoss an die Masse sein. Für uns wäre dies: +\begin{equation} +Q = \left( + \begin{array}{ccc} +{\sigma_x }^2& 0& 0 \\ +0 & {\sigma_v }^2& 0\\ +0 & 0& {\sigma_f }^2\\ +\end{array}\right) +\end{equation} + +Die Standabweichungen müssten Statistisch ermittelt werden, da der Fehler nicht vom Sensor kommt und somit nicht vom Hersteller gegeben ist. Das Bedeutet wiederum dass $Q$ die Unsicherheit des Prozesses beschreibt, und die Messung. + +\subsubsection*{Messmatrix $H$} +Die Messmatrix gibt an, welcher Parameter gemessen werden soll. In unserem Fall ist es nur die Position der Massen. + +\[ H = (1, 0, 0) \] + + +\subsubsection*{Messrauschkovarianz $R$} +Die Messrauschkovarianzmatrix beinhaltet, wie der Name es schon sagt, das Rauschen der Positionssensoren. In unserem Fall wird nur die Position der Masse gemessen. Da wir keine anderen Sensoren haben ist dies lediglich: +\begin{equation} +R= ({\sigma_x}^2). +\end{equation} +Diese Messrauchen wird meistens vom Sensorhersteller angegeben. Für unsere Theoretische Apparatur wird hier ein kleiner Fehler eingesetzt. + +\subsection{Fiter Algorithmus} +Nachdem alle Parameter aufgestellt sind, wird der Filter initialisiert und wird den Zustand der Feder vorherzusagen, die Messung zu präzisieren und laufend zu aktualisieren. Das Filter berechnet aufgrund der aktuellen Schätzung eine Vorhersage. Diese wird, sobald verfügbar, mit der Messung verglichen. Aus dieser Differenz und den Unsicherheiten des Prozesses ($Q$) und der Messung ($R$) wird der wahrscheinlichste, neue Zustand geschätzt. + + +\subsubsection*{Vorhersage} +Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. Dies funktioniert ganz Trivial mit dem Rechenschritt: +\begin{equation} +{x_{t+1}}=A\cdot{x_t}. +\end{equation} + + +Die Kovarianz $P_{pred}$ wird ebenfalls neu berechnet, da die Unsicherheit im Vorhersage grösser wird als im Aktuellen. Da wir ein mehrdimensionales System haben, kommt noch die Messunsicherheit $Q$ dazu, so dass die Unsicherheit des Anfangsfehlers $P$ immer grösser wird. Dies funktioniert durch multiplizieren der Systemmatrix, deren Ableitung und mit dem aktualisierten Anfangsfehler. Dazu wird noch die Messunsicherheit addiert, somit entsteht die Gleichung + +\begin{equation} +{P_{pred}}=APA^T+Q. +\end{equation} + +wird dieser Vorgang wiederholt, schaut der Filter wie genau die letzte Anpassung von $P$ zur Messung stimmt. Ist der Unterschied klein, wird die Kovarianz $P$ kleiner, ist der Unterschied gross, wird auch die Kovarianz grösser. Das Filter passt sich selber an und korrigiert sich bei grosser Abweichung. + +\subsubsection*{Messen} +Der Sensor wurde noch nicht benutz, doch genau der liefert Werte für den Filter. Die aktuellen Messwerte $z$ werden die Innovation $w$ mit dem Zustandsvektor $x$ und der Messmatrix $H$ zusammengerechnet. +Hier bei wird lediglich die Messung mit dem Fehler behaftet, und die Messmatrix $H$ +\begin{equation} +w=Z-(H\cdot x) +\end{equation} +Die Innovation ist der Teil der Messung, die nicht durch die Systemdynamik erklärt werden kann. +Innovation = Messung - Vorhersage. Dies ist Intuitiv logisch, eine Innovation von 0 bedeutet, dass die Messung nichts Neues hervorbrachte. + +Im nächsten Schritt wir analysiert, mit welcher Kovarianz weiter gerechnet wird. + +\subsubsection*{Korrigieren} +Udpdate +\section{Anfügen der Schwingung} + +Ein Erdbeben breitet sich im Boden wellenartig aus und bringt Objekte, wie zum Beispiel ein Gebäude, in Schwingung. +Diese Schwingungen pflanzen sich im Gebäude mit gleicher Amplitude, Geschwindigkeit und Beschleunigung in horizontaler und vertikaler Bewegung fort. +Wir möchten herauszufinden, wie gross die Massenbeschleunigung infolge eines Erdbeben ist. +Mit Hilfe von fiktiven Sensoren, die eine Ortsveränderung des Gebäude messen, können wir mit Anwendung von Matrizen und dem Kalman-Filter die Beschleunigung berechnen. + \begin{equation} \int_a^b x^2\, dx = @@ -21,35 +245,20 @@ voluptatem sequi nesciunt \frac{b^3-a^3}3. \label{erdbeben:equation1} \end{equation} -Neque porro quisquam est, qui dolorem ipsum quia dolor sit amet, -consectetur, adipisci velit, sed quia non numquam eius modi tempora -incidunt ut labore et dolore magnam aliquam quaerat voluptatem. - -Ut enim ad minima veniam, quis nostrum exercitationem ullam corporis -suscipit laboriosam, nisi ut aliquid ex ea commodi consequatur? -Quis autem vel eum iure reprehenderit qui in ea voluptate velit -esse quam nihil molestiae consequatur, vel illum qui dolorem eum -fugiat quo voluptas nulla pariatur? - -\subsection{De finibus bonorum et malorum -\label{erdbeben:subsection:finibus}} -At vero eos et accusamus et iusto odio dignissimos ducimus qui -blanditiis praesentium voluptatum deleniti atque corrupti quos -dolores et quas molestias excepturi sint occaecati cupiditate non -provident, similique sunt in culpa qui officia deserunt mollitia -animi, id est laborum et dolorum fuga \eqref{000tempmlate:equation1}. - -Et harum quidem rerum facilis est et expedita distinctio -\ref{erdbeben:section:loesung}. -Nam libero tempore, cum soluta nobis est eligendi optio cumque nihil -impedit quo minus id quod maxime placeat facere possimus, omnis -voluptas assumenda est, omnis dolor repellendus -\ref{erdbeben:section:folgerung}. -Temporibus autem quibusdam et aut officiis debitis aut rerum -necessitatibus saepe eveniet ut et voluptates repudiandae sint et -molestiae non recusandae. -Itaque earum rerum hic tenetur a sapiente delectus, ut aut reiciendis -voluptatibus maiores alias consequatur aut perferendis doloribus -asperiores repellat. + +\section{Erreger-Schwingung} +Wir möchten mit einer gedämpften harmonischen Schwingung ein einfaches Erdbeben simulieren, die im Kalman Filter eingespeist wird. +Die Gleichung lautet + +\begin{equation} +x(t)=Ae^{t/2}sin(t). +\end{equation} + +Mit dieser Schwingung können wir ein einachsiger Seismograph simulieren, der eine Ortsverschiebung auf der x-Achse durchführt. +Die Dämpfung der Schwingung ist relevant, da das System beim Schwingungsvorgang durch die Federkonstante und der Reibung, Energie verliert. + +Die Ergebnisse dieser Schwingung setzen wir in die Messmatrix ein und können den Kalman-Filter starten. + + diff --git a/buch/papers/mceliece/example_code/mceliece_simple.py b/buch/papers/mceliece/example_code/mceliece_simple.py new file mode 100644 index 0000000..bac3b42 --- /dev/null +++ b/buch/papers/mceliece/example_code/mceliece_simple.py @@ -0,0 +1,327 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Apr 18 08:58:30 2021 + +@author: terminator +""" +import numpy as np +from numpy.polynomial import Polynomial as Poly + + +def gen_perm_M(n): + ''' + generating random binary permutation matrix: https://en.wikipedia.org/wiki/Permutation_matrix + + Parameters + ---------- + n : int + size of output matrix. + + Returns + ------- + perm_M : np.ndarray + binary permutation matrix. + + ''' + perm_M=np.zeros((n, n), dtype=int) + perm_V=np.random.default_rng().permutation(n) + perm_M[perm_V, np.r_[0:n]]=1 + return perm_M + + + +def create_linear_code_matrix(n, k, g): + ''' + create generator matrix for linear encoding + use this matrix to create code_vector: matrix @ data_vector=code_vector + + Parameters + ---------- + n : int + len of code_vector (matrix @ data_vector). + k : int + len of data_vector. + g : numpy.Polynominal + polynomina wich defines the the constellations of . + + Returns + ------- + np.ndarray + Generator matrix + + ''' + if n != k+len(g)-1: + raise Exception("n, k not compatible with degree of g") + rows=[] + for i in range(k): #create potences of p + row=np.r_[np.zeros(i), g.coef, np.zeros(n-k-i)] + rows.append(row) + return np.array(rows, dtype=int).T + + +def gf2_inv(M, get_full=False): + ''' + create inverse of matrix M in gf2 + + Parameters + ---------- + M : numpy.ndarray + input Matrix. + get_full : Bool, optional + if Ture, return inverse of G with I on the left (if gaussian inversion successful) If True, vaidity is not proven. The default is False. + + Returns + ------- + np.ndarray or None + returns inverse if M was not singular in gf2, else None. + + ''' + size=len(M) + G=np.hstack((M, np.eye(size))) + G=np.array(G, dtype=int) + for n in range(size): #forward reduction + if G[n, n] == 0: #swap line if necessary + for i in range(n+1, size): + if G[i, n]: + G[i, :], G[n, :] = G[n, :].copy(), G[i, :].copy() #swap + + for i in range(n+1, size): #downward reduction + if G[i, n]: + G[i, :] = G[i, :] ^ G[n, :] + #reached buttom_right with pivo + for n in range(size-1, -1, -1): #backwards + for i in range(n): + if G[i, n]: + G[i, :]= G[i, :] ^ G[n, :] + + + + if get_full: + return G + else: + valid=np.sum(np.abs(G[:, :size]-np.eye(size))) == 0 #reduction successfull when eye left on the left of G + if not valid: + return None + else: + return G[:, size:] + +def create_rand_bin_M(n, with_inverse=False): + ''' + create random binary matrix + + Parameters + ---------- + n : int + size. + with_inverse : bool, optional + if False, return only random binary matrix. + if True, return also inverse of random martix. for this purpose, random matrix will not be singular. + The default is False. + + Returns + ------- + M : TYPE + if with_inverse is True: return tuple(random_matrix, inverse_of_random_matrix) + else: return random_matrix + + ''' + inv=None + while type(inv) == type(None): #do it until inversion of m is successful wich means det(M)%2 != 0 + M=np.random.randint(0,2, (n,n)) + inv=gf2_inv(M) + if with_inverse: + return(M, inv) + else: + return M + +def create_syndrome_table(n, g): + ''' + create syndrome table for correcting errors in linear code + + Parameters + ---------- + n : int + len of linear-code code_vector. + g : numpy.Polinominal + generator polynominal used by linear-code-encoder. + + Returns + ------- + list of error_vectors, one per syndrome. + get the corresponding error_vector by using the value represented by your syndrome as index of this list. + + ''' + zeros=np.zeros(n, dtype=int) + syndrome_table=[0 for i in range(n+1)] + syndrome_table[0]=zeros #when syndrome = 0, no bit-error to correct + for i in range(n): + faulty_vector=zeros.copy() + faulty_vector[i]=1 + q, r=divmod(Poly(faulty_vector), g) + r=np.r_[r.coef%2, np.zeros(len(g)-len(r))] + index=np.sum([int(a*2**i) for i, a in enumerate(r)]) + syndrome_table[index]=faulty_vector + return np.array(syndrome_table) + + +def decode_linear_code(c, g, syndrome_table): + ''' + function to decode codeword encoded with linear_code + + Parameters + ---------- + c : list or np.ndarray + code_vector. + g : numpy.Polynominal + generator polynominal, used to create code_vector. + syndrome_table : list of error_vectors + if codeword contains an error, syndrome_table is used to correct wrong bit. + + Returns + ------- + numpy.ndarray + data_vector. + + ''' + q, r=divmod(Poly(c), g) + q=np.r_[q.coef%2, np.zeros(len(c)-len(q)-len(g)+1)] + r=np.r_[r.coef%2, np.zeros(len(g)-len(r))] + syndrome_index=np.sum([int(a*2**i) for i, a in enumerate(r)]) + while syndrome_index > 0: + c=c ^ syndrome_table[syndrome_index] + q, r=divmod(Poly(c), g) + q=np.r_[q.coef%2, np.zeros(len(c)-len(q)-len(g)+1)] + r=np.r_[r.coef%2, np.zeros(len(g)-len(r))] + syndrome_index=np.sum([int(a*2**i) for i, a in enumerate(r)]) + return np.array(q, dtype=int) + +def encode_linear_code(d, G): + ''' + uses generator matrix G to encode d + + Parameters + ---------- + d : numpy.ndarray + data_vector. + G : numpy.ndarray + generator matrix. + + Returns + ------- + c : numpy.ndarray + G @ d. + + ''' + c=(G @ d)%2 + return c + +def encrypt(d, pub_key, t): + ''' + encrypt data with public key, adding t bit-errors + + Parameters + ---------- + d : numpy.ndarray or list of numpy.ndarray + data_vector or list of data vectors to encrypt. + pub_key : numpy.ndarray + public key matrix used to encrypt data. + t : int + number of random errors to add to code_vector. + + Returns + ------- + numpy.ndarray or list of numpy.ndarray (depending on d) + encrypted data. + + ''' + if type(d) in (list,): + encrypted_list=[encrypt(data, pub_key, t) for data in d] + return encrypted_list + else: + c = np.array((pub_key @ d)%2, dtype=int) + + indexes_of_errors=np.random.default_rng().permutation(pub_key.shape[0])[:t] #add t random errors to codeword + e=np.zeros(pub_key.shape[0], dtype=int) + e[indexes_of_errors]=1 + c=c ^ e + return np.array(c, dtype=int) + +def decrypt(c, P_inv, linear_code_decoder, S_inv): + ''' + decrypt encrypted message + + Parameters + ---------- + c : numpy.ndarray or list of numpy.ndarray + code_vector or list of code_vectors to decrypt. + P_inv : numpy.ndarray + inverted permutation matrix. + linear_code_decoder : function(x) + function to use to decode linear code. + S_inv : numpy.ndarray + inverted random binary matrix. + + Returns + ------- + numpy.ndarray or list of numpy.ndarray (depending on d) + decrypted data. + + ''' + if type(c) in (list,): + decrypted_list=[decrypt(codew, P_inv, linear_code_decoder, S_inv) for codew in c] + return decrypted_list + else: + c=np.array(P_inv @ c, dtype=int)%2 + d=linear_code_decoder(c) + d=(S_inv @ d)%2 + return np.array(d, dtype=int) + +def str_to_blocks4(string): + blocks=[] + for char in string: + bits=[int(value) for value in np.binary_repr(ord(char), 8)[::-1]] #lsb @ index 0 + blocks.append(np.array(bits[:4], dtype=int)) #lower nibble first + blocks.append(np.array(bits[4:], dtype=int)) + return blocks + +def blocks4_to_str(blocks): + string='' + for i in range(0, len(blocks), 2): + char=np.sum([b*2**i for i, b in enumerate(blocks[i])]) + \ + np.sum([b*2**(i+4) for i, b in enumerate(blocks[i+1])]) + string+=chr(char) + return string + +if __name__ == '__main__': + + #shared attributes: + n=7 + k=4 + t=1 + + #private key(s): + g=Poly([1,1,0,1]) #generator polynom for 7/4 linear code (from table, 1.0 + 1.0·x¹ + 0.0·x² + 1.0·x³) + P_M=gen_perm_M(n) #create permutation matrix + G_M=create_linear_code_matrix(n, k, g) #linear code generator matrix + S_M, S_inv=create_rand_bin_M(k, True) #random binary matrix and its inverse + P_M_inv=P_M.T #inverse permutation matrix + + syndrome_table=create_syndrome_table(n, g) #part of linear-code decoder + linear_code_decoder=lambda c:decode_linear_code(c, g, syndrome_table) + + #public key: + pub_key=(P_M @ G_M @ S_M)%2 + + + msg_tx='Hello World?' + + blocks_tx=str_to_blocks4(msg_tx) + encrypted=encrypt(blocks_tx, pub_key, t) + + blocks_rx=decrypt(encrypted, P_M_inv, linear_code_decoder, S_inv) + msg_rx=blocks4_to_str(blocks_rx) + + print(f'msg_rx: {msg_rx}') + +
\ No newline at end of file |
