Hilbert Transformation im uC berechnen

[Georg S. (randy)]

Hallo,

ich versuche eine Hilbert-Transformation in einem AVR Mega contoller zu 
berechnen, dazu wollte ich das im Bild dargestellte Flussdiagramm 
umsetzen:

Beitrag "Re: QO-100 und Schmalband-Digimodes"

Erstmal in GNU Octave um zu überprüfen ob das Ergebniss so stimmt. Daran 
scheitert es schon. Da ich von Signalverarbeitung eher keine Ahnung habe 
und einfach nur das Flussdiagramm versucht habe umzusetzen, habe ich 
auch keinen Ansatz wie ich Zwischenergebnisse oder Teilschritte prüfen 
könnte.
Könnt ihr mir das weiter helfen?

Im Skript:

"block2_calc" berechnet einen der rechteckigen Blöcke

"hilbert_algo_1" berechnet den Algo an sich und verwaltet die 
Verzögerungen

das main Programm läd eine wav, berechnet die Hilbert Trafo mit der 
eingebauten Funktion und mit dem eigenen Algorithmus und plottet jeweils 
das Spektrum

[Christoph db1uq K. (christoph_kessler)]

Hallo Georg
Ich habe mit meinem Vorhaben noch nicht weitergemacht, aber die vier 
Allpässe lassen sich ja einzeln testen.
Es gibt auch noch eine andere Veröffentlichung aus der gleichen Zeit vor 
20  Jahren, die acht Allpässe benutzt. Die habe ich in der 
Artikelsammlung 2009 beschrieben und ein Programm auf einem AVR 
angefangen, aber nie zu Ende gebracht.
Schön, dass sich noch jemand für das Thema interessiert.

[Georg S. (randy)]

Hallo Christoph,

vielen Dank für den Tipp. Ich hab einen Block einzeln getestet (hätte 
ich aus selber drauf kommen können), jetzt weiß ich auch dass das ein 
Allpass-Block ist ;-)
Erkenntnis war, dass ich die Variablendefinition der persistenten 
Variablen (=wie static in einer C-Funktion) in Octave so schreiben muss:
  persistent q2a=0;          %  quadratur 0.6628 block
  persistent q2b=0;          %  quadratur 0.6628 block

nicht so:
  persistent q2a=0,q2b=0;    %  quadratur 0.6628 block

Immer diese Details...
Jetzt funktioniert es so wie es soll... zumindest ab ca. 150Hz bei den 
8kHz Abtastrate die ich benutzt habe (siehe Bild, rot ist die FFT der 
Hilbert-Trafo, ca. 30dB Unterdrückung bei 150Hz). Wenn ich das Prinzip 
richtig verstanden habe bedeutet das dass es bei 16kHz Abtastrate des 
Audio dann erst ab 300Hz funktionieren würde?

Und: In der Tat plane ich etwas ähnliches wie du in deinem Beitrag den 
ich verlinkt habe. Es soll ein SSB-Audio-Modulator für 2,4GHz werden um 
ihn für AO-100 Uplink zu benutzen. Mit einem ARM-Blupill-Board, einem 
AD4351-Board und einem digital-Attenuator 0-31dB (alles Aliexpress). Mal 
sehen wie weit ich komme...


>es gibt auch noch eine andere Veröffentlichung aus der gleichen Zeit vor
> 20  Jahren, die acht Allpässe benutzt. Die habe ich in der
> Artikelsammlung 2009 beschrieben

Leider bin ich kein Experte für Signalverarbeitung, so dass ich nur in 
der Lage bin so ein direktes Blockschaltbild nachzuprogrammieren. Mit 
der Beschreibung im pdf alleine hätte ich nichts anfangen könne.
Hast du den Hilbert mit den 8 Allpässen auch so für Laien nachbaubar 
irgendwo erklärt?

Edit: 3x das gleiche Bild weil ich bei editieren nur Bilder hinzufügen 
kann, aber keine löschen...

[Christoph db1uq K. (christoph_kessler)]

Hier die beiden Original-Fundstellen. Der AVR hat ja Multiplizierer, 
wenn auch nur 8 Bit breit. In dieser Literatur geht es speziell um einen 
Aufbau in einem CPLD ohne Multiplizierer, ein Thema, dem sich die beiden 
Professores vom Balkan öfters gewidmet haben. Die Zahlen im 
Blockschaltplan sind noch nicht auf Binärbrüche geändert.

Hatte ich meinen alten Artikel schon verlinkt?
https://www.mikrocontroller.net/articles/Hilbert-Transformator_(Phasenschieber)_mit_ATmega
da ist auch der Link für die Version mit 8 Allpässen
https://web.archive.org/web/20070814013543/http://yehar.com/ViewHome.pl?page=dsp/hilbert/011729.html
da stehen die 8 Parameter. Olli Niemitalo hat sich damals hier im Forum 
auch gemeldet, nachdem ich seinen Namen erwähnt hatte. Google findet 
alles.
Beitrag "Re: Breitband-NF-Phasenschieber, DSP mit AVR"

