Abgetastete Sinuskurve filtern

[Third E. (third-eye)]

Hallo zusammen,

ich taste mit dem ADC eines 8Bit-µCs jeweils eine komplette Periode 
einer Sinus-Kurve ab und speichere die Werte in einem Array mit 80 
Elementen.
Es wird immer eine komplette Periode gespeichert, allerdings ist die 
Amplitude und die Phase beliebig. Die Frequenz hingegen ist konstant und 
proportional zur Oszillatorfrequenz des µCs, da er selber den Sinus 
erzeugt.
Leider sind auch ab und zu Störspitzen/Ausreißer in den Abtastungen 
enthalten. Um die Ausreißer gering zu halten, summiere ich 64 Perioden 
in dem Array auf und bilde Anschließend pro Arrayelement den Mittelwert. 
Die Zahl 64 ergibt sich übrigens daraus, dass max. 64 10Bit-Werte in ein 
u16 passen.
Leider kann es auch Störungen geben, die immer zum selben Phasenwinkel 
auftreten. In dem Fall bringt mir auch die Mittelung nichts.
Hat jemand eine Idee für einen robusten Algorithmus, wie man diese 
Ausreißer trotzdem los wird?

Ein Software-Tiefpass hoher Ordnung würde wohl helfen. Allerdings traue 
ich mir das nicht wirklich zu, so ein Filter selber in C zu 
implementieren.

Im Artikel https://www.mikrocontroller.net/articles/Filter steht 
folgendes:
"Im einfachsten Fall multipliziert man einen eingehenden Datenstrom mit 
einer COS-Schwingung (im komplexen Fall auch zusätzlich mit der 
Sinus-Schwingung) und erhält eine Resonanz: Frequenzen im Datenstrom, 
die genau der Filterfrequenz entsprechen, werden vollständig 
durchgelassen."

Könnte ich das für meinen Fall irgendwie anwenden, da ich ja kein 
größeres Passband haben will, sondern im Idealfall nur exakt die eine 
Frequenz?

Grüße,
Third Eye

[H.Joachim S. (crazyhorse)]

Fürs Verständnis: wofür muss man ein Signal, welches man gerade ausgibt, 
wieder messen?

[Tom (Gast)]

H.Joachim S. schrieb:
> Fürs Verständnis: wofür muss man ein Signal, welches man gerade ausgibt,
> wieder messen?

Third E. schrieb:
> Es wird immer eine komplette Periode gespeichert, allerdings ist die
> Amplitude und die Phase beliebig.
     ^^^             ^^^

[Yalu X. (yalu) (Moderator)]

Du multiplizierst das Eingangssignal mit sin und cos und bildest jeweils
das Integral über eine Periode. Daraus kannst da dann direkt die
Amplitude und die Phase des Signals berechnen.

[Der Andere (Gast)]

Third E. schrieb:
> Leider kann es auch Störungen geben, die immer zum selben Phasenwinkel
> auftreten. In dem Fall bringt mir auch die Mittelung nichts.

Wenn du tatsächliche ständig auftretende Signalanteile rausfiltern 
willst. Was ist dann das Ziel?
Das ungestörte Sinussignal? Das kennst du doch wenn du es selbst mit dem 
µC erzeugst?

Ich verstehe den Sinn hier gerade auch nicht.

[Pandur S. (jetztnicht)]

Woher kommen denn die Ausreisser ? EMV, Software Glitches, tatsaechliche 
Stoersignale ?
Ein softwaretiefpass ist eher einfach, laesst sich aber stark beirren. 
Allenfalls waere ein medianfilter besser.

Aber : .. das Ganze soll ?

[Third E. (third-eye)]

Es handelt sich um eine Sensorschaltung. Die Amplitude und Phase möchte 
ich messen.

Die Ausreißer kommen durch Störeinstrahlung (Schaltnetzteile, 
umschaltende µC-Pins, etc.).

[stromtuner (Gast)]

ich brauche das auch: Lautsprecher gibt aus (1kHz), Mikrofon nimmt auf. 
Ich habe ein antialiasing Filter als einziges drinnen.
Ein Sweep 100Hz -10kHz mit mitlaufendem Filter wäre optimal...


StromTuner

[Third E. (third-eye)]

Yalu X. schrieb:
> Du multiplizierst das Eingangssignal mit sin und cos und bildest jeweils
> das Integral über eine Periode. Daraus kannst da dann direkt die
> Amplitude und die Phase des Signals berechnen.

Das hört sich interessant an. Kannst Du das bitte weiter erläutern?
Meine Sinuswerte haben ein Offset von 1/2 der Versorgungsspanung, also 
z.B. der Wert 512. Dieses Offset vorher abziehen?

Wie ergeben sich aus den beiden Ergebnissen dann Amplitude und Phase?

[Dergute W. (derguteweka)]

Moin,

Mitlaufendes Filter ist immer schwierig und laeuft selten so mit, wie 
man's gerne haette. Speziell da und auch bei konstanter Frequenz ist das 
"runtermischen" mit 'nem sin und cos der Frequenz das Mittel der Wahl. 
Danach kann man dann beide Signale jeweils tiefpassfiltern (man kanns 
auch aufsummieren oder integrieren nennen). Beim tiefpassfiltern kann 
man evtl. noch Rechenbumms sparen, wenn man geeignet unterabtastet. Nach 
den Tiefpaessen kriegt man dann 2 Werte, die nur noch so wenig zappeln 
sollten, dass es nicht mehr stoert. Wenn nicht: Dann mehr Tiefpaesse :). 
Aus diesen 2 Werten kann man dann Amplitude + Phase bestimmen.

Gruss
WK

