Implementierung eines Low pass Filters

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

Das kann man auswählen. Und was np.vectorize macht einfach nachlesen in der Dokumentation.
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Achso okay danke! Aber wie kann man gleichzeitig mit 1 Hz und 5Hz filtern, bzw. aus der Kombination aus den Parametern?
Ich dachte du hättest mit 1 Hz gefiltert.
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

Kann man nicht. Wenn man in Python an den selben namen etwas Neues bindet, ist das alte Objekt darüber nicht mehr erreichbar. Und in diesem Fall einfach weg.
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Gut zu wissen vielen Dank!

Ich beschäftige mich gerade mit der Segmentierung. Ich mache es über das Verfahren der Peak detection.

Code: Alles auswählen

filtered_data = filter(sensor_data)
Die gefilterten Daten sind ja in einem mehrdimensionalem Array. Ich habe versucht das zu vereinfachen in dem ich das Array mit

Code: Alles auswählen

filtered_data2=filtered_data.ravel()
zu einem eindimensionalem Array zusammenfüge.

Denn was ich so an Peak detection Verfahren finde, sind sie meist mit einem 1d Array.

Problem ist nur das keine Peaks gefunden werden.

Code: Alles auswählen

filtered_data = filter(sensor_data)
filtered_data2=filtered_data.ravel()


peaks = find_peaks(numbers, height = 300, threshold = 1, distance = 5)

peaks = find_peaks(filtered_data2, height = 1, threshold = None, distance=200)
height = peaks[1]['peak_heights'] #list of heights of peaks
peak_pos = peaks[0]
print(peaks)

# plot the peaks
fig = plt.figure()
ax = fig.subplots()
ax.plot(filtered_data2)
ax.scatter(peak_pos, height,color = 'r', s = 25, label = 'Maxima')
ax.legend
plt.show()
Da kommt ein leeres Array raus für die Peaks.

https://ibb.co/3zYVzfy

Das ist nur ein Beispiel von 3 Peak detection die wo nichts angezeigt wurde und ich rumprobiert habe. Ich weiß jetzt nur nicht, ob ich anders herangehen muss strukturell, da die Daten in einem mehrdimensionalen Array sind oder ich einen falschen Gedankengang habe?
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Da ich verständlicherweise meinen Beitrag nicht bearbeiten kann, stelle ich den kompletten Code hier rein und kann konkretesieren, was mich beschäftigt.

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt
 
# BW 20Hz cutoff at 100Hz samplerate
# 2.376⋅yi = (1⋅xi + 1⋅xi-1) + 0.376⋅yi-1

class BW:
    
    def __init__(self, a, b):
        self._a = a
        self._b = b
        self._xn_1 = None
        self._yn_1 = 0
        
    def __call__(self, x):
        if self._xn_1 is None:
            self._xn_1 = x
            
        y = self._yn_1 = ((x + self._xn_1) + self._a * self._yn_1) / self._b
        self._xn_1 = x
        return y
            
# 5Hz
filter = BW(5.376, 7.376)
# 1Hz
filter = BW(30.821, 32.821)
filter = np.vectorize(filter)

data=pd.read_csv('samplesdata.csv',sep=";", decimal=",",encoding='latin-1')
sensor_data=data[['Euklidische Norm']]
sensor_data=np.array(sensor_data)
      

filtered_data = filter(sensor_data)
Das ist der Butterworth filter und ich möchte nun aufgrund dessen die Daten durch Peak Detection (Hochpunkte) die Daten segmentieren.

Ich habe versucht die Daten mit dieser Zeile

Code: Alles auswählen

filtered_data2=filtered_data.ravel()
zu vereinfachen, denn dadurch sind die Zahlen in einem Array dargestellt, denn sonst stellt jede Zahl ein eigenes Array dar oder?
Denn so was ich gesehen habe sind viele Peak Detection mit einem 1D Array.

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt
 
# BW 20Hz cutoff at 100Hz samplerate
# 2.376⋅yi = (1⋅xi + 1⋅xi-1) + 0.376⋅yi-1

