Forum: Digitale Signalverarbeitung / DSP / Machine Learning Hilbert Transformation im uC berechnen


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
von Georg S. (randy)


Angehängte Dateien:

Lesenswert?

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

von Christoph db1uq K. (christoph_kessler)


Lesenswert?

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.

von Georg S. (randy)


Angehängte Dateien:

Lesenswert?

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...

: Bearbeitet durch User
von Christoph db1uq K. (christoph_kessler)


Angehängte Dateien:

Lesenswert?

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"

: Bearbeitet durch User
von Georg S. (randy)


Angehängte Dateien:

Lesenswert?

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...

von Georg S. (randy)


Lesenswert?

Ich hab es jetzt mal in C gegossen. Ist für Prozessoren gedacht bei 
denen ein int 32bit breit ist.
Fall es jemanden interessiert.
1
// hilbert
2
3
// calculates an allpass for the hilbert transform 
4
// allplass block according to the picture here: https://www.mikrocontroller.net/topic/480404#7103779
5
// delay : 2-stage delay storage
6
// coeff : coefficient (16 bits before the decimal, 16 bits after the decimal)
7
  
8
int allpass(int x, int coeff, int *delay)
9
{
10
  int signal_top_right;
11
  int y;
12
13
  signal_top_right = ( coeff * (delay[0]-x) ) >> 16;
14
  y                = signal_top_right + delay[0];
15
  delay[0]         = delay[1];
16
  delay[1]         = signal_top_right + x;
17
  
18
  return y;
19
}
20
21
22
23
// Hilbert tranformation for microcontroller according to: https://www.mikrocontroller.net/topic/480404#7103779
24
// coefficients according to:
25
// https://www.mikrocontroller.net/topic/92630#839661
26
// https://web.archive.org/web/20070814013543/http://yehar.com/ViewHome.pl?page=dsp/hilbert/011729.html
27
28
void hilbert(int sample_in,int *i_out, int *q_out)
29
{
30
  int i1,i2,i3,q1,q2,q3;
31
    // delay storage
32
  static int delay_i0;      // inphase Z^-1 block
33
  static int delay_i1[2]={0,0};  // inphase 1. block
34
  static int delay_i2[2]={0,0};  // inphase 2. block
35
  static int delay_i3[2]={0,0};  // inphase 3. block
36
  static int delay_i4[2]={0,0};  // inphase 4. block
37
  static int delay_q1[2]={0,0};  // quadrature 1.block
38
  static int delay_q2[2]={0,0};  // quadrature 2.block
39
  static int delay_q3[2]={0,0};  // quadrature 3.block
40
  static int delay_q4[2]={0,0};  // quadrature 4.block
41
42
43
  // inphase 1. to 4. Block
44
  i1    =allpass(delay_i0,31418,delay_i1);   // coeff = ( ( 0.6923877778065 ) ^ 2 ) << 16
45
  i2    =allpass(i1,      57434,delay_i2);   // coeff = ( ( 0.9360654322959 ) ^ 2 ) << 16
46
  i3    =allpass(i2,      64002,delay_i3);   // coeff = ( ( 0.9882295226860 ) ^ 2 ) << 16
47
  *i_out=allpass(i3,      65372,delay_i4);   // coeff = ( ( 0.9987488452737 ) ^ 2 ) << 16
48
  // inphase Z^-1 block
49
  delay_i0=sample_in;
50
  
51
  // quadrature 1. to 4. Block
52
  q1    =allpass(sample_in,10601,delay_q1);   // coeff = ( ( 0.4021921162426 ) ^ 2 ) << 16
53
  q2    =allpass(q1       ,48040,delay_q2);   // coeff = ( ( 0.8561710882420 ) ^ 2 ) << 16
54
  q3    =allpass(q2       ,61954,delay_q3);   // coeff = ( ( 0.9722909545651 ) ^ 2 ) << 16
55
  *q_out=allpass(q3       ,64920,delay_q4);   // coeff = ( ( 0.9952884791278 ) ^ 2 ) << 16
56
 
57
}

Beitrag #7153423 wurde von einem Moderator gelöscht.
von Christoph db1uq K. (christoph_kessler)


Angehängte Dateien:

Lesenswert?

>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

: Bearbeitet durch User
von Georg S. (randy)


Lesenswert?

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;

: Bearbeitet durch User
von Christoph db1uq K. (christoph_kessler)


Angehängte Dateien:

Lesenswert?

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.

: Bearbeitet durch User
von Pepe (Gast)


Angehängte Dateien:

Lesenswert?

Der löschmod hat was gegen Hilbert transformationen.

von Christoph db1uq K. (christoph_kessler)


Lesenswert?

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)

: Bearbeitet durch User
von Christoph db1uq K. (christoph_kessler)


Lesenswert?

>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

von Christoph db1uq K. (christoph_kessler)


Angehängte Dateien:

Lesenswert?

>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.

von ccooss (Gast)


Lesenswert?


von Christoph db1uq K. (christoph_kessler)


Lesenswert?

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.

von Christoph db1uq K. (christoph_kessler)


Angehängte Dateien:

Lesenswert?

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.

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.