[Georg S. (randy)]

Im "hilbert_algo_2" habe ich die Koeffizienten für die 8pol IIR Filter 
ausprobiert:
https://www.mikrocontroller.net/attachment/33904/hilbert_Olli_Niemitalo.pdf

Man muss beachten dass die Zahlenwerte der Koeffizienten noch quadriert 
werden müssen, im Gegensatz zum ersten Beispiel. Besonders elegant ist 
der Code nicht, aber das sollte nur ein Test werden ob der Algorithmus 
funktioniert wie gedacht.

Funktioniert bei tiefen Frequenzen fast schon besser als erforderlich 
(das Besipiel im Bild ist auch noch 16kHz, in Gegensatz zum ersten 
Beispiel mit 8KHz Sample Rate). Die zusätzliche Anzahl Multiplikationen 
soll mir aber egal sein, ich habe sowieso vor eine ARM CPU zu benutzen. 
Erstens weil ich keine Lust habe ein Projekt mit einem Prozessor 
anzufangen der von Anfang an schon eher zu klein ist, und zweitens weil 
ich das Projekt als Einstieg in die ARM-Prozessorwelt benutzen will.

Vielen Dank für die Hilfe beim Hilbert.

Das nächste auf dem Programm ist dann der Cordic, da melde ich mich 
vermutlich auch mal im DSP Forum...

[Georg S. (randy)]

Ich hab es jetzt mal in C gegossen. Ist für Prozessoren gedacht bei 
denen ein int 32bit breit ist.
Fall es jemanden interessiert.

// hilbert

// calculates an allpass for the hilbert transform 
// allplass block according to the picture here: https://www.mikrocontroller.net/topic/480404#7103779
// delay : 2-stage delay storage
// coeff : coefficient (16 bits before the decimal, 16 bits after the decimal)
  
int allpass(int x, int coeff, int *delay)
{
  int signal_top_right;
  int y;

  signal_top_right = ( coeff * (delay[0]-x) ) >> 16;
  y                = signal_top_right + delay[0];
  delay[0]         = delay[1];
  delay[1]         = signal_top_right + x;
  
  return y;
}



// Hilbert tranformation for microcontroller according to: https://www.mikrocontroller.net/topic/480404#7103779
// coefficients according to:
// https://www.mikrocontroller.net/topic/92630#839661
// https://web.archive.org/web/20070814013543/http://yehar.com/ViewHome.pl?page=dsp/hilbert/011729.html

void hilbert(int sample_in,int *i_out, int *q_out)
{
  int i1,i2,i3,q1,q2,q3;
    // delay storage
  static int delay_i0;      // inphase Z^-1 block
  static int delay_i1[2]={0,0};  // inphase 1. block
  static int delay_i2[2]={0,0};  // inphase 2. block
  static int delay_i3[2]={0,0};  // inphase 3. block
  static int delay_i4[2]={0,0};  // inphase 4. block
  static int delay_q1[2]={0,0};  // quadrature 1.block
  static int delay_q2[2]={0,0};  // quadrature 2.block
  static int delay_q3[2]={0,0};  // quadrature 3.block
  static int delay_q4[2]={0,0};  // quadrature 4.block


  // inphase 1. to 4. Block
  i1    =allpass(delay_i0,31418,delay_i1);   // coeff = ( ( 0.6923877778065 ) ^ 2 ) << 16
  i2    =allpass(i1,      57434,delay_i2);   // coeff = ( ( 0.9360654322959 ) ^ 2 ) << 16
  i3    =allpass(i2,      64002,delay_i3);   // coeff = ( ( 0.9882295226860 ) ^ 2 ) << 16
  *i_out=allpass(i3,      65372,delay_i4);   // coeff = ( ( 0.9987488452737 ) ^ 2 ) << 16
  // inphase Z^-1 block
  delay_i0=sample_in;
  
  // quadrature 1. to 4. Block
  q1    =allpass(sample_in,10601,delay_q1);   // coeff = ( ( 0.4021921162426 ) ^ 2 ) << 16
  q2    =allpass(q1       ,48040,delay_q2);   // coeff = ( ( 0.8561710882420 ) ^ 2 ) << 16
  q3    =allpass(q2       ,61954,delay_q3);   // coeff = ( ( 0.9722909545651 ) ^ 2 ) << 16
  *q_out=allpass(q3       ,64920,delay_q4);   // coeff = ( ( 0.9952884791278 ) ^ 2 ) << 16
 
}

[Christoph db1uq K. (christoph_kessler)]