class BW:
    
    def __init__(self, a, b):
        self._a = a
        self._b = b
        self._xn_1 = None
        self._yn_1 = 0
        
    def __call__(self, x):
        if self._xn_1 is None:
            self._xn_1 = x
            
        y = self._yn_1 = ((x + self._xn_1) + self._a * self._yn_1) / self._b
        self._xn_1 = x
        return y
            
# 5Hz
filter = BW(5.376, 7.376)
# 1Hz
filter = BW(30.821, 32.821)
filter = np.vectorize(filter)

data=pd.read_csv('samplesdata.csv',sep=";", decimal=",",encoding='latin-1')
sensor_data=data[['Euklidische Norm']]
sensor_data=np.array(sensor_data)
      

filtered_data = filter(sensor_data)
filtered_data2=filtered_data.ravel()


peaks = find_peaks(numbers, height = 300, threshold = 1, distance = 5)

peaks = find_peaks(filtered_data2, height = 1, threshold = None, distance=200)
height = peaks[1]['peak_heights'] #list of heights of peaks
peak_pos = peaks[0]
print(peaks)

# plot the peaks
fig = plt.figure()
ax = fig.subplots()
ax.plot(filtered_data2)
ax.scatter(peak_pos, height,color = 'r', s = 25, label = 'Maxima')
ax.legend
plt.show()
Das ist ein Beispiel was ich probiert habe an Peak detection algorithmen und probiert habe.
Jedoch findet es keine Peaks
https://ibb.co/3zYVzfy

Liegt es daran das ich von vornerein es falsch mache in dem ich die gefilterten Daten ändere oder was für eine Struktur haben sie, dass es mit Peak detection-Verfahren nicht funktioniert?
Benutzeravatar
Dennis89
User
Beiträge: 1123
Registriert: Freitag 11. Dezember 2020, 15:13

Hallo,

aber der Code läuft doch so gar nicht. Es gibt kein 'pd', kein 'find_peaks', und keine 'numbers'.
Es wurde ja schon gesagt, dass das was an Namen gebunden wird überschrieben wird, wenn man an den Namen wieder was bindet.

Irgendwie kann ich dir nicht recht folgen.

Grüße
Dennis
"When I got the music, I got a place to go" [Rancid, 1993]
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Stimmt du hast recht, ich habe was ausprobiert gehabt und es nicht zurück geändert und die Importe wurden davor schon gemacht, daher war es nicht dabei.

Komischerweise funktioniert es jetzt :shock:

Code: Alles auswählen

import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
import numpy as np

 
# BW 20Hz cutoff at 100Hz samplerate
# 2.376⋅yi = (1⋅xi + 1⋅xi-1) + 0.376⋅yi-1

class BW:
    
    def __init__(self, a, b):
        self._a = a
        self._b = b
        self._xn_1 = None
        self._yn_1 = 0
        
    def __call__(self, x):
        if self._xn_1 is None:
            self._xn_1 = x
            
        y = self._yn_1 = ((x + self._xn_1) + self._a * self._yn_1) / self._b
        self._xn_1 = x
        return y
            
# 5Hz
filter = BW(5.376, 7.376)
# 1Hz
filter = BW(30.821, 32.821)
filter = np.vectorize(filter)

data=pd.read_csv('samplesdata.csv',sep=";", decimal=",",encoding='latin-1')
sensor_data=data[['Euklidische Norm']]
sensor_data=np.array(sensor_data)
      

filtered_data = filter(sensor_data)
filtered_data2=filtered_data.ravel()


peaks = find_peaks(filtered_data2, height = 300, threshold = 1, distance = 5)

peaks = find_peaks(filtered_data2, height = 1, threshold = None, distance=200)
height = peaks[1]['peak_heights'] #list of heights of peaks
peak_pos = peaks[0]
print(peaks)