Edit - So ungefaehr:
 Amplitude=sqrt(x*x+y*y);
 Phase=atan2(y,x);
  x=tiefpassgefiltertes(Eingangssignal*cos);
  y=tiefpassgefiltertes(Eingangssignal*sin);

[Third E. (third-eye)]

Dergute W. schrieb:
> Edit - So ungefaehr:
>  Amplitude=sqrt(x*x+y*y);
>  Phase=atan2(y,x);
>   x=tiefpassgefiltertes(Eingangssignal*cos);
>   y=tiefpassgefiltertes(Eingangssignal*sin);

Das probier ich doch glatt mal im Excel aus. Mal sehen, was raus kommt.

[Yalu X. (yalu) (Moderator)]

Third E. schrieb:
> Das hört sich interessant an. Kannst Du das bitte weiter erläutern?

x(t) sei das ungestörte Sinussignal mit der Amplitude a und der
Phasenverschiebung φ, also das, was der ADC lesen würde, wenn es keine
Störungen gäbe:

$$x(t)=a\cdot\cos(2\pi ft+\varphi)$$

s(t) seit das Störsignal, das sich aus Rauschen, Störspitzen und einem
DC-Offset zusammensetzt. y(t) ist das mit der Störung beaufschlagte
x(t):

$$y(t)=x(t)+s(t)$$

Es sollen nun aus y(t) die Amplitude und die Phase bestimmt werden. Dazu
werden die folgenden beiden Integrale über jeweils eine Signalperiode
gebildet:

$$I = \int\limits_0^{1/f}\cos(2\pi ft)\,y(t)\,\mathrm dt$$

$$Q = \int\limits_0^{1/f}\sin(2\pi ft)\,y(t)\,\mathrm dt$$

Dabei ist f und damit auch die Integrationsdauer bekannt. Die
Integration erfolgt einfach durch Aufsummieren des in jedem Abtastzyklus
errechneten Produkts. Aus den Ergebnissen können die Amplitude und die
Phase des ungestörten Signals wir folgt abgeschätzt werden:

$$a'=2f\sqrt{I^2+Q^2}$$

$$\varphi'=-\operatorname{atan2}(Q, I)$$

Die Schätzung x'(t) des ungestörten Signals (Rekonstruktion) sieht also
folgendermaßen aus:

$$x'(t)=a'\cdot\cos(2\pi ft+\varphi')$$

Hier steht auch noch ein Bisschen was dazu:

  https://de.wikipedia.org/wiki/I%26Q-Verfahren

Ich habe das mal an einem Beispiel mit folgenden Parametern durchrechnen 
lassen:

  Anzahl der Abtastungen pro Periode:   80
  Amplitude:                             3
  Phasenverschiebung:                  +21°
  Maximale Rauschamplitude:              0,5
  Störspitze:                           -3
  Offset:                               +1

Die nach obiger Methode geschätzten Signalparameter sind:

  Amplitude (a'):                        2,92
  Phasenverschiebung (φ'):             +22,11°

Man kann die Integration auch über mehrere Perioden laufen lassen, was
die Empfindlichkeit gegenüber stochastischen Störungen reduziert.

Des Weiteren kann man mit dem Verfahren eine erste Schätzung bestimmen,
mit dessen Hilfe die Ausreißer im Signal bestimmen, diese Ausreißer aus
dem Signal entfernen und das Verfahren ein zweites Mal durchlaufen
lassen. Dann wird das Ergebnis noch genauer.

[Third E. (third-eye)]

@Yalu: Vielen Dank! Das nenne ich mal eine ausführliche Erläuterung 
Daumenhoch!
Ein kleiner Schönheitsfehler hat sich aber eingeschlichen. Es müsste 
heißen:
y(t) = x(t) + s(t)

Ich habe inzwischen mit Excel das besagte ausprobiert und es 
funktioniert hervorragend!
Bei dieser Methode fällt auch dankenswerterwerise das Offset der 
abgetasteten Sinuskurve heraus.

Allerdings habe ich zur Berechnung der Amplitude noch eine Frage: Ich 
bilde die beiden Integrale durch "dummes" Aufsummieren der einzelnen 
I-Werte bzw. Q-Werte.
Wenn ich dann berechne
Amplitude = Wurzel(I^2 + Q^2)
ist der berechnete Wert exakt um den Faktor "Anzahl der Abtastungen pro 
Periode" zu groß.

Wie ist da der Zusammenhang zwischen 2*f und "Anzahl der Abtastungen pro 
Periode"?

[Mampf F. (mampf)]

Third E. schrieb:
> Hat jemand eine Idee für einen robusten Algorithmus, wie man diese
> Ausreißer trotzdem los wird?

Probier mal Median statt Mittelwert ...

[Yalu X. (yalu) (Moderator)]

Third E. schrieb:
> Ein kleiner Schönheitsfehler hat sich aber eingeschlichen. Es müsste
> heißen:
> y(t) = x(t) + s(t)

Nicht nur das, die dx in den Integralen sollten dt heißen. Ich hab's
oben geändert, damit nicht noch jemand darüberstolpert.

> ist der berechnete Wert um exakt um den Faktor "Anzahl der Abtastungen
> pro Periode" zu groß.
> Wie ist da der Zusammenhang zwischen 2*f und "Anzahl der Abtastungen pro
> Periode"?

Wenn man das Integral durch eine Summe ersetzt, muss man das dt durch
die Abtastperiode 1/f_a ersetzen, da man ja, um das Integral anzunähern,
die Flächen einzelner Rechtecke mit dieser Breite aufsummiert.

Zieht man den Konstanten Faktor 1/f_a erst vor die Summe und dann auch
noch vor die Wurzel in der Gleichung für a', ergibt sich dort (zusammen
mit dem 2·f) ein Faktor von 2·f/f_a. f/f_a ist aber gerade die Anzahl
der Abtastungen während einer Sinusperiode.