scipy.signal - Filter

[FelixW]

Hallo,

ich möchte gerne einen Tiefpassfilter auf einen diskreten Datensatz anwenden und habe in scipy.signal die passenden Routinen dafür gefunden. Leider sind sie praktisch nicht dokumentiert, was das wirkliche Anwenden etwas schwierig macht...

Dokumentation: http://docs.scipy.org/doc/scipy/reference/signal.html

die besten Beispiele, die ich fand:
http://mpastell.com/2010/01/18/fir-with-scipy/
http://mpastell.com/2009/11/05/iir-filt ... and-scipy/

Die Idee ist glaube ich:
1. Filterkern erstellen:

Code: Alles auswählen

import scipy.signal as signal
kern = signal.firwin(n, cutoff = 0.3, window = "hanning")

2. Filter auf den Datensatz anwenden:

Code: Alles auswählen

filtered = signal.lfilter(kern,1,Datensatz)

Kann mir jemand erklären, wie ich bei einer bestimmten Frequenz filtere? Mir fehlt irgendwie der Zusammenhang zwischen cutoff=x und z.B. cutoff = 250 Hz
x=1 ist übrigens die Nyquist Frequenz = 1 / 2 * Samplingtime (warum das so ist weiß ich leider auch nicht).

Für ein paar Tipps wäre ich überaus dankbar!

Gruß!

[CM]

Hoi,

die Frequenz Deiner Daten kennen numpy/scipy-Routinen natürlich nicht, die kennst Du (hoffentlich ). D. h. wenn Deine Daten mit einer Frequenz von x gesampled wurden, dann ist ein cutoff von 1 die Nyquistfrequenz von x/2. Der Rest ist Dreisatz.

Die Normalisierung auf die Nyquistfrequenz rührt daher, daß ein Lowpassfilter bei Fraktionen dieser Frequenz ein sinnvoller Ansatz ist.

HTH
Christian

[FelixW]

Hallo CM,

Danke für deine Antwort. Tatsächlich ist mir wirklich erst eben klar geworden, wie ich meine cutoff Frequenz bestimme - ganz schön peinlich
Du hast recht, der Filterkern kann meine Sampling Frequenz nicht kennen und die Nyquistfrequenz ist ein sehr sinnvoller Ansatz für einen Filterkern - ich habe mich an der Stelle teilweis missverständlich ausgedrückt (vor dem ersten Kaffee...)

Da ich Forenbeiträge mag, wo auch mal eine Lösung steht, habe ich ein kleines bischen Code geschrieben, dass zeigt, dass der Filter funktioniert und wo "Probleme" sind - sprich, wo ich noch Probleme habe
Ich habe reichlich kommentiert, ich hoffe nicht zu viel für euch Profis:

(Gegeben sei eine Messung eines 10 Hz Sinus, aber unser Messgerät ist schlecht geschirmt und wir haben außerdem 50 Hz Netzbrummen...)

Code: Alles auswählen

# -*- coding: utf-8 -*-

import numpy as N
import pylab as P
import scipy.signal as signal

# Messrate = 1kHz --> Sampling Rate = 1ms
# wir messen 1 Sekunde lang
timescale = N.linspace(0, 1, 1000)

# echtes Signal - 10 Hz Sinus
A = 1.0
Sig = A * N.sin(2 * N.pi * 10 * timescale)

# Stoersignal - die boesen 50 Hz aus der Steckdose
B = 0.1
Noise = B * N.sin(2 * N.pi * 50 * timescale)

# Also messen wir
Messung = Sig + Noise

# Aus der Sampling Rate ergibt sich unsere Nyquist Frequenz
# 1/(2 * 1ms) = 500 Hz
# Aus einer vorhergegangenen FFT Analyse wissen wir, dass wir 50 Hz Stoersignal haben
# der maechtige Dreisatz sagt also 500 Hz : 50 Hz --> 1 : 0.1 = cutoff
# sodass wir unseren Filterkern erstellen koennen 
# (Stuetzstellen und Fenster sind willkuerlich gewaehlt - Verbesserungen gerne!)

n = 31 # Anzahl Stuetzstellen
kern = signal.firwin(n, cutoff = 0.1, window = "hamming")

# Filtern
Smooth = signal.lfilter(kern, 1, Messung)
# macht die ersten n Eintraege unbrauchbar... ggf. Constant / Zero padding der Eingangsdaten?

