Integer (1 Byte) Datentransformation DCT (FFT) gesucht

[Möbius]

Guten Tag Forenteilnehmer

Dies ist mein erster Beitrag in diesem Forum überhaupt und ich begrüsse Euch herzlich. Ehrlich gesagt mein erster Forenversuch überhaupt.https://www.delphipraxis.net/images/.../icon_cool.gif
Das Forum ist gut und hat mir schon viel geholfen. Hoffentlich kann ich anderen auch mal helfen.

Das Problem. Ich bin am Verzweifeln. Schon drei (unbrauchbare) Bücher gewälzt aber noch keine Lösung. Immer ist die Rede von Butterfly-Diagrammen (also Zerteilungen der Daten) und irgendwelchen "Lifting up Faktoren" etc.

Ich bräuchte für ein Projekt Delphi/Pascal Quellcode der mir für eine beliebige Datenmenge (kann auch auf Potenzen von 2, verfüllt mit restlichen Nullen sein) eine Integer Cosinus Transformation durchführt.
Auch eine Quelle mit lesbarem/übersetzbaren Pseudocode oder so wäre hilfreich.
Nach tagelanger Suche im Netz stosse ich immer nur auf reelwertige Funktionen (DCT) bzw. komplexwertige Funktionen (FFT) mit Floating Point Zahlen. Habe ich selber schon implementiert, runden nach der Operation führt aber bei der Rücktransformation zu Fehlern.

Dies wie schon im Titel vermerkt für Zeichengrössen eines Bytes, jedoch nicht für irgendwelche 8x8 oder 16x16 oder sonstige Matrizen sondern, mehr oder weniger, beliebige eingegebene Datengrössen.

Also der Wunsch ist eigentlich eine Datenmenge vom Zeitbereich in den Frequenzbereich zu transformieren. Dies auf Integer (Byte) Basis welche auch ohne Fehler Rücktransformierbar ist.

Eine passende Integer FFT ist auch okay, sofern im Übetragungskanal nur der reelwertige Anteil übergeben werden muss und die Rücktransformation auch mit einem (imaginären) Phasenanteil von Null klappt. Daher denke ich aber eine Cosinus Transformation eignet sich besser da kein komplewertiger Anteil.

Kennt jemand von Euch so eine Transformation und Rücktransformation oder hat sie schon selber implementiert?

Hoffe hab mich verständlich ausgedrückt
Bin für jede Hilfe dankbar

Reto

[TiGü]

Sowas?

https://github.com/jcodec/jcodec/blo...ct/IntDCT.java
https://javadoc.io/static/org.jcodec...ct/IntDCT.html

https://github.com/libjpeg-turbo/lib...ain/jfdctint.c
https://github.com/libjpeg-turbo/lib...ain/jidctint.c

[Möbius]

Hi TiGü

Ganz herzlichen Dank.
Nach sowas habe ich tagelang gesucht.
Reto

[TiGü]

Darf ich fragen, wie und womit du gesucht hast?

[Möbius]

Ja sicher gerne
Natürlich mit Google.
Dabei bin ich auf einige Bücher gestossen welche ich z.T. bestellt habe.
-> Besslich diskrete Orthogonaltransformationen.
-> K. R. Rao Fast Fourier Transform: Algorithms and Applications.

Das Problem bei den Büchern ist, dass sie sehr viel mathematischen Background enthalten aber wenn es um Integer Implementationen geht dann wird nur kurz von irgendeinem "Leverage Faktor" gesprochen. Zudem ist in den sog. "Butterfly" Diagrammen nicht ganz klar wie eine beliebige Datenmenge bearbeitet werden kann. Also auf eine Konstruktion eines Algorhitmuses in einer beliebigen Sprache gehen die Bücher nicht ein. Schon Pseudocode würde helfen.

Weitere Suche im Internet hat ein bisschen Code z.B. bei Rosetta Code zu Tage gefördert. Jedoch nur teilweise brauchbar. Auch in den Büchern steht, dass die DCT über den Umweg einer DFT gemacht werden kann. Aber eben auch wenn Integer ist diese komplexwertig und es ist mir nicht gelungen ohne den Imaginärteil eine Rücktransformation zu machen. Ist ja auch irgendwie klar. Der Realteil einer DFT stellt das Frequenzspektrum dar während der Imaginärteil die Phasenlagen des Signals darstellt.