# plot the peaks
fig = plt.figure()
ax = fig.subplots()
ax.plot(filtered_data2)
ax.scatter(peak_pos, height,color = 'r', s = 25, label = 'Maxima')
ax.legend
plt.show()
Also für diejenigen die ähnlich wie ich einen Peak detection suchen für einen 1d Array, das ist der Code der funktioniert zumindest bei meinen Daten.
Anscheinend hat es geholfen die gefilterten Daten in einem eindimensionalen Array zusammenzufügen.

Ich weiß nur tatsächlich nicht wofür die Parameter distance und Height stehen. Bei distance habe ich gedacht das man da die erwartbare Distanz zwischen zwei Peaks angeben soll, aber ich bin mir nicht sicher.

Weiß da vielleicht jemand mehr, was man bei den beiden Parameter angibt?
Benutzeravatar
Dennis89
User
Beiträge: 1123
Registriert: Freitag 11. Dezember 2020, 15:13

Klaro und derjenige hat das zum Glück aufgeschrieben:
https://docs.scipy.org/doc/scipy-1.1.0/ ... peaks.html

Grüße
Dennis
"When I got the music, I got a place to go" [Rancid, 1993]
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Ich hätte eine Frage bzgl. des Low-Pass-filters.

Code: Alles auswählen

class BW:
    
    def __init__(self, a, b):
        self._a = a
        self._b = b
        self._xn_1 = None
        self._yn_1 = 0
        
    def __call__(self, x):
        if self._xn_1 is None:
            self._xn_1 = x
            self._yn_1 = x
        y = self._yn_1 = ((x + self._xn_1) + self._a * self._yn_1) / self._b
        self._xn_1 = x
        return y
            
# 5Hz
filter = BW(5.376, 7.376)
# 1Hz
filter = BW(30.821, 32.821)
Dieser Code ist ja angelehnt auf die Parameter mit dem Order=1 und dem Cutoff= 1 Hz.
Was mir aufgefallen ist das bei wenn man den Order erhöht auf 2, dann hat man eine komplett andere Formal mit mehr variablen
Wie müsste man den Code umändern, wenn man einen höheren Order haben will?

Das ist mir aufgefallen, als ich rumspielen wollte mit den Werten, mit dem Link: https://www.meme.net.au/butterworth.html
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

Wenn du dir klar machst, was die Werte sind, und wie das mit der Ordnung zusammenhaengt, dann ist es trivial, den Code entsprechend abzuaendern.
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Ich habe es mit Order=2 probiert.

Code: Alles auswählen

class BW:
    
    def __init__(self, a, b,c):
        self._a = a
        self._b = b
        self._c = c
        self._xn_1 = None
        self._yn_1 = 0
        self._yn_2 = 0

        
    def __call__(self, x):
        if self._xn_1 is None:
            self._xn_1 = x
            self._yn_1 = x
            self._yn_2 = x

        y = self._yn_1 = ((x + 2*x + self._xn_1) + self._a * self._yn_1+ self._c*self._yn_2) / self._b
        self._xn_1 = x
        return y
            
# 1Hz
filter = BW(-0.948,1.789, 4.841)
Ich denke ich muss noch was am Ende beim dividieren hinzufügen. Jedoch erschließt sich mir nicht durch was dividiert werden soll.
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

Das ist falsch, weil du nicht nur fuer y zwei Werte brauchst. Und warum du da was dividieren willst, verstehe ich nicht.

Gibt es einen echten Grund, warum du auf einen BW von 2ter Ordnung gehst? Oder ist das nur wieder so eine "das macht man halt so"-Nummer vom werten Herr Professor?
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

Oh, und da ist noch mehr falsch, weil du die y n nicht korrekt updatest. Was sind die y_n?
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Ich wollte in zweiter Ordnung gehen und schauen wie der Graph aussieht, denn ich wurde tatsächlich gefragt, ob ich es in einer anderen Ordnung probiert habe^^

Am Ende hast du ja durch self.b dividiert, daher meine Frage.

y_n stehen für die y Variablen in der Formel, denke ich.
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

Ich habe dadurch dividiert. Wenn du dir die Formel auf der Seite anschaust, dann siehst du ja, warum. Denn die sieht anders aus. Und was ist dann da ein Unterschied zur 2ten Ordnung?