# Ausgabe zum Angucken...
P.clf()
P.subplot(211)
P.plot(timescale, Sig, timescale, Messung)
P.ylabel('Amplitude [a.u.]')
P.xlabel('Zeit [s]')
P.legend(['echtes Signal', 'Messsignal'])

P.subplot(212)
P.plot(timescale, Smooth)
P.ylabel('Amplitude [a.u.]')
P.xlabel('Zeit [s]')
P.legend(['gefiltertes Signal'])

P.subplots_adjust(hspace=0.5)
P.show()

Verständnissprobleme habe ich folgende:
- wie wähle ich ein passendes n?
- sollte ich meinen Datensatz ähnlich wie bei einer FFT mit Nullen oder Konstanten auffüllen, sodass mir nicht Daten verloren gehen? (Vgl. die ersten n Punkte des gefilterten Blocks)
- woher weiß ich was die zu meinem Problem passende Fensterfunktion ist?

Ich vermute, wenn man sich im k-Raum eher zu Hause fühlt als ich sind die Antworten auf meine Fragen eher trivial, aber für sowas ist das Forum ja da


Edit:
Ich glaube, dass das Problem mit den n Punkten am Anfang dadurch auftritt, dass diese Art von FFT Filterung voraussetzt, dass das Signal a) periodisch ist (haben wir gegeben (; ) und b) dass man es für unendlich lange Zeiten aufnimmt.
Wenn ich mein Numerical Recipies richtig verstehe empfehlen sie großzügiges Zeropadding m Randeffekte zu minimieren - ein Randeffekt ist natürlich das unstetige Ende / der Beginn der Messung. Ich werds nochmal in aller Ruhe durcharbeiten, aber wenn ihr auch Ideen habt, immer her damit!

[CM]

Könnte mich ohrfeigen: Es gab mal eine richtig gute, leider längst veraltete, scipy-Dokumentation, wo Deine Fragen beantwortet waren (glaube ich zumindest). Es war an der Zeit das zu aktualisieren, doch finde ich nur noch dies und das (und vergleichbare Dokumente). Aber irgendwo gab es doch mal einen so schönen Vergleich der Filter ...

Bzgl. N: So klein wie möglich, so groß wie nötig.

Bzgl. Filterwahl: Das kann man so nicht sagen - wir wissen ja noch nicht mal welche Art von Daten Du hast und was Dein Störsignal (oder Rauschen) ist. Wirklich die Steckdose? Ggf. ist das was für Dich. Oder probiere mal a) Gauss, b) Bartlett, c) Wiener. Am besten plottest Du Deine Daten mit versch. Filtern (echte Daten nehmen) und kontrolliere die Ergebnisse. (Trivial, aber effektiv.) Aus dem Bauch heraus würde ich meinen, Hamming ist nicht ideal.

Das ich die Doku jetzt nicht parat habe und selber eingerostet bin fuchst mich. Ich schaue heute Abend nochmal auf meinen Rechner daheim. Vielleicht kann ich dann qualifizierter Antworten.

Gruß,
Christian

[CM]

PS Bin kein Freund von Zero-Padding bei solchen Daten: Wenn Du kannst nimm als Konstante eine Art von Mittelwert links und rechts. Aber das ist u. U. bei periodischen Daten keine gute Wahl, wenn linke und rechte Seite unterschiedlich sein können, weil die Meßperiode nicht vorwählbar ist. Am einfachsten dürfte sein N Punkte mehr zu samplen - wenn das bei Dir möglich ist und nicht zu viel Zeit kostet. Evtl. hilft auch die Samplingfrequenz zu erhöhen - aber da bin ich mir ebenso wenig sicher - weil N dann relativ kleiner wird.

[DaMutz]

ich weiss nicht wie du diese Signale erhälst, aber wenn du einen kontinueirlichen Datenstrom hast, so kannst du beim Filter jeweils die Anfangswerte schon mitgeben, damit hast du das Einschwingen nur beim ersten Block den du verarbeitest. In deinem Beispiel kannst du einfach 2mal Filtern einmal um das Filter zu initialisieren und das 2. mal um das Signal zu filtern:

Code: Alles auswählen

zf = [0]*(len(kern)-1)