Tja Deine Hilfestellung hat Hoffnungen geweckt, aber die Codes übersetzt in Delphi brauchen ca. 10 Minuten für 32 KiloByte :-(.
Es ist eben auch nur eine DCT und keine FCT.
Die Bücher sagen aber eine solche soll möglich sein.

Was ich wirklich brauche ist ein Algorithmus der mir ein Signal (Daten beliebiger Grösse N -> also wenn möglich ohne Nullen die auf 2^n verfüllen) vom Zeitbereich in den Frequenzbereich transformiert und wieder zurück. Dies aber nur reelwertig und auf Integer (Bytes) basierend.

Ach ja und mit den Codes aus dem Netz stimmt etwas nicht. Weil sie scheinen kein erwartetes Spektrum darzustellen.

Um was geht es überhaupt:
Ich arbeite an einem Daten-Kompressionsalgorithmus. Der ist schon sehr weit fortgeschritten und sieht derzeit mit den Daten aus dem ->Calgary Corpus bzw. dem ->Canterbury Corpus oder auch dem ->Large Corpus sehr gut aus.
Jedoch wie üblich mit relativ unstrukturierten Daten mit hoher Entropie (also auf Deutsch gesagt eine Sauordnung), da schneidet er permanent mit rund 100% (und nur sehr wenig mehr bei moderat grossen Files) ab! Wohlgemerkt gemessen an RAD-Daten (randomisierten Files).
Die Sachlage ist klar, wenn ich das nur um ein paar % runterdrücken kann, dann habe ich einen Kompressor der sich selber komprimieren kann! Die Konsequenzen bei einer iterativen Anwendung dürften klar sein.
Somit dürfte schon klar sein, dass mir eine FFT sehr wenig nützt. Ob nun integer oder nicht. Der Frequenzanteil (reel) (bei meiner eigenen Float64 basierten FFT) sieht wir erwartet gut aus bei solchen Daten. Ich bin ja von Haus aus Elektroniker und Toningenieur (ist schon eine Weile her ;-)). Aber eben solche RAD-Daten bezeichnet man im Audio-Jargon als statistisches weisses Rauschen. Und dieses Zeichnet sich aus durch sehr kleine Schwankungen der Amplituden (infinitesimal wohl gar keine) im Spektrum.

Und genau diese kleinen Schwankungen könnten dazu führen, dass mein Algorithmus die Daten besser verarbeiten kann.

Wie gesagt. Ist reine Theorie und Forschungsarbeit.

Ach ja, die DCT aus dem Netz von Deinem Tipp rechnet nicht nur ewig (eine FCT wäre natürlich besser), Sie zeigt auch ein Spektrum von sehr stark schwankenden Amplituden der Frequenanteile - Da kann etwas nicht stimmen. Die müssten eigentlich näher beisamen liegen.
Ganz klar ist, dass am Anfang des Spektrums sehr hohe Amplituden auftreten wenn ich die Daten als rein positive Werte (Bytes) einspeise. Dieser DC-Anteil schlägt hier natürlich zu Buche. Wenn ich die Daten als Int8 einspeise sieht das schon besser aus, die realtiv hohen Werte am Anfang der Sequenz dürften Verzerrungen sein, die sich ergeben weil ich kein kontinuierliches Signal einspeise. Damit liesse sich leben wenn die Rücktransformation diese Verzerrungen richtig heraus rechnet.

So suche ich also weiter...

Hoffe Dir aber den Werdegang meiner Idee und Suche etwas näher gebracht zu haben.

[TiGü]

Aber es ist doch auch möglich, dass du dich beim Übersetzen vertan hast?!

Welches Codebeispiel hast du denn nach Delphi übersetzt?
Das Java- oder das C-Beispiel?!

[Möbius]

Ja könnte natürlich auch sein.

Hier die Übersetzungen (non integer):

Delphi-Quellcode:

class procedure TCompressOSCT.DCT(const Input: TFloat64Array; out Output: TFloat64Array);  //1 double[] DCT (double[] g) { // forward DCT on signal g
var
  M, mCnt, uCnt: Int64;
  s, cm, sums, Phi: Float64;
begin
  M := Length(Input); //2 int M = g.length;
  s := sqrt(2.0 / M); //3 double s = Math.sqrt(2.0 / M); // common scale factor
  SetLength(Output, M); //4 double[] G = new double[M];
  for mCnt := 0 to High(Input) do //5 for (int m = 0; m < M; m++) {
  begin
    cm := 1; // 6 double cm = 1.0;
    if mCnt = 0 then // 7 if (m == 0)
    begin
      cm := 1 / sqrt(2); // 8 cm = 1.0 / Math.sqrt(2);
    end;
    sums := 0; // 9 double sum = 0;
    for uCnt := 0 to High(Input) do // 10 for (int u = 0; u < M; u++) {
    begin
      Phi := Pi * mCnt * (2 * uCnt + 1) / (2 * M); //11 double Phi = Math.PI * m * (2 * u + 1) / (2 * M);
      sums := sums + (Input[uCnt] * cm * cos(Phi)); // 12 sum += g[u] * cm * Math.cos(Phi); }
    end;
    Output[mCnt] := s * sums; // 14 G[m] = s * sum; }
  end;
end;

class procedure TCompressOSCT.InvDCT(const Input: TFloat64Array; out Output: TFloat64Array);  // 20 double[] iDCT (double[] G) { // inverse DCT on spectrum G
var
  M, mCnt, uCnt: Int64;
  s, cm, sums, Phi: Float64;
begin
  M := Length(Input); // 21 int M = G.length;
  s := sqrt(2.0 / M); // 22 double s = Math.sqrt(2.0 / M); //common scale factor
  SetLength(Output, M); // 23 double[] g = new double[M];
  for uCnt := 0 to High(Input) do // 24 for (int u = 0; u < M; u++) {
  begin
    sums := 0; // 25 double sum = 0;
    for mCnt := 0 to High(Input) do //26 for (int m = 0; m < M; m++) {
    begin
      cm := 1; // 27 double cm = 1.0;
      if mCnt = 0 then // 28 if (m == 0)
      begin
        cm := 1 / sqrt(2); // 29 cm = 1.0 / Math.sqrt(2);
      end;
      Phi := Pi * mCnt * (2 * uCnt + 1) / (2 * M); // 30 double Phi = Math.PI * m * (2 * u + 1) / (2 * M);
      sums := sums + (Input[mCnt] * cm * cos(Phi)); // 31 sum += G[m] * cm * Math.cos(Phi); 32 }
    end;
    Output[uCnt] := s * sums; // 33 g[u] = s * sum; 34 }
  end;
end;

[freimatz]

"0 to High(Input)" - müsste das nicht -1 sein? Passt vermutlich. "Low(Input) to High(Input)" wäre schöner.

[TiGü]

Warum zwei geschachtelte Schleifen?
In der Wikipediadefintion von DCT-I bis DCT-IV sehe ich immer nur ein Summenzeichen?

[Blup]

Ich habe mal versucht zu optimieren.
Unnötige doppelte Operationen mit Hilfe von Zwischenvariablen aus Schleifen ziehen.
Multiplikation durch Addition ersetzen, wenn das möglich ist.
Ob das eine merkliche Verbesserung der Laufzeit bringt, müsste man in der Praxis testen.

Delphi-Quellcode:

class procedure TCompressOSCT.DCT(const Input: TFloat64Array; out Output: TFloat64Array);
var
  M, mCnt, uCnt: Int64;
  s, sums, Phi, Pi2M, Pi2MmCnt: Float64;
begin
  M := Length(Input);
  SetLength(Output, M);
  if M > 0 then
  begin
    s := sqrt(2.0 / M);
    Pi2M := Pi / (2 * M);
    Pi2MmCnt := 0;
    for mCnt := 0 to M - 1 do
    begin
      sums := 0;
      Phi  := Pi2MmCnt;
      for uCnt := 0 to M - 1 do
      begin
        sums := sums + (Input[uCnt] * cos(Phi));
        Phi := Phi + Pi2MmCnt + Pi2MmCnt;
      end;
      Output[mCnt] := s * sums;
      Pi2MmCnt := Pi2MmCnt + Pi2M;
    end;
    Output[0] := Output[0] / sqrt(2);
  end;
end;

Warum wird eigentlich ausgerechnet nur Output[0] durch sqrt(2) dividiert?