>Der mkfisher war ein britischer prof
Das war ja einfach, kleine Internet-Schnitzeljagd mit den ungenauen 
Angaben
Eine Literaturangabe in ISBN 9780080977683

Fisher, T. (2010). Interactive Digital Filter Design. Online Calculator. 
http://www-users.cs.york.ac.uk/~fisher/mkfilter/

Tony Fisher has died, after a short battle with cancer, at the age of 
43. February 29 2000

Auf Github gibt es noch Dateien dazu
https://github.com/university-of-york/cs-www-users-fisher

[Georg S. (randy)]

Eines muss ich noch nachschieben. In der Zeile 13 sollte man die 
Multiplikation explizit auf 64 bit casten, damit man mehr Dynamikbereich 
bekommt.

  signal_top_right = ( (int64_t)coeff * (delay[0]-x) ) >> 16;

[Christoph db1uq K. (christoph_kessler)]

Hallo Georg

Im Webarchiv kann man sich noch Tony Fishers Onlinerechner anschauen:
https://web.archive.org/web/20000815222811/http://www-users.cs.york.ac.uk/~fisher/mkfilter/
Aber rechnen kann man nichts mehr. Es waren anscheinend auch nur 
FIR-Hilbert möglich, kein IIR:

Please choose a filter type from the list:
Infinite Impulse Response (IIR)
* Butterworth  Bessel  Chebyshev
* Resonator
* Proportional-Integral
Finite Impulse Response (FIR)
* Raised Cosine
* Hilbert Transformer

Auf Github gibt es eine Anleitung (6 Seiten) als Postscript-Text. Hier 
in PDF gewandelt.

[Pepe (Gast)]

Der löschmod hat was gegen Hilbert transformationen.

[Christoph db1uq K. (christoph_kessler)]

Die Softwaresammlung "Iowa Hills" hat auch Hilbert-Berechnungen, ist 
aber derzeit "under construction"
die letzte Version im Webarchiv vom Oktober 2021:
https://web.archive.org/web/20210618095158/http://iowahills.com/A5HilbertPhasePage.html

"To start, we show the frequency response for a 65 tap, 90 degree, 
Hilbert Transform Filter."

Es gibt hier noch weitere spezielle Filter, z.B.
"the phase of a 31 tap band pass filter can be adjusted from -135 to 
+135 degrees" (hier als +/- 45 Grad Phasenschieber benutzt, das ergibt 
90 Grad Differenz, ebenfalls zur SSB-Erzeugung)

[Christoph db1uq K. (christoph_kessler)]

>Das nächste auf dem Programm ist dann der Cordic

Andere Baustelle... Auch dazu habe ich mal einen Artikel verfasst:
https://www.mikrocontroller.net/articles/AVR-CORDIC

[Christoph db1uq K. (christoph_kessler)]

>explizit auf 64 bit
Ich habe versucht, auch die Koeffizienten von Olli Niemitalo in 
Binärbrüche umzuwandeln. Dabei sollte man auch möglichst breit rechnen, 
denn die Rundungsfehler wachsen mit jeder Zeile.

Binärbruchzerlegung kann man auch mit dem Taschenrechner. Die Zahl wird 
immer wieder mit 2 multipliziert und wenn sie die Eins überschreitet 
eins abgezogen. Diese Eins ist auch das Ergebnis für die 
Binärbruchstelle, sonst Null.
Ich habe es mit LibreOffice Calc berechnet, aber das benutzt 64 bit 
Gleitkomma. Man sieht ganz unten, dass schon in der siebten oder achten 
Stelle Abweichungen auftreten, das scheint aber erst am Schluß in der 
Kontrollrechnung zu passieren. Der gnome-calculator in Ubuntu rechnet 
mit bis zu 64 bit, da liefert die Kontrollrechnung mit dem berechneten 
Binärbruch auf 12 oder 13 Dezimalstellen identische Ergebnisse.

[ccooss (Gast)]

seine alte Website:

https://web.archive.org/web/20060626003940/http://www.biochem.oulu.fi/~oniemita/dsp/

[Christoph db1uq K. (christoph_kessler)]

von 2003.
Ollis Veröffentlichungen hatte ich schon 2010 im oben genannten Artikel 
verlinkt:
https://www.mikrocontroller.net/articles/Hilbert-Transformator_(Phasenschieber)_mit_ATmega#Scilab-Berechnung_des_Halbbandfilters

https://web.archive.org/web/20070719141658/http://yehar.com/ViewHome.pl?page=dsp/hilbert/
https://web.archive.org/web/20070814013543/http://yehar.com/ViewHome.pl?page=dsp/hilbert/011729.html
http://yehar.com/blog/
http://yehar.com/blog/?p=368
https://dsp.stackexchange.com/questions/37411/iir-hilbert-transformer/59157#59157