# Filter initialisieren
Smooth, zf = signal.lfilter(kern, 1, Messung, zi=zf)
# Filtern
Smooth, zf = signal.lfilter(kern, 1, Messung, zi=zf)

Die Filterordnung ist von deinen Bedürfnissen abhängig. Eine grössere Ordnung erzeugt eine grössere Verzögerung (braucht auch länger zum rechnen), filtert aber besser. Du kannst auch statt eines FIR Filter ein IIR Filter nehmen, dieser Filtertyp filtert bei gleicher Ordnung besser. Hat aber natürlich andere Nachteile.

# edit: nimm als cutoff 0.05 dann wird der 50Hz brumm sehr gut herausgefiltert.

[FelixW]

Hallo,

Danke für die Antworten! Der Link http://docs.scipy.org/doc/scipy/referen ... ignal.html von CM sieht brauchbar aus, ich frage mich wie ich das nicht finden konnte.

Ich habe heute leider nicht so viel Zeit, aber werde versuchen einige von euren Vorschlägen umzusetzen. Ich nehme mir mal vor:
- Constant Padding (ich verweile nicht an einer Position, ich messe einen sich verändernden Abstand, daher ist constant besser als ein Durchschnittswert; ob Zeropadding nicht doch vorzuziehen ist muss ich mir nochmal überlegen bzw. im obigen Beispielprogramm einfach ausprobieren. Ich betrachte die gepaddeten Daten an den Anfang als "Opferlamm" der Filter lässt davon nichts übrig und nach dem Filtern schmeißt man sie wieder raus)
- FIR mit relativ großem Fenster für ein gutes Filterergebnis
- Spielen mit doppeltem Filter anwenden
- vergleiche verschiedene Fensterfunktionen miteinander

Da die Frage aufkam, kann ich euch ja mal kurz umreißen, wie meine Messaufgabe aussieht; wahrscheinlich hiflts euch mir zu helfen

Also ich habe:
- einen Abstandssensor
- eine Doppelblattfeder (als Kraftsensor)

Der Sensor misst die Auslenkung der Feder. Messungen sind im Moment nur Kraft-Abstand Kurven. Glücklicherweise messe ich nur und regle nicht, sodass ich meine Messung gemütlich nachbearbeiten kann und muss keine Rücksicht auf die Prozessgeschwindigkeit nehmen. Die Feder ist makroskopisch groß (anders als beim AFM, aber von der Idee her ähnlich) und schwingt thermisch getrieben. In der FFT erkennt man die Eigenfrequenz der Feder und sie schwingt bei 49 Hz und bei 54 Hz (laterale und vertikale Mode) (ihr könnt euch vorstellen, dass ich verzweifelt nach der Erdschleife gesucht habe ).
Dieses thermische Rauschen überlagert oft die Messung vollständig und somit muss ich etwas dagegen tun. Ich habe zur Auswahl:
- Smoothing (Numerical Recipes: "Smoothing is art, not science" - fällt somit für mich irgendwie raus)
- Filtering (da arbeite ich gerade dran)
- geschickt filtern: Signal FFT, durch die FFT des Rauschspektrums teilen, rücktransformieren - "funktioniert" aber nur bei ausreichend langer Phasenkohärenz des thermischen Rauschens, gleicher Sampling Rate und so weiter.

Ich hoffe, mit den Informationen könnt ihr mir bei der passenden Filter-/Fenster-/Ordnungswahl helfen?
Ich würde mich auch sehr über einen guten Artikel über die verschiedenen Filter FIR/IIR usw. freuen, da ich auf diesem Gebiet eher unerfahren bin (ich hab bis jetzt auf Wikipedia gelesen ).

PS: heute keine Zeit mehr, morgen wahrscheinlich schon!

[FelixW]

Hallo mal wieder.

Ich habe gepaddet, und egal wie man es macht, man gewinnt nichts; es liegt einfach daran, dass man einfach keine "Informationen" hinzufügt. Vielleicht sollte man einfach einige Daten anhängen, die zeitlich zurückliegen (spiegeln), aber ich vermute auch hier lässt sich der Filter bzw. die Theorie dahinter nicht überlisten

Eine weitere Idee von mir war dann, den Filter erst langsam zuzuschalten, indem man im obigen Code einfach diese Funktion aufruft:

Code: Alles auswählen