Was die y_n angeht - ja gut. Das war jetzt die simple Antwort. Denn offensichtlich hat das was mit der Formel zu tun. Viel mehr gibt's ja nicht. Aber diese y_n (genauso wie die x_n) stehen fuer etwas ganz spezifisches, wenn man das als diskretes Signal betrachtet. Und das musst du dir klar machen, und dann sollte hoffentlich auch klar sein, wie yn_2 und xn_2 zustande kommen. Und das das, was du da gemacht hast, nicht reicht.

Wenn du keine Ahnung von digitalen Filtern hast (weiss nicht, ob du solltest, oder ob deinem Prof reicht, etwas zu benutzen, dass "einfach funktioniert"), dann solltest du ggf. eher investieren, da eine "offizielle" Implementierung wie zB aus scipy zum laufen zu bekommen. Oder vielleicht eine andere, das wird denke ich nicht die einzige sein. Du weisst jetzt ja, wie es aussehen soll. Denn mir scheint, dass deine Programmierkenntnisse eher wenig ausgepraegt sind - eher auf der Ebene, gegebenes zu verstoepseln. Darum wird eine allgemeine Implementierung zu erstellen dich wahrscheinlich ueberfordern.

Wobei es IMHO so ganz ohne Verstaendnis deren Wirkungsweise auch nicht geht, denn woher sonst soll eine schluessige Beurteilung der Ergebnisse kommen. "Sieht anders aus" ist ja nicht so richtig gut als Antwort auf die Frage, warum man denn einen 2ter Ordnung gewaehlt hat.
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

So, ich habe mir nicht vorstellen koennen, dass die scipy BW implementierung wirklich nicht geht. Wie gesagt, letzte Woche hatte ich keinen Zugriff darauf - weil nur iPad. Also habe ich mir gerade mal angeschaut, was denn da mit deinem Code eigentlich so vor sich geht.

Und im Grunde hattest du die Loesung schon selbst halb hergestellt - das Problem ist die Array-Repraesentation deiner Daten. Statt einem 1-dim Array (was die sind), ist das ein daemliches 2-dimensionales Array, mit einem shape (N, 1). Wieso auch immer.

Mit deiner selbst genutzten ravel funktioniert dann der Filter ganz wunderbar. Und wenn man dann noch den DC-Offset wegnimmt, indem man den Durchschnitt der Daten berechnet, und den abzieht, eleminiert man auch das Anfangswert-Problem. Ich habe keinen Weg gefunden, wie man den Filter anders irgendwie "vorwaermt", wobei man natuerlich auch einfach kuenstlich Daten bis zum Einschwingen erzeugen kann.

Damit sollte dir das ganze deutlich leichter fallen, und zum Vergleich der verschiedenen Filter-Ordnungen und Cutoffs bietet sich eine FFT-Darstellung an. Denn sonst ist eben nur schoene Kurven diskutieren.

Code: Alles auswählen

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
import pandas as pd


def butter_lowpass(cutoff, fs, order=5):
    nyq = fs/2
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a


def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y


data=pd.read_csv('samplesdata.csv',sep=";", decimal=",",encoding='latin-1')
sensor_data=data[['Euklidische Norm']]
sensor_data = np.array(sensor_data).ravel()
sensor_data = sensor_data - np.average(sensor_data) # eliminate DC

# Filter requirements.
order = 10
fs = 100  # sample rate, Hz
cutoff = 1

y = butter_lowpass_filter(sensor_data, cutoff, fs, order)

plt.figure()
plt.subplot(211)
plt.plot(sensor_data, color='green')
plt.grid(color='green', linestyle='--', linewidth=0.5)
plt.title('original signal')

plt.subplot(212)
plt.plot(y, color='red')
plt.title('Filtered')
plt.grid(color='green', linestyle='--', linewidth=0.5)
plt.show(block=True)
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Du hast recht, das Filtern dient zur Datenvorverarbeitung und wenn ich damit segmentieren kann sollte es reichen. Die Daten sollten nicht zu sehr alles wegfiltern, dann sind die Daten zu sehr verzerrt.

