Verfasst am: 12.06.2015, 15:54
Titel: Signal Filterung
Hallo allerseits,
es geht um den Entwurf eines Filters um das Ursprungssignal aus einem gestörten zu extrahieren. Dabei wurde das Signal durch 3 Störquellen gestört.
Die Abtastfrequenz des Signals ist 200Hz.
Anbei sind das gestörte Signal und das originale skizziert.
Das gestörte signal tue ich in einer Zip Datei mit
Welche FIR und IIR Filtern kommen Euer Meinungen nach in Frage?
Ich habe noch nie mit Filtern gearbeitet und bräuchte daher ein Hinweis um in die richtige Richtung zu gehen
So ohne weitere Infos kann man wohl keinen vernünftigen Filter vorschlagen, da aus den Bildern nicht ersichtlich wird, welchen Frequenzbereich Nutz- und Störsignal haben. Man führt vor dem Filterentwurf daher immer eine Analyse des Frequenzspektrums vor. Da du offensichtlich das Nutzsignal auch ungestört vorliegen hast, macht es die Sache etwas einfacher.
Du brauchst zunächst folgende Funktion zur Erstellung eines Betragsspektrums. Die Phase des Signal interessiert uns eigentlich nicht.
Code:
function[mag, mag_dB, fv] = FFT_betragsspektrum( signal, nfft, fa, scale) ;
% Input: % Signal im Zeitbereich % nfft = Anzahl Messwerte für fft % wenn nfft > length(sig) -> fft(sig,nfft) führt Zeropadding durch % fa = Abtastfreq. % scale: 0 = keine Impulsantwort als Eingang, 1 = Impulsantwort % Output: % Magnitude des Spektrums linear und dB skaliert % Frequenzvektor fv in [Hz] von 0...fa/2
% un-,gerade Anzahl Messwerte? ifmod(nfft,2) == 0;
k = (nfft/2) + 1;
else
nfft = nfft + 1;
k = (nfft/2) + 1;
end
Vielen Dank für die ausführliche und kompetente Antwort !
Ich habe mich vielleicht nicht richtig ausgedruckt!
Ich habe nur das gestörte Signal(Bitte Anhang sehen) . Die Aufgabe ist das orginale Signal rauszufiltern! Das was wir wissen ist nur, dass das ein EKG Signal ist, das durch 3 Störquellen gestört wurde und mit einer Abtastrate von 200 Hz abgetastet wurde. Es soll dabei die Detektion von Herzschlägen und die genaue Form des Herzschlagsignals ermöglicht werden. Die Filterkette sollte also die Form des Herzschlages möglichst originalgetreu beibehalten.
Das gestörte Signal habe ich nach Ihren Anweisungen im Frequenzbereich dargestellt, anbei ist das Plot dafür. Für nfft habe ich length(signal)*10 also 4*10^5 genommen.
da nur das gestörte Signal vorhanden ist, kann ich momentan nur raten. Ich habe bisher derartigte Signale nicht untersucht und bin mir daher nicht sicher was in dem Spektrum zum Nutzsignal gehört.
Ich tippe mal darauf das der Peak im Spektrum bei ca. 50Hz zu der Halbwelle gehört oder aber die Frequenz von dem Peak+Halbwelle (beziehe mich auf Bild original.png) ist. Daher darf dieser Anteil nicht gefiltert werden. Die Peaks im EKG Signal haben mehrere Frequenzanteile, da die Fouriertransformation (wird zur Erstellung des Frequenzspektrums genutzt -> siehe Funktion FFT) solche Signale nur aus sinus und cosinus Wellen zusammensetzen kann. Im Spektrum ist hier vor allem der niederfrequente Anteil bis ca. 20 Hz entscheidend.
Ich würde zunächst einmal einen Tiefpass mit einer Grenzfrequenz > 50Hz einsetzen um die Störungen oberhalb von 50Hz zu entfernen.
Code:
Fa = 200; % Hz
Fn = Fa/2;
Fg = 55; Hz Grenzfrequenz des Tiefpass
[B,A] = butter(4,Fg/Fn,'low'); % IIR Tiefpassfilter
EKGfiltered_1 = filter(B,A,EKG);
% Alternativ filtfilt = doppelte Filterung und Aufhebung des Filterdelays
EKGfiltered_1 = filtfilt(B,A,EKG);
Möglicherweise muss aber auch der Bereich ab ca. 20Hz bis zum 50Hz Peak gefiltert werden. Dieser Bereich kann mit einem Bandstop Filter gefiltert werden.
Code:
Fg1 = 20Hz; % untere Grenze
Fg2 = 45Hz; % obere Grenze [B,A] = butter(4,[Fg1/Fn, Fg2/Fn],'stop'); % IIR Bandsperre
EKGfiltered_2 = filter(B,A,EKGfiltered_1);
% Alternativ filtfilt = doppelte Filterung und Aufhebung des Filterdelays
EKGfiltered_2 = filtfilt(B,A,EKGfiltered_1);
Ich habe den Code aber nicht getestet und auch nicht überprüft welche Dämpfung und Dämpfungsverlauf die beiden Filter erzeugen noch deren Stabilität. Sonst mal die Filterordnung von derzeit 4 ändern. Ob es Filter mit geeigneteren Eigenschaften (z.B. IIR Bessel, Chebyshev oder ein FIR Window Sinc) für dein Signal gibt, würde ich dann im zweiten Teil testen. Erstmal schauen ob das Prinzip Tiefpass + Bandsperre überhaupt sinnvoll ist.
Schau dir erstmal das gefilterte Signal im Zeit- und Frequenzbereich an. Nach wie vor wäre es hilfreich mal ein Frequenzspektrum eines ungestörten Signals zu sehen wie in original.png abgebildet. Es müss ja nicht exakt passen...nur damit man mal der Verlauf im Spektrum sieht und somit sagen kann, was man nicht filtern darf.
Allerdings bitte nochmal das Frequenzspektrum des gestörten Signals plotten mit einer Fensterfunktion um den Leakage Effekt zu vermeiden. Nicht das wir auf Signalanteile reinfallen, die eigentlich gar nicht vorhanden sind.
Code:
signal_orig = ...% Nutzsignal
nfft = 2^nextpow2(length(signal_orig)); % nächste 2er Potenz
fa = 200; % Abtastfrequenz in Hz des Signals
nfft=2^nextpow2( length(signal_orig) ); ist hier keine gute Idee, da für die Fensterung mit Amplitudenkorrektur, das Signal damit multipliziert wird und deswegen beide die gleiche lenge haben sollen. Ich habe also nfft=length(signal) gesetzt.
Die plots für die 3 Signale ( gestörtes Signal, gefiltert1 mit IIR Tiefpassfilter und gefiltert2 mit IIR Bandsperre ) sind anbei.
Wäre das eigentlich nicht sinnvoll ein Herzschlag zu bestimmen, da ein EKG eine Folge von mehreren ist ? Und dann versuchen dieses zu filtern?
Bei der Filterung hast du vermutlich etwas missverstanden. Es wird NICHT das Signal signal_win, also Signal mit Fensterung verwendet. Das Fenster wird nur für die Darstellung des Frequenzspektrums verwendet. Für die Filterung nutzt zu weiterhin nur dein Rohsignal.
Die Plots sind aber zum Vergleich so nicht zu gebrauchen. Nutze die Funktion
um das Signal auf eine oder wenige Perioden des EKG/Herzschlags einzugrenzen. So erkennt man jedenfalls nicht viel. Außerdem wäre es sinnvoll die unterschiedlichen Signale in einem Plot darzustellen.
Code:
N = 1:length(signal);
plot(N,signal,[0.6510.6510.651], N,EKGfiltered_1,'b', N,EKGfiltered_2,'r') legend('Rohdaten','Tiefpassfilter','Bandpassfilter') xlabel('Sample N') ylabel('Amplitude') xlim([1 ???])% nur Ausschnitt darstellen (wenige Perioden) grid on;
Für die Darstellung des Frequenzspektrums wäre es ebenfalls sinnvoll, Rohdaten, EKGfiltered_1 + 2 in einem Plot darzustellen und am Besten auch die selben Farben zu verwenden. Im Frequenzbereich wird aber nicht mit xlim gearbeitet .
Zitat:
Wäre das eigentlich nicht sinnvoll ein Herzschlag zu bestimmen, da ein EKG eine Folge von mehreren ist ? Und dann versuchen dieses zu filtern?
Ob nun das ganze Signal oder nur ein Ausschnitt gefilter wird, ist erstmal egal. I.d.R ist es besser etwas mehr Messwerte zu haben, damit sich das Filter einschwingen kann. Für die Darstellung der Daten ist es wie oben erwähnt schon sinnvoll nur einen Ausschnitt zu zeigen.
Als weiteren Schritt wäre es sinnvoll dann mal die Ordnung des Filters zu erhöhen um eine höhere Dämpfung zu erhalten.
Prinzipiell glaube ich aber nicht, dass es zu deinem Thema nicht schon Lektüre gibt. Du bist mit Sicherheit nicht der Erste, der EKG Signale filtern muss. Man muss das Rad ja nicht neu erfinden...sich bei anderen Input zu holen ist nie verkehrt. Ich würde also mal im Inet auf Suche zum Thema EKG und Filterung gehen. Es muss ja kein Matlab Code sein...Hinweise welche Filtertypen sich dafür besonders eignen, wäre schon eine Hilfe.
hier sind die plots!
Was kann man eigentlich der Frequanzanalyse entnehmen?
Und woran erkennt man wie weit man noch filtern soll bzw. darf umd das gewollte Signal zu haben ?
ich habe mal selbst ein wenig mit deinem Signal getestet. In Matlab gibt es eine Funktion ecg zur Erzeugung eines idealierstierten EGK Signals. Somit haben wir schon mal einen Vergleich, wie denn in etwa das Nutzsignal im Zeit und Frequenzbereich aussehen muss.
Mit den Filtern haben ich lediglich kurz getestet, muss mir das aber noch genauer ansehen. Einfach nur den m-file ausführen. Das Notch Filter wie auch dein Signal werden dann geladen, so fern sie im selben Verzeichnis wie EGC_filtering.m liegen. Zumindest sieht man nun schon mal halbwegs das EGK Signal und nichtz nur Rauschen.
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
.
RSS