def filter(Datensatz, Stuetzstellen=61, CutOffFreq=0.08, WindowFunc="hanning"):
# diese rekursive Funktion "schaltet" den Filter langsam ein, bis die gewünschte Filterbreite (Stuetzstellen) erreicht ist.
# Also wird der erste Teil rekursiv verkleinert und mit weniger Stuetzstellen bearbeitet.

# dringend noch zu implementieren: Phasenversatz ausgleichen

	if len(Datensatz) == 1:
		return Datensatz
	if len(Datensatz) < Stuetzstellen:
		print "Denk nochmal drueber nach..."
		return Datensatz
	sub1 = Datensatz[:Stuetzstellen]
	
	kern = signal.firwin(Stuetzstellen, cutoff = CutOffFreq, window = WindowFunc)
	Smoothed = signal.lfilter(kern, 1, Datensatz)
	# rekursives Ersetzen der "verlorenen" Daten
	Smoothed[:Stuetzstellen] = filter( sub1, N.ceil(Stuetzstellen/2.0), CutOffFreq, WindowFunc) 
	
	return Smoothed

Leider haben Filter die unpraktische Eigenschaft, einen Phasenversatz zu machen (vgl. Bode-Diagramm); dieser macht dieses Zusammenbasteln nicht besser als die n Stellen direkt wegzulassen.
Ich weiß zum einen nicht so recht, wie sich der Phasenversatz in Abhängigkeit von der Anzahl der Stützstellen, Nyquistfrequenz und Fensterfunktion abhängen; zum anderen weiß ich nicht, wie man den Phasenversatz sinnvoll korrigieren sollte.
Daher werde ich nun erstmal einfach nach Gefühl filtern und auf diesem Wege die Stützstellenzahl sowie meine Fensterfunktion bestimmen und meine Datenerfassung entsprechend anpassen.

[DaMutz]

ich weiss jetzt nicht was dein Problem ist. Jedes Filter erzeugt eine Verzögerung (je grösser die Ordnung um so grösser die Verzögerung), auch ein analoges Tiefpassfilter (z.B. RC-Glied) erzeugt eine Verzögerung. Diese Verzögerung ist bei einem Filter nicht weg zu bringen.
Oder willst du für die Darstellung das Eingangssignal synchron mit dem Ausgangssignal des Filters haben. In diesem Fall könntest du auf den grössten Wert oder die Nulldurchgänge synchronisieren. Aber was dir das genau bringt weiss ich nicht.

In der Praxis / Messung kannst du doch nicht mit so schönen Sinus Schwingungen rechnen.

über die Window Funktion würde ich mir nicht all zu grosse Sorgen machen, lass diesen Wert doch einfach weg.

[FelixW]

Ja DaMutz, du hast recht. Die Phase interessiert bestenfalls bei Schwingungen und wenn sie bekannt sind könnte man sie nachträglich synchronisieren. Ich war bis vor Kurzem was Filter angeht leider sehr unvorbelastet und erarbeite mir das Thema gerade erst und einige Stellen, die sich vermutlich aus der Fouriertheorie ergeben springen mir leider (noch) nicht ins Auge
Ich habe in meinen Messungen auch nichts periodisches und habe einfach mal einen gewöhnlichen Filter drüberrennen lassen und war mit dem Ergebnis sehr zufrieden. Die ersten N Stellen schmeiße ich hinterher einfach raus, weil sie etwas sch... im plot aussehen und keine Informationen mehr tragen.
Einen kleinen "Trick" habe ich noch; nämlich bei vielen Messungen wird der selbe Datenstrom von der Messsoftware in z.B. vorwärts und rückwärts aufgespalten; die kann man für den Filtervorgang einfach aneinander hängen (zeitliche Reihenfolge beachten) und verliert somit nur einmal N Datenpunkte

Danke nochmal für eure Hilfe und ein schönes WE!

[DaMutz]

die Messpunkte aneinander zu hängen ist eine Möglichkeit. Eleganter ist es aber den Filterstate von der vorherigen Messungen weiter zu verwenden. Dies geschieht indem du den Filter folgendermassen aufrufst:

Code: Alles auswählen

# Filtern
Smooth, zf = signal.lfilter(kern, 1, Messung, zi=zf)

zf muss eine Liste mit der Länge der Ordnung des Filters sein.

Schau dir das Beispiel von mir weiter oben an. Dann sind meine Erklärungsversuche vielleicht auch klarer.