Danke übrigens, für die Lösung mit Scipy! Das hilft mir sehr fürs Verständnis!

Ich habe es mit meinem Peakdetection ausprobiert und irgendwie, findet es garkeine Peaks.

Code: Alles auswählen

[
import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
import pandas as pd


def butter_lowpass(cutoff, fs, order=5):
    nyq = fs/2
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a


def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y


data=pd.read_csv('samplesdata.csv',sep=";", decimal=",",encoding='latin-1')
sensor_data=data[['Euklidische Norm']]
sensor_data = np.array(sensor_data).ravel()
sensor_data = sensor_data - np.average(sensor_data) # eliminate DC

# Filter requirements.
order = 10
fs = 100  # sample rate, Hz
cutoff = 1

y = butter_lowpass_filter(sensor_data, cutoff, fs, order)

peaks = find_peaks(y, height = 1, threshold = None, distance=200)
height = peaks[1]['peak_heights'] #list of heights of peaks
peak_pos = peaks[0]
print(peaks)

# plot the peaks
fig = plt.figure()
ax = fig.subplots()
ax.plot(y)
ax.scatter(peak_pos, height,color = 'r', s = 25, label = 'Maxima')
ax.legend
plt.show()
Eigenlich müsste es doch die selben Peaks sogar finden bei einer Ordnung von von 1. Aber es kommt allgemein nur ein leeres Array raus.

https://ibb.co/WfwgZn8

Also das Array ist sogar 1d, aber irgendwie hat es mit der Scipy implementierung ein Problem.
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Okay , ich habe das Problem gelöst. Ich muss den Height Paramater natürlich anpassen, denn jetzt hat sich ja die Skalierung geändert. Ich habe den Height auf 0,00 gestellt.
Den Order habe ich jetzt auf 3 gestellt, denn sonst wären nicht alle Peaks korrekt erkannt worden.

Code: Alles auswählen

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
import pandas as pd


def butter_lowpass(cutoff, fs, order=5):
    nyq = fs/2
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a


def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y


data=pd.read_csv('salim_acc3.csv',sep=";", decimal=",",encoding='latin-1')
sensor_data=data[['Euklidische Norm']]
sensor_data = np.array(sensor_data).ravel()
sensor_data = sensor_data - np.average(sensor_data) # eliminate DC

# Filter requirements.
order = 3
fs = 100  # sample rate, Hz
cutoff = 1

y = butter_lowpass_filter(sensor_data, cutoff, fs, order)

peaks = find_peaks(y, height = 0.00, threshold = None, distance=200)
height = peaks[1]['peak_heights'] #list of heights of peaks
peak_pos = peaks[0]
print(peaks)

# plot the peaks
fig = plt.figure()
ax = fig.subplots()
ax.plot(y)
ax.scatter(peak_pos, height,color = 'r', s = 25, label = 'Maxima')
ax.legend
plt.show()
Leider nimmt es einen Peak zu viel, egal was ich einstelle ist der 11. Peak nicht richtig gesetzt und der letzte wird nicht erkannt. Das ist das Beste Ergebnis was ich erreichen konnte bei diesen Daten. Ich habe noch mehr Daten aufgenommen, da werden alle Peaks richtig gefunden bis oft auf den letzten.

https://ibb.co/xLTZsXx

So sehen die meisten Daten aus, da wird der letzte Peak nicht erkannt. Ich werde wahrscheinlich den letzen Peak rausschneiden müssen immer.
Oder hat jemand einen besseren Vorschlag?
https://ibb.co/94k2LLZ


Der nächste Schritt zum segmentieren wäre ja die Werte zwischen den Peaks zu bekommen oder?
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

Ich würde ja eher die negativen peaks nehmen, die sehen viel besser aus. Und für die Segmentierung ist das egal.
Risingsun
User
Beiträge: 43
Registriert: Dienstag 8. November 2022, 14:56

Da eine Bewegung wie ein "M" aussieht, muss ich dann jeden zweiten negativen Peak nehmen oder?
Antworten