www.mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Notch-Filter


Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Moin, ich benötige einen digitalen 50 bzw 60Hz notchfilter.
Den hat doch garantiert schon jemand programmiert. Wer kann ihn mir
online zur Verfügung stellen?

Danke, Thomas

Autor: Alexander (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
http://www.nauticom.net/www/jdtaft/Notch.htm

Entwirf ihn dir selbst...

Hast du die Abtastfrequenz nicht genannt weil du keine Ahnung hast oder
weil du es vergessen hattest? Digitale Filter sind nämlich schwerlich
portabel. Programmierst du in Assembler, C, C++, Basic, ...?

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
habs vergessen ... 200 Hz !!

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
achso ..und ich programmiere in c !

die abtastfrequenz wird doch eh als konstante angegeben und cih dachte
an einen algor. der sich daran anpasst ...
ich habe noch nie einen notch filter entworfen.
gibts dazu nen literatur tip ?

Thomas

Autor: Alex (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Die Abtastzeit als Konstante - interessanter Ansatz, wenn man den
Filterentwurf während der Laufzeit auf dem Controller machen möchte.
Warum verwendest du nicht das verlinkte Applet, um an die nötigen
Koeffizienten zu kommen, kannst alternativ auch Matlab oder sonstige
Filtertools dafür nehmen. Die Differenzengleichung musst du dann
einfach implementieren. 200Hz wird aber meiner Meinung nach für eine
ordentliche Arbeitsweise unter Umständen nicht ausreichen, aber ich
kenne ja deine Anforderungen nicht.

Helfen kann dir jedes gängige DSP-Buch, was der Filter macht ist für
die Implementierung eigentlich egal, die Implementierung ist für alle
FIR bzw. IIR Filter gleich. Man braucht halt einen bzw. zwei
Ringspeicher und die entsprechenden Arrays mit den Koeffizienten. Dann
läuft einfach in einer Schleife ein multiply&accumulate ab und das
wars.

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
schon klar, ich will aber nen IIR, der link lässt mich nur nen FIR
berechnen ob hab ich was übersehen ?

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
.... dann brauch ich auch nur nen filter 2ter ordnung ..
übrigens ist bei der filterfrequenz von 50 bzw 60 hz das abtasttheorem
eingehalten. 200Hz sind ausreicehnd ! zummindest theoretisch !

Autor: Detlef A (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
haste doch schon mal gefragt, habe ich auch schon mal drauf geantwortet,
oder gibts mehr als einen Thomas?

http://www.mikrocontroller.net/forum/read-10-23227...

Cheers
Detlef

Autor: Alexander (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hast du einen analogen Filter (Anti-Aliasing) mit einer Schnittfrequenz
von 100Hz vorgeschalten? Ansonsten ist deine Aussage Quark.

Hier ein Entwurf mit Matlab:

w0 = 50/200; bw = w0/35; [num,den] = iirnotch(w0,bw)

num =
    0.9889   -1.3985    0.9889

den =
    1.0000   -1.3985    0.9778

Im Anhang hast du den Frequenzgang. Oben sind die Koeffizienten des
Zählers (num) und unten die des Nenners (den).

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
wie sieht das in C aus ?

Autor: Alexander (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Zeig mir, dass du es verstanden hast und schreibe mir einfach die
Differenzengleichung hin, welche zu implementieren ist. Dann lasse ich
mich vielleicht zu einem C-Programm hinreißen :-)

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
y(k)=a1 x(k) + a2 x(k-1) + a3 x(k-2) + b2 y(k-1) + b3 y(k-2)

das ist die funktion die ich implementieren möchte ! muss ich jetzt
einfach die genannten koeffizienten verwenden, damit ich den filter
realisiert habe ?

Autor: Alex (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Genau, du legt dir zunächst zwei kleine Ringpuffer für die Ein- und
Ausgangswerte an. Das kann man einfach realisieren, indem man Arrays
nimmt und sich mittels eines Pointer (oder einfach Index) merkt,
welches Element gerade das k-te ist. Dann brauchst du z.B. zwei
for-Schleifen, in welchen du jeweils alle Elemente der Arrays (Ein-
bzw. Ausgangswerte) mit den jeweiligen Koeffizienten (a's bzw. b's)
multiplizierst und aufsummierst.

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
ok, dann werde ich das so machen:

a1=0.9889
a2=-1.3985
a3=0.9889

b1=1 (entfällt)
b2=-1.3985
b3=0.9778

die eigentliche realisierung in C weiss ich ! ich hab das nur noch nie
mit nem filter gemacht !!

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
kann mir jmd noch kurz beschrieben wie ich auf die koeffs komme ?
das wäre zu nett ..danke !!

Autor: Alex (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
So sollte es klappen, ist halt immer vorteilhaft, ein Werkzeug wie
Matlab zur Hand zu haben :)

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Halt, nicht weglaufen Alex ;) wie komme ich auf die Koeffizienten ?

Autor: Alexander (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Wie gesagt habe ich Matlab verwendet:
http://www.technion.ac.il/guides/matlab/toolbox/fi...

Interessiert dich wirklich der Hintergrund, wie der Entwurf abläuft? Da
wird es denke ich ein wenig komplexer.

Im wesentlichen ist es ja "nur" ein Bandpass. Die bildet man, indem
man einen Hoch- und einen Tiefpaß überlagert. Die Bandbreite im
Entwurfswerkzeug (bw) gibt an, wie weit deren Schnittfrequenzen
auseinander liegen. Stellt sich die Frage, mit welchem Ansatz man an
den Entwurf dieser beiden Filter geht - da gibt es dutzende?!

Das m-file macht im wesentlichen das folgende:

BW = BW*pi;
Wo = Wo*pi;

Gb   = 10^(-Ab/20);
beta = (sqrt(1-Gb.^2)/Gb)*tan(BW/2);
gain = 1/(1+beta);

num  = gain*[1 -2*cos(Wo) 1];
den  = [1 -2*gain*cos(Wo) (2*gain-1)];

Auf welchen Theorien dieser Ansatz basiert kann ich dir jedoch nicht
verraten. Als Referenz ist angegeben:
[1] Sophocles J. Orfanidis, Introduction To Signal Processing
Prentice-Hall 1996.

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
hmmm..das hilft nicht wirklich weiter ....
gibts andere Lösungsansätze ?

Autor: Andreas Auer (aauer1) Benutzerseite
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hi

Also, wenns hier jetzt darum geht, was hinter dem Entwurf von digitalen
Filtern steckt, dann kann ich das Buch "Signalverarbeitung" von
Oppenheim und Schäfer empfehlen. Wir haben das Buch in einer Vorlesung
verwendet an der TU verwendet.

Bezüglich dem Filterdesign... ich hab vor kurzem mal so richtig "from
the ground up" ein digitales Butterworth Filter 6. Ordnung auf einem
Mega128 implementiert (IIR).

Es gibt da verschiedene Ansätze... für IIR Filter werden meist nach
analogen Filtern entworfen. D.h. du hast irgenwelche Kenngrößen, die
deine Filterkennlinie (im Frequenzbereich) einhalten soll (z.B. im
Durchlassbereich eine Verstärkung von 0dB... im Sperrbereich z.B. -30dB
und einen mehr oder weniger breiten Übergangsbereich).
Danach gibts einige Ansätze...
* Beim ersten (weiß leider momentan den Begriff dafür nicht mehr - war
irgendwas mit "linearer...") wird die Gleichung für ein bestimmtes
Filter (in meinem Fall wäre es ein Butterworth Tiefpass gewesen)
aufgelöst und man erhält z-transformierte Übertragungsfunktion. Aus
dieser Übertragungsfunktion kannst du dann die Koeffizienten für deine
Filterstruktur herauslesen.

* Problem bei diesem Ansatz ist, dass durch diese "Annäherung" des
Filters Aliasing Effekte auftreten. D.h. bei einem Tiefpassfilter
überlagern sich die Kennlinien im Frequenzbereich. Im Durchlassbereich
des TP ist das weniger schlimm. Im Sperrbereich könnte es aber sein,
dass dadurch die geforderte Unterdrückung nicht mehr eingehalten wird.

* Der 2. Ansatz (den ich verwendet habe) nennt sich "Bilineare
Transformation". Diese Transformation umgeht das oben genannte Problem
indem es die gesamte Frequenzachse so skaliert, dass sie nur ein
einziges Mal auf dem Einheitskreis in der z-Ebene abgebildet wird.
Dadurch entsteht im Endergebnis kein Aliasing mehr.
Die vorgehensweise entspricht nach der Skalierung in etwa dem obigen
Schema... erst wird die Übertragungsfunktion gelöst und danach daraus
die Koeffizienten berechnet.

Wie gesagt... in dem oben gennanten Buch ist das alles sehr gut
beschrieben! Und... sorry für den etwas langen Beitrag, aber das Thema
Signalverarbeitung ist einfach sehr interessant, wenn auch nicht
unbedingt einfach!!!

mfg
Andreas

Autor: Dieter (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
hallo, ich steige gerade in Matlab ein und wollte ebenfalls einen
Notchfehler realisieren.

könnte das hinhaun?

clear;
clc;
% gestörtes Signal laden
load EKG_gestoert.mat;
% ===========================================================
% =             Notchfilter                                 =
% ===========================================================
% Koeffizienten für Notchfilter berechnen
w0 = 50/200;    %Annahme 200Hz Abtastfrequenz
bw = w0/35;
[num,den] = iirnotch(w0,bw);

a0=num(1);
a1=num(2);
a2=num(3);
b1=den(2);
b2=den(3);
y=zeros(1,5002);
% Filtereingangswerte
x(1)=0;
x(2)=0;
for n=1:1:5000
    x(n+2)=EKG_gestoert(n,2);
end
y(1)=0;
y(2)=0;
for n=3:1:5000
    y(n)=a0*x(n)+a1*x(n-1)+a2*x(n-2)-b1*y(n-1)-b2*y(n-2);
end
for n=1:1:5000
    EKG_notchfilter(n)=y(n+2);
end

Autor: Detlef A (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Das könnte hinhaun, wenn die Koeffizienten in der richtigen Reihenfolge
sind. iirnotch kenn ich nicht als Standard Matlab, ich weiß nicht was
das zurückliefert.

Aber: Extrem ungewöhnlich, die Vorwärtskoeffizienten heißen überall
'b', die Rückwärtsk. heißen 'a', bei Dir ist es umgekehrt.

Das ganze läßt sich auch so auf eine Zeile zusammenschnurren:

EKG_notchfilter=filter(num,den,EKG_gestoert);

das ist in Matlab eingebaut.

BTW: '1:1:5000' ist äquivalent zu '1:5000', Vektoren kannse in eins
zuweisen, z.B. EKGbla=y(3:length(y));, ohne 'for', Matlab ohne 'for'
is besser!

Cheers
Detlef

Autor: Dieter (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
ich dank dir vielmals Detlef,
wenn das mit einer Zeile geht, dann brauch ich mir je keine weiteren
Sorgen um die Implementierung der Differenzengleichung zu machen.

die iirnotch-Funktion wurde weiter oben im thread geklärt.
mittlerweile habe ich schon ein paar andere Filter umgesetzt,
an denen man die wirkung besser am graph erkennen kann.
da weiss man wenigstens, ob man's richtig gemacht hat oder nicht.
der notchfilter macht irgendwie ziemlich wenig.

danke für den tip mit dem vektor zuweisen,
kann ich auch mit einem kurzen befehl eine spalte oder zeile einer
matrix einem vektor zuweisen?

Autor: Dieter (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
hmm, ich glaub das mit dem Notchfilter haut doch nicht hin.
durch die anderen Filterrungen hab ich fast alles raus bis auf diese
schöne 50Hz Sinusfrequenz die das EKG noch überlagert.
die müsst er eigentlich glätten, aber passiert fast garnix.

Autor: Thomas W. (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Poste doch einfach mal die .mat Datei, dann sehen wir auch was du
meinst.

Autor: Detlef A (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hi,

ah, iirnotch hatte ich nicht gesehen. die krypten da etwas rum 'where
1.0 corresponds to pi radians per sample in the frequency range.' und
meinen vermutlich, daß 1 der halben samplefrequenz entspricht, Deine
Zeile also zu 50/100 ändern. BTW: freqz zeigt Dir den Frequenzgang.
Vorsicht, Bandbreite nich zu schmal,sonst kriegste die Frequenz nicht
,wenn etwas daneben. Zitiere mich selbst hier zu zweiten Mal:
http://www.mikrocontroller.net/forum/read-10-23227..., ist
nen notch per hand.

z.B.
>> a=ones(10,10);
>> v=1:10
>> a(3,:)=v
>> a(:,4)=v'
>> v=a(2,:);
>> v=a(:,7);
baut Vektoren inne Matrix ein und aus.
'If you got the ':', you got Matlab' nen Sruch von Mathworks
selbst.

Cheers
Detlef

Autor: Dieter (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
puuh es lag tatzächlich an dem
w0 = 50/200;
ich hab blöderweise die falsche abtastfrequenz eingesetzt. ich war der
festen überzeugung, die wäre 200Hz gewesen. dann hätte es tatsächlich
w0 = 50/100; heissen müssen, da wo = f0/(fA/2);
sie ist allerding 500Hz
also
w0 = 50/250;
und das ganze funktioniert.

jetzt nach dem FIR-Tiefpass(Hamming), Notchfilter und Siganal-Averaging
sieht das EKG schon ganz gut aus.

jetzt muss ich nur noch ein Wienerfilter hinbekomm.
doch da muss ich mich erstmal belesen, was das nun wieder ist.

Autor: Dieter (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
so, hab alles hinbekommen.
ich hätte noch eine matlabfrage, wo ich schonmal dabei bin.

ich habe jetzt jeden filter als funktion
und ein m-file in dem ich diese aufrufe.
der hintergrund ist, dass ich die reihenfolge der filter schnell
verändern kann.
nun frage ich mich gerade, ob ich das ganze auch mit einer grafischen
oberfläche realisieren kann.

dass man sich den filter oder die reihenfolge sozusagen aussuchen kann?

Autor: Markus (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Mal eine ganz andere Frage:

Wie finde ich solche tollen Seiten? Wurde gepostet von Alexander:

http://www.technion.ac.il/guides/matlab/toolbox/fi...

Dort gibt es ja einige schoene Beschreibungen zu Matlab. Sehr toller
Link. Vielen Dank ;).

Autor: Andreas Schwarz (andreas) (Admin) Benutzerseite Flattr this
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Das ist die offizielle Matlab-Hilfe, die findest du auch direkt bei
Mathworks.

Autor: jojo (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hi ich will gerne vom analogen filter ins digitalen bereich
tranzformnieren. was soll ich machen?? also der analoge befehl für ein
cheby-filter lautet [B,A]=cheby1(grad,welligkeit,grenzfrequenz,'s').
mein problem ist es jetzt ein digitales Filter zuentwerfen. Bittel
helft mir

Autor: Michael (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Lass einfach das 's' bei derm Befehl weg. Aber kannst du auch einfach 
nochmals in der Hilfe unter -> Help butter <- nachschauen.

Autor: Matthias (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Der Befehl gibt dir doch bereits ein digitales Filter zurück 
(Koffezienten A und B). Tschebyscheff Filter sind digitale 
Approximationen der gänigen analogen Filter.

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]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [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.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

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