In meiner Tabelle habe ich die quadrierten Koeffizienten benutzt, die 
für das Filter benötigt werden, deshalb sind die Zahlen auch so breit. 
Tatsächlich sollten 32 Bit bei weitem ausreichen, das Filter ist mit 16 
Bit Arithmetik vermutlich gut genug.

In der ersten Veröffentlichung vom 4.9.2001 nennt er die Randbedingungen 
seiner Berechnung:
Sampling frequency: 44100Hz. Band: 20-22030Hz. Maximum absolute phase
difference deviation from 90 degrees in that band: 0.7 degrees.
Also sozusagen HiFi-Qualität, die Samplerate einer Audio-CD. Die 22,03 
kHz obere Grenze sind die Nyquist-Frequenz der AD-Wandlung.
Für meine Anwendung im SSB-Sender ist das völlig übertrieben.

[Christoph db1uq K. (christoph_kessler)]

Ich habe noch Ollis (quadrierte) Koeffizienten (auf vier 
Nachkommastellen gerundet) in die Zeichnung des IIR-Hilbert neunter 
Ordnung von Milic/Lutovac eingetragen.

Für die multiplizererlose Version im CPLD werden die durch meine 
Näherungen durch Binärbruchzerlegung ersetzt. Dann sind nur noch 
Schiebungen und Additionen nötig.

[R. M. (n_a_n)]

Ich habe ein wenig experimentiert um die Hilbert Transformation im 
Frequenzbereich zu programmieren.
Es ist als Anregung gedacht, nicht als fertiges Werkzeug.
Ich bin dabei eine DSP Library zu entwerfen.
Als µC verwende ich den P2 von Parallax. Der hat einen CORDIC - Solver 
in Hardware eingebaut und lässt sich in C programmieren.
Anbei das Testergebniss, Codeschnipsel und das Datenblatt.
Vielleicht ist der Beitrag für jemanden nützlich.

[T.U.Darmstadt (Gast)]

Christoph db1uq K. schrieb:
> Ich habe noch Ollis (quadrierte) Koeffizienten (auf vier
> Nachkommastellen gerundet)

Sind die nur in der o.g. Darstellung gerundet oder auch in der 
Applikation?

[Christoph db1uq K. (christoph_kessler)]

Oben hatte ich die Originalfundstellen im Webarchiv verlinkt.
Die Zahlen dort müssen noch quadriert werden, damit werden sie noch 
länger:

"The a coefficients are:
filter1:
 0.6923877778065, 0.9360654322959, 0.9882295226860, 0.9987488452737
filter2:
 0.4021921162426, 0.8561710882420, 0.9722909545651, 0.9952884791278"

Der Gnu-Taschenrechner sagt beispielsweise zum ersten Faktor:
0,6923877778065 * 0,6923877778065 = 0,47940083485582321395144225

[He. (Gast)]

Christoph db1uq K. schrieb:
> ür die multiplizererlose Version im CPLD werden die durch meine
> Näherungen durch Binärbruchzerlegung ersetzt. Dann sind nur noch
> Schiebungen und Additionen nötig.

Ist das noch zeitgemäß? Moderne FPGA-Typen werfen doch mit 
Multiplizieren nur so um sich.

[Christoph db1uq K. (christoph_kessler)]

Ja man kann sich auch Weltschmerz als Hobby wählen. Alles doch schon 
dagewesen. Nichts neues unter der Sonne. Kann man doch alles fertig 
kaufen.

Ich habe hier ein paar Schrottplatinen mit einem CPLD drauf, die kann 
ich auslöten. Aus meinem Bohrständer und einem Heißluftfön habe ich eine 
Entlötstation gebastelt. Dazu mehrere Laborplatinen für TQFP144 vom 
Flohmarkt Friedrichshafen. Die aktuelle Quartus Version installiert, 
jetzt muss ich mich nur wieder an Verilog erinnern, ist schon ein paar 
Jahre her.
FPGAs müsste ich kaufen, damit habe ich noch praktisch keine Erfahrung 
(zwei Eval-Boards liegen schon lange rum). Einen Byteblaster habe ich 
auch um die .pof-Datei einzuflashen. 16-bit-ADC und DAC habe ich auch, 
dazu muss ich dann auch eine SPI-Bus-Logik in Verilog einbauen.

[-gb- (Gast)]

Christoph db1uq K. schrieb:
> 0,6923877778065 * 0,6923877778065 = 0,47940083485582321395144225

Damit werden aber hinten einige Stellen dazu erfunden. Grundsätzlich 
wird man nicht genauer, als die Ausgangswerte, hier also ca 13 Stellen 
-> 12.5.
Durch den forgesetzten Fehler sind es sogar nur 12.25 und weil 0.47 nur 
70% von 0.69 sind, bleiben nur knapp 12 Stellen!