Normalerweise nimmt man ja die euklidische Distanz / Metrik. Bei der DFT hat man einen Vektor von N Werten, und man bekommt N Frequenzamplituden & Phasen, für Frequenzen 0,1,2,3,...,N-1 als Vektor heraus. Das ist auch wieder umkehrbar. Bei der normalen DFT geht das ja so:
Normale Diskrete Fourier Transformation
Wenn ich bei einem Punkt die Distanz zum Nullpunkt haben will, dann ist das sqrt(x^2+y^2).
Bei sqrt(x^2+y^2)=1 hat man einen Einheitskreis in der Euklidischen Metrik.
Mit sin(w) / cos(w) kriegt man die x y Koordinaten von einem Punkt, angefangen mittig oben, und wenn man w linear erhöht, geht der Punkt linear entlang dem Kreis herum.
Mit atan2(y,x) bekommt man den Winkel wieder.
Wir haben also diese N Eingabe Werte.
Und von jeder der Frequenzen nehmen wir N sample, mit sin() und cos() können wir die erzeugen:
Wir haben also eine Matrix mit den Wellenformen, und einen Vektor mit den Eingabewerten. Die multipliziert man einfach, und man bekommt zurück, wie gut die korrelieren. Jeweils für Sin und Cos, die sind ja 90deg Phasenversetzt.
Bei diesen 2 Vektoren muss man nur noch die Distanz (sqrt(x^2+y^2)) und die Phase (atan2(y,x)) ermitteln, x/y sind jeweils die korrespondierenden Werte dieser Vektoren.
Und schon hat man die gesuchten Frequenzanteile und Phasen.
Das Umkehren ist dann ja einfach. Man nimmt N Samples von jeder Frequenz, mit der jeweiligen Phase, mit sin, multipliziert es mit der Amplitude, und addiert die Resultate der unterschiedlichen Frequenzen zusammen.
Das unschöne an euklidisch ist halt, mit sin / cos bekommt man irrationale Zahlen da braucht man Floats und verliert Genauigkeit. Darum dachte ich mir:
Diskrete Fourier Transformation mit Taxicab Metrik
Das vorgehen müsste ja gleich sein. Aber statt dass die Distanz durch sqrt(x^2+y^2) definiert ist, ist es in der Taxicab Metrik |x|+|y|.
Ein Einheitskreis ja durch eine Distanz zum Zentrum von 1 definiert.
|x|+|y|=1 ergibt ein um 45deg gedrehtes Quadrat. In der Taxicab Metrik ist das der Einheitskreis.
Statt Sin / Cos hat man dann eine Dreieck Welle. Das kommt dabei raus, wenn man einem Punkt dem "Kreis" entlang folgt.
Die Funktion, um aus X/Y wieder den Winkel zu kriegen, muss auch angepasst werden.
Alles andere müsste ja gleich bleiben.
Das Problem
Also habe ich das mal Ausprobiert (siehe Anhang). DFT, und wieder umkehren. Ich dachte schon, es hätte Funktioniert. 1 beliebiger Wert geht. 2 beliebige Werte gehen. 4 beliebige Werte gehen. Einfache Sachen wie 5 4 3 4 5 4 3 4 gehen. Sachen wie 5 4 3 4 5 4 7 4 auch.
Aber bei so ziemlich allen anderen Situationen, liege ich mit Taxicab doch etwas daneben... (euklidisch geht einwandfrei)
Woran könnte das liegen? Mache ich was falsch? Geht das einfach nicht? Sind meine floats nicht genau genug?
Wozu das ganze
Bei Taxicab habe ich nur Additionen, Multiplikationen, und Lineare Funktionen. Es müsste eigentlich machbar sein, sich auszurechnen, wie Gross ein integer sein muss, um alle Zustände abzubilden und die müssten da dann ja alle genau da rein passen. Dann bräuchte ich keine Floats mehr, und könnte die Daten verlustfrei übertragen.
Ich wollte das auf ein Bild anwenden. Die Idee war, statt wie bei anderen Bildformaten nur Blöcke begrenzter Grösse so umzuwandeln, das mit dem ganzen Bild zu machen, und zwar verlustfrei.
Und das wiederum wollte ich für das progressive Laden von Bildern in einer Anwendung verwenden. Die dft Daten hätte ich deflate komprimiert übertragen. Es hätte viele Bilder geben können, deren Daten ich abwechselnd in kleinen Stückchen über die selbe Verbindung übertragen hätte. Das hätte ich priorisiert, je nachdem ob die Bilder in der nähe des Sichtbereichs sind oder nicht.
Dann würden die Bilder währen dem Laden immer schärfer, biss ich das Original rekonstruiert hätte.
Bei Webseiten werden normalerwiese Tricks verwendet. z.B. ein "Blur Hash" (also ein weiteres, kleines, unscharfes vorgeladenes Bild). Ok z.B. JPEG kann progressives laden, aber mit den Blöcken sieht es nicht so schön aus. Irgendwie eckig.
Mit dem Ansatz wäre das ohne Hacks gegangen...
Kontinuierliche und diskrete Fourier-Transformationen beruhen beide auf Skalarprodukten. (gewöhnliches) Skalarprodukt eines Vektors mit sich selbst ergibt - euklidische Norm zum Quadrat. D. h. in dem Fall stammt die Metrik von einem Skalarprodukt. Die 1-Norm (|x|+|y|) kann man dagegen nicht durch ein Skalarprodukt erhalten, denn die erfüllt die sog. Parallelogrammidentität nicht.
"Die multipliziert man einfach, und man bekommt zurück, wie gut die korrelieren." Naja, das ist halt nix anderes als ein Skalarprodukt zu bilden ...
Ich habe leider gerade wenig Zeit, daher kann ich nicht auf deine DFT-Fragestellung eingehen. Eine Sache fällt mir aber direkt auf bzw. ein: Warum nutzt du keine DWT mit einem Haar-Wavelet? Wenn ich mich nicht täusche, dann kriegst du da mehr oder minder das, was du willst.
Ich muss mir das mit der Parallelogrammidentität usw. mal genauer anschauen, und das alles nochmal schritt für schritt Durchgehen.
Von DWT + Haar-Wavelet hatte ich bisher noch nichts gehört. Sieht vielversprechend aus, ich muss mir mal anschauen, ob das für mein Vorhaben geeignet ist.
Bezüglich Lambda-Kalkül... Naja, ich weiss ja nicht. Ich glaube das wäre mir zu umständlich. Da probiere ich lieber einen theorem prover wie z.B. Coq.
Nachdem ich diverses über Orthogonale Funktionen, Matrizen, etc. gelesen habe, ist mir einiges klar geworden. Mit einer DFT Matrix, ist die transponiert-konjugierte Matrix die Inverse der Matrix (zumindest, wenn man sie richtig skaliert).
Eine Matrix mal ihre Inverse gibt eine Identity Matrix.
Die transponiert-konjugierte Matrix der Dreiecksfunktion war hingegen nicht die Inverse der Matrix, also hat es sich nicht mehr aufgehoben.
Soweit ich das bisher ausprobiert habe, hat die Matrix aus der Dreiecksfunktion aber in der Regel eine Inverse, kann also umgekehrt werden. Ich hoffe mal dass sie immer eine hat, aber ich wüsste nicht, warum das nicht der Fall sein sollte. Sieht so aus, als bestünde die Inverse Matrix aus Rationalen Zahlen, aber ich muss noch genauer nachsehen, welche Eigenschaften diese hat, und wie diese zustande kommt.
Ich bin auch noch auf die Hadamard Matrizen / Transformationen gestossen. Die haben einige interessante Eigenschaften. Bestehen nur aus 1 und -1, sind ihre eigene Inverse (abgesehen von der Skalierung). Und richtig sortiert hat es ähnliche Eigenschaften wie eine Fourier Transform. Es gibt sogar auch eine FHT ähnlich wie bei der FFT.
Ich habe mit ein paar Python Skripten mal etwas ausprobiert, um zu sehen, wie gut diese geeignet sind, um ein Bild mit zunehmender Auflösung zu übertragen, und ob die Unschärfe dabei auch gut aussieht. Ich habe einen FFT shift direkt auf die Matrizen angewendet, in der Mitte ein Rechteck, Kreis, oder Rhombus ausgeschnitten, und dann mit der Inversen Matrix Multipliziert.
Hadamard Pixelliert das Bild, und man sieht am Wenigsten mit der Transformation.
Fourier sieht, wie zu erwarten, am besten aus.
Dreiecksfunktion hat ein paar Streifen drin, sieht aber eigentlich auch sehr gut aus. Muss ich definitiv genauer anschauen, die kann ich sicher so skalieren, dass ich nur Integer in der Matrize und der Inversen habe.
Im Anhang noch die Skripte und ein paar Testbilder zum Vergleich. Das s<nummer> gibt die Anzahl Samples an, die genutzt wurden.
Noch ein weiterer Vergleich. Bei Interlacing für progressives Laden wird ja jedes n-te Pixel verwendet. Also habe ich mal 32x32 Pixel vom Bild oben genommen, und wieder auf 512 Pixel vergrössert. Beim vergrössern, einmal nearest (Pixelliert), und einmal Kubisch.
Ich würde sagen, besser als Hadamard, aber weniger gut als die anderen DFTs.
Sieht so aus, als bestünde die
Inverse Matrix aus Rationalen Zahlen, aber ich muss noch genauer
nachsehen, welche Eigenschaften diese hat, und wie diese zustande kommt.
Ich habe gerade mal nachgesehen, wie Gross die DFT Matrix Werte mit der Dreiecksfunktion würden, wenn ich diese mit einem Faktor multipliziere, so dass ich nur Integer habe. Die grösste Zahl einer N x N Matrize wäre A003418(N), das least common multiple aller Zahlen bis und mit N.
Für das 512 x 512 Bild, wäre das 37502641566909080909774372402973685284974612296925670251719288153975177600052224658419106726814633548114580543484744358444505293773294084252491062491589656945636456083649455529431070527198260042812115280192232693996980672000.
Der Ansatz ist also nicht besonders praktikabel.
Bei Hadamard wäre es zwar machbar, aber da das Ergebnis schlechter als bei simplem Interlacing wäre, lohnt es sich nicht.
Eventuell könnte ich eine normale FFT nehmen, und am Ende die Differenz zum Original berechnen, und wieder hinzufügen. Aber float Berechnungen sind nicht besonders deterministisch. Selbst wenn 2 Platformen beide IEEE754 floats verwenden, kann es geringe Unterschiede bei Berechnungen geben, besonders auch bei Grafikkarten, das kann ich hier nicht brauchen...
Edit: Ich habe gerade von der Number-theoretic transform (NTT) gelesen. Die muss ich mal ausprobieren, vielleicht ist das ja, wonach ich gesucht habe.
Aber ich weiss jetzt, wie ich es wohl am besten umsetze. Ich nehme eine gewöhnliche DFT, aber verwende Fixed Point statt Floating Point für die Berechnungen. Damit habe ich die Genauigkeit bei jedem Schritt ja genau im griff, und müsste die so wählen können, dass am Ende alles wieder exakt passt. Unterschiede zwischen Platformen gibt es auch keine Relevanten.
Manchmal sind die einfachsten & offensichtlichsten Lösungen die besten...
Damit habe ich die Genauigkeit bei jedem Schritt ja genau
im griff,
Ja genau, die FFTs rechnest Du mit integer und führst eine Skalierung mit. Bei einer 512er FFT hast Du log2(512)=9 stages, jedesmal kuckst Du auf den Wertebereich und skalierst ggf. runter. Modernere Proz. benötigen für 32Bit float aber genauso lange wie für 32 Bit integer, die Zeite wo man mit integer sparen konnte sind imho vorbei.
Noch ein Vergleich, diesmal DCT-IV statt eine DFT. Da gibt es einige interessante Unterschiede. Bei allen anderen Varianten habe ich einen FFT-Shift gemacht, und einen teil in der Mitte verwendet. Bei der DCT kann man hingegen einfach den Teil Oben Links ausschneiden, ohne FFT-shift. Hier hat man auch keine Komplexe Zahlen. Qualität ist vergleichbar mit der DFT, abgesehen vom schwarzen Streifen, den man unten und rechts bekommt.
Heute habe ich von der DHT (Discrete Hartley transform) gelesen. Resultate sehen genau gleich aus wie bei der DFT (macht ja auch Sinn, da DHT(X)*(1+i)=DFT(x) ). Aber man hat keine komplexen Zahlen.
Modernere Proz.
benötigen für 32Bit float aber genauso lange wie für 32 Bit integer, die
Zeite wo man mit integer sparen konnte sind imho vorbei.
Dann hast du aber die Möglichkeit unterschlagen, diese Operation zu parallelisieren, bspw. in FPGAs. Und versuche dich mal mit Floats im FPGA ...
Ich hätte ansonsten auch vorgeschlagen, es mit Wavelets zu versuchen. Bei der Transformation kann man statt mit Tap-Filtern (immer verlustbehaftet) auch ein digitales, verlustfreies Verfahren wählen. Werden nicht alle Planes codiert/decodiert ist bei der eigentlich verlustfreien Variante allerdings das S/N schlechter als bei den Tap-Filtern. Letztere sind in der Performance deutlich überlegen. Bei den hier gezeigten Bildern ist aber Verlustfreiheit kein Thema.
Vor etlichen Jahren hatten wir für Space-Anwendungen das SPIHT-Verfahren umgesetzt. Zum Einsatz kamen DSP mit 40 Bit, von denen wir aber nur 32 Bit genutzt hatten. Das Float-Format hatte dann nur 23 oder 24 Bit in der Mantisse, war für diese Anwendung aber seinerzeit das Optimum, was man unter Weltraumbedingungen (Leistung/Strahlung) machen konnte. Bildgrößen waren 2k x 2k.
Ich habe jetzt mal eine simple DHT in JS implementiert: https://temp-s.s.abrecht.li/ft/ft2.html (Ist aber nur für kleine Bilder geeignet, und geht nur in chrome).
Ich habe aber die Reihenfolge der Zeilen und Spalten der DHT Matrix etwas angepasst. Im python Program weiter oben hatte ich die einfach um die hälfter der Bildbreite verschoben (fftshift), so das die niederfrequentesten Anteile in der Mitte sind. Aber hier nehme ich stattdessen 0 1 3 5 7 ... 6 4 2. Dadurch sind die Anteile oben links, statt in der Mitte, wie bei einer DCT.
In der web Demo hab ich noch einen Slider eingebaut, wo man einstellen kann, wie viele Werte der DHT bei der Umkehrung verwendet werden. Ist schon erstaunlich, wie wenig Daten es braucht, bis man was erkennen kann.
Eventuell schaue ich mir Wavelets später nochmal an. Das SPIHT zeug tönt auch interessant, aber da scheinen Patente drauf zu sein, sind die noch gültig?
Ich bin auch noch auf ein paar andere Sachen gestossen, die ich mir noch genauer ansehen will.
Offenbar hatten auch schon andere die Idee mit einer Integer Version von DFT, DHT, usw. Da habe ich ein Paper dazu gefunden "The integer transforms analogous to discrete trigonometric transforms" (DOI: 10.1109/78.886998). Mal sehen, ob ich damit was anfangen kann.
Die DTFT & co. will ich mir auch mal noch anschauen.
Ich muss mir auch noch überlegen, was ich mit wirklich grossen Bildern mache.
Große Bilder kannst du splitten (Tiles oder Slices) mit etwas Überlappung. Die Ränder werden für die Faltung ohnehin entsprechend gespiegelt, um Artefakte zu reduzieren.
