Obwohl der Entwurf von FIR-Filtern mit linearer Phase eine einfache Aufgabe ist, gilt dies sicherlich nicht für IIR-Filter, die in der Regel einen stark nichtlinearen Phasengang aufweisen, insbesondere um die Grenzfrequenzen des Filters. In diesem Artikel werden die Eigenschaften besprochen, die ein digitales Filter benötigt, um eine lineare Phase zu haben, und wie die Durchlassphase eines IIR-Filters modifiziert werden kann, um eine lineare Phase mit Allpass-Entzerrungsfilter zu erreichen.

Warum brauchen wir linearphasige Filter?

Digitale Filter mit linearer Phase haben den Vorteil, dass sie alle Frequenzkomponenten um den gleichen Betrag verzögern, d. h. sie erhalten die Phasenbeziehungen des Eingangssignals. Diese Erhaltung der Phase bedeutet, dass das gefilterte Signal die Form des ursprünglichen Eingangssignals beibehält. Diese Eigenschaft ist für Audioanwendungen essentiell, da die Signalform entscheidend für die Beibehaltung einer hohen Wiedergabetreue im gefilterten Audiosignal ist. Ein weiterer Anwendungsbereich, der dies erfordert, ist die biomedizinische EKG-Wellenformanalyse, da jegliche Artefakte, die durch den Filter eingeführt werden, als Herzanomalien fehlinterpretiert werden können.

Das folgende Diagramm zeigt die Filterleistung eines Tschebyscheff-Typ-I-Tiefpasses auf EKG-Daten – die Eingangswellenform (in blau dargestellt) wurde um 10 Samples (\(\small \Delta=10\)) verschoben, um die Gruppenlaufzeit des Filters annähernd zu kompensieren. Beachten Sie, dass das gefilterte Signal (in rot dargestellt) abgeschwächt, verbreitert und mit Oszillationen um den EKG-Peak versehen ist, was unerwünscht ist.

Abbildung 1: Ergebnis der IIR-Tiefpassfilterung mit Phasenverzerrung

Damit ein digitales Filter eine lineare Phase hat, muss seine Impulsantwort eine konjugiert-gerade oder konjugiert-ungerade Symmetrie um ihren Mittelpunkt haben. Dies ist für ein FIR-Filter leicht zu erkennen,

\(\displaystyle H(z)=\sum\limits_{k=0}^{L-1} b_k z^{-k}\tag{1} \)

Mit der folgenden Einschränkung für seine Koeffizienten,

\(\displaystyle b_k=\pm\, b^{\ast}_{L-1-k}\tag{2} \)

was dazu führt,

\(\displaystyle z^{L-1}H(z) = \pm\, H^\ast (1/z^\ast)\tag{3} \)

Die Analyse von Gl. 3 zeigt, dass die Wurzeln (Nullstellen) von \(\small H(z)\) auch die Nullstellen von \(\small H^\ast (1/z^\ast)\)sein müssen. Das bedeutet, dass die Wurzeln von \(\small H(z)\) in konjugierten reziproken Paaren auftreten müssen, d. h. wenn \(\small z_k\) eine Nullstelle von \(\small H(z)\) ist, dann muss auch \(\small H^\ast (1/z^\ast)\) eine Nullstelle sein.

Warum IIR-Filter keine lineare Phase haben

Ein digitales Filter wird als “bounded input, bounded output stable” oder “BIBO stable” bezeichnet, wenn jede gebundene Eingabe zu einer gebundenen Ausgabe führt. Alle IIR-Filter haben entweder Pole oder sowohl Pole als auch Nullstellen und müssen BIBO-stabil sein, d. h.

\(\displaystyle \sum_{k=0}^{\infty}\left|h(k)\right|<\infty \tag{4}\)

Dabei ist \(\small h(k)\) die Impulsantwort des Filters. Bei der Analyse von Gl. 4 sollte klar sein, dass das BIBO-Stabilitätskriterium nur dann erfüllt ist, wenn die Pole des Systems innerhalb des Einheitskreises liegen, da die ROC (Region der Konvergenz) des Systems den Einheitskreis einschließen muss. Folglich ist es ausreichend zu sagen, dass ein begrenztes Eingangssignal immer ein begrenztes Ausgangssignal erzeugt, wenn alle Pole innerhalb des Einheitskreises liegen.

Die Nullstellen hingegen sind nicht durch diese Anforderung eingeschränkt und können folglich überall auf der z-Ebene liegen, da sie die Systemstabilität nicht direkt beeinflussen. Daher kann eine Systemstabilitätsanalyse durchgeführt werden, indem zunächst die Wurzeln der Übertragungsfunktion (d. h. die Wurzeln der Zähler- und Nennerpolynome) berechnet werden und dann die entsprechenden Pole und Nullstellen auf der z-Ebene aufgetragen werden.

Wendet man die entwickelte Logik auf die Pole eines IIR-Filters an, kommt man nun zu einer sehr wichtigen Schlussfolgerung, warum IIR-Filter keine lineare Phase haben können.

Ein BIBO-stabiles Filter muss seine Pole innerhalb des Einheitskreises haben, und um eine lineare Phase zu erhalten, müsste ein IIR-Filter konjugierte reziproke Pole außerhalb des Einheitskreises haben, was es BIBO-instabil macht.

Basierend auf dieser Aussage scheint es nicht möglich zu sein, einen IIR so zu konstruieren, dass er eine lineare Phase hat. Wie unten beschrieben, können jedoch Phasenentzerrungsfilter verwendet werden, um den Phasengang des Durchlassbereichs zu linearisieren.

Phasenlinearisierung mit Allpassfiltern

All-pass phase linearisation filters (equalisers) (Equalizer) sind eine gut etablierte Methode, um den Phasengang eines Filters zu ändern, ohne den Betragsgang zu beeinflussen. Ein Allpassfilter zweiter Ordnung (Biquad) ist definiert als:

\( A(z)=\Large\frac{r^2-2rcos \left( \frac{2\pi f_c}{fs}\right) z^{-1}+z^{-2}}{1-2rcos \left( \frac{2\pi f_c}{fs}\right)z^{-1}+r^2 z^{-2}}\tag{5} \)

Dabei ist\(\small f_c\) die Mittenfrequenz, \(\small r\) der Radius der Pole und \(\small f_s\) die Abtastfrequenz. Beachten Sie, dass die Koeffizienten von Zähler und Nenner als spiegelbildliches Paar zueinander angeordnet sind.  Die Spiegeleigenschaft ist es, die dem Allpassfilter seine wünschenswerte Eigenschaft verleiht, nämlich dem Konstrukteur zu ermöglichen, den Phasengang zu verändern, während der Betragsgang über das gesamte Frequenzspektrum konstant oder flach bleibt.

Wenn man eine APF-Entzerrungskaskade (bestehend aus mehreren APFs) mit einem IIR-Filter kaskadiert, besteht die Grundidee darin, dass man nur den Phasengang im Durchlassbereich linearisieren muss. Die anderen Bereiche, wie z. B. das Übergangsband und das Sperrband, können ignoriert werden, da jegliche Nichtlinearitäten in diesen Bereichen für das Gesamtergebnis der Filterung von geringem Interesse sind.

Die Herausforderung

Die APF-Kaskade klingt wie ein idealer Kompromiss für diese Herausforderung, aber in Wahrheit ist ein erheblicher Zeitaufwand und eine sehr sorgfältige Feinabstimmung der APF-Positionen erforderlich, um ein akzeptables Ergebnis zu erzielen. Jeder APF hat zwei Variablen \(\small f_c\) und \(\small r\) , die optimiert werden müssen, was die Lösung verkompliziert. Erschwerend kommt hinzu, dass die Gruppenverzögerung (Latenz) des Gesamtfilters umso höher wird, je mehr APF-Stufen zur Kaskade hinzugefügt werden. Letzteres kann für schnelle Echtzeit-Regelsysteme, die sich auf die niedrige Latenzzeit eines IIR-Filters verlassen, problematisch werden.

Trotz dieser Herausforderungen ist der APF-Equalizer ein guter Kompromiss für die Linearisierung der Durchlassband-Phaseneigenschaften eines IIR-Filters.

Der APF-Entzerrer

Der ASN Filter Designer bietet Entwicklern eine sehr einfach zu bedienende grafische All-Phase-Equalizer-Oberfläche für die Linearisierung der Durchlassbandphase von IIR-Filtern. Wie unten zu sehen ist, ist die Oberfläche sehr intuitiv und erlaubt es dem Designer, APF-Filterpositionen schnell mit der Maus zu platzieren und fein abzustimmen. Das Werkzeug berechnet automatisch \(\small f_c\) und \(\small r\), basierend auf der Markerposition.

APF-Equalizer ASN Filter Designer

APF equaliser ASN Filter Designer

Wenn Sie mit der rechten Maustaste auf das Frequenzgangdiagramm oder auf einen vorhandenen Allpass-Design-Marker klicken, wird ein Optionsmenü angezeigt, wie links dargestellt.

Sie können bis zu 10 Biquads hinzufügen (nur professionelle Version).

Ein IIR mit linearer Phase im Durchlaßbereich

Wenn wir einen Equalizer entwerfen, der aus drei APF-Paaren besteht, und ihn mit dem Tschebyscheff-Filter aus Abbildung 1 kaskadieren, erhalten wir eine Filterwellenform, die eine viel schärfere Spitze mit weniger Dämpfung und Oszillation als der ursprüngliche IIR hat – siehe unten. Diese Verbesserung geht jedoch auf Kosten von drei zusätzlichen Biquad-Filtern (der APF-Kaskade) und einer erhöhten Gruppenlaufzeit, die nun auf 24 Samples im Vergleich zu den ursprünglichen 10 Samples angestiegen ist.

IIR lowpass filtering result with three APF phase equalisation filters
(minimal phase distortion)
Ergebnis der IIR-Tiefpassfilterung mit drei APF-Phasenausgleichsfiltern
(minimale Phasenverzerrung)

Der Frequenzgang sowohl des ursprünglichen IIR als auch des entzerrten IIR ist unten dargestellt, wobei die Gruppenverzögerung (in lila dargestellt) die durchschnittliche Verzögerung des Filters ist und eine einfachere Methode zur Beurteilung der Linearität darstellt.

IIR without equalisation cascade
IIR ohne Entzerrungskaskade

IIR with equalisation cascade
IIR mit Entzerrungskaskade

Beachten Sie, dass die Gruppenlaufzeit des entzerrten IIR-Durchlassbereichs (rechts dargestellt) fast flach ist, was bestätigt, dass die Phase tatsächlich linear ist.

Automatische Code-Generierung für Arm-Prozessorkerne über CMSIS-DSP

Die automatische Code-Generierung des ASN Filter Designers erleichtert den Export eines entworfenen Filters auf Cortex-M Arm basierte Prozessoren über das CMSIS-DSP Software-Framework. Die integrierten Analyse- und Hilfefunktionen des Tools unterstützen den Designer bei der erfolgreichen Konfiguration des Designs für den Einsatz.

Vor der Generierung des Codes müssen die IIR- und Entzerrungsfilter (d. h. H1- und Heq-Filter) zunächst zu einer H1-Filterstruktur (Hauptfilter) für den Einsatz reoptimiert (zusammengeführt) werden. Das Optionsmenü befindet sich unter der Registerkarte P-Z in der Haupt-Benutzeroberfläche.

Alle Entwürfe von Fließkomma-IIR-Filtern müssen auf Single-Precision-Arithmetik und entweder auf einer Direct Form I- oder Direct Form II-Transposed-Filterstruktur basieren. Die Direct Form II Transposed-Struktur wird aufgrund ihrer höheren numerischen Genauigkeit für die Fließkomma-Implementierung befürwortet.

Die Einstellungen für Quantisierung und Filterstruktur finden Sie unter der Registerkarte Q (wie links dargestellt). Durch Einstellen von Arithmetic auf Single Precision und Structure auf Direct Form II Transposed und Klicken auf die Schaltfläche Apply wird die hier betrachtete IIR für das Software-Framework CMSIS-DSP konfiguriert.

Wählen Sie das Arm CMSIS-DSP-Framework aus der Auswahlbox im Filterübersichtsfenster aus:

Der automatisch generierte C-Code auf Basis des CMSIS-DSP-Frameworks für die direkte Implementierung auf einem Arm-basierten Cortex-M-Prozessor ist unten dargestellt:

The ASN Filter Designer's automatic code generator generates all initialisation code, scaling and data structures needed to implement the linearised filter IIR filter via Arm's CMSIS-DSP library.

Der automatische Code-Generator des ASN Filter Designers generiert den gesamten Initialisierungscode, die Skalierung und die Datenstrukturen, die für die Implementierung des linearisierten IIR-Filters über die CMSIS-DSP-Bibliothek von Arm benötigt werden.

Was wir gelernt haben

Die Wurzeln eines linearphasigen Digitalfilters müssen in konjugierten reziproken Paaren auftreten. Obwohl dies für ein FIR-Filter kein Problem ist, wird es für ein IIR-Filter undurchführbar, da die Pole sowohl innerhalb als auch außerhalb des Einheitskreises liegen müssten, was das Filter BIBO-instabil macht.

Der Durchlassbereich-Phasengang eines IIR-Filters kann mit Hilfe einer APF-Entzerrungskaskade linearisiert werden. Der ASN Filter Designer bietet dem Entwickler über eine sehr einfach zu bedienende, grafische Allpass-Phasenentzerrer-Oberfläche alles, was er braucht, um eine geeignete APF-Kaskade nur mit der Maus zu entwerfen!

Das linearisierte IIR-Filter kann über den automatischen Code-Generator unter Verwendung der optimierten CMSIS-DSP-Bibliotheksfunktionen von Arm für den Einsatz auf jedem Cortex-M-Mikrocontroller exportiert werden.

 

 

Download demo now

Licencing information

Bei der EKG-Signalverarbeitung ist die Rauschunterdruckung von 50/60Hz-Powerline-Interferenzen aus empfindlichen, informationsreichen biomedizinischen EKG-Wellenformen eine herausfordernde Aufgabe! Die Herausforderung wird durch die Anpassung an die Auswirkungen des EMG, wie z.B. die Bewegung von Gliedmaßen/Rumpf des Patienten oder sogar die Atmung, noch komplizierter. Ein traditioneller Ansatz, der von vielen angenommen wird, ist die Verwendung eines IIR-Kerbfilters zweiter Ordnung:

\(\displaystyle H(z)=\frac{1-2cosw_oz^{-1}+z^{-2}}{1-2rcosw_oz^{-1}+r^2z^{-2}}\)

wobei \(w_o=\frac{2\pi f_o}{fs}\) die Mittenfrequenz, \(f_o\) die Kerbe und \(r=1-\frac{\pi BW}{fs}\) die Bandbreite (-3dB-Punkt) der Kerbe steuert.

Wo liegt die Herausforderung?

Wie oben gesehen, ist \(H(z) \) einfach zu implementieren, aber die Schwierigkeit liegt darin, einen optimalen Wert für \(r\), zu finden, da eine wünschenswerte scharfe Kerbe bedeutet, dass die Pole nahe am Einheitskreis liegen (siehe rechts).

In Gegenwart von stationären Störungen, z.B. wenn der Patient absolut ruhig ist und die Auswirkungen der Atmung auf die Sensordaten minimal sind, ist dies möglicherweise kein Problem.

Betrachtet man jedoch die Auswirkungen von EMG auf die erfasste Wellenform (eine viel realistischere Situation), verursacht die Rückkopplung (Pole) des IIR-Filters ein Klingeln auf der gefilterten Wellenform, wie unten dargestellt:

Contaminated ECG with non-stationary 50Hz powerline interference (FIR filtering), ECG sigal processing, ECG DSP, ECG measurement

Verschmutztes EKG mit nicht-stationären 50Hz-Stromleitungsstörungen (IIR-Filterung)

Wie oben zu sehen ist, wurde zwar ein Großteil der 50Hz-Stromleitungsstörungen beseitigt, aber es gibt immer noch ein deutliches Ringing um die Hauptspitzen herum (gefilterter Ausgang in rot dargestellt). Dieses Ringing ist für viele biomedizinische Anwendungen unerwünscht, da wichtige kardiale Informationen wie das ST-Segment nicht eindeutig analysiert werden können.

Der Frequenzgang des IIR, der zur Filterung der obigen EKG-Daten verwendet wird, ist unten dargestellt.

IIR notch filter frequency response, ECG signal processing, ECG DSP, ECG  measurement

Frequenzgang des IIR-Kerbfilters

Die Analyse des Diagramms zeigt, dass die Gruppenlaufzeit (group delay, oder durchschnittliche Verzögerung) des Filters nichtlinear ist, aber in den Durchlassbereichen fast Null ist, was bedeutet, dass keine Verzerrung vorliegt. Die Gruppenlaufzeit bei 50 Hz steigt auf 15 Abtastwerte an, was die Quelle des Schwingens ist – wobei die Gruppenlaufzeit umso größer ist, je näher sich die Pole am Einheitskreis befinden.

ASN FilterScript bietet Designern die Funktion notch(), die eine direkte Implementierung von H(z), ist, wie unten gezeigt:

ClearH1;  // clear primary filter from cascade
ShowH2DM;   // show DM on chart

interface BW={0.1,10,.1,1};

Main()

F=50;
Hd=notch(F,BW,"symbolic");
Num = getnum(Hd); // define numerator coefficients
Den = getden(Hd); // define denominator coefficients
Gain = getgain(Hd); // define gain

Savitzky-Golay FIR-Filter

Eine Lösung für das oben erwähnte Ringing sowie die Rauschunterdrückung kann durch einen Savitzky-Golay-Tiefpass-Glättungsfilter erreicht werden. Diese Filter sind FIR-Filter und haben daher keine Rückkopplungskoeffizienten und kein Ringing!

Savitzky-Golay-(Polynom-)Glättungsfilter oder Glättungsfilter der kleinsten Quadrate sind Verallgemeinerungen des FIR-Durchschnittsfilters, die den Hochfrequenzanteil des gewünschten Signals besser erhalten können, auf Kosten der Entfernung von nicht so viel Rauschen wie ein FIR-Durchschnitt. Die besondere Formulierung von Savitzky-Golay-Filtern bewahrt verschiedene Momentenordnungen besser als andere Glättungsmethoden, die dazu neigen, Spitzenbreiten und -höhen besser als Savitzky-Golay zu erhalten. Als solche sind Savitzky-Golay-Filter sehr gut für biomedizinische Daten, wie z.B. EKG-Datensätze, geeignet.

Eliminierung der 50Hz-Powerline-Komponente

Beim Entwurf eines Savitzky-Golay-Filters 18. Ordnung mit einer Polynomanpassung 4. Ordnung (siehe den Beispielcode unten) erhalten wir einen FIR-Filter mit einer Nullverteilung, wie rechts dargestellt. Da wir jedoch die 50Hz-Komponente vollständig eliminieren möchten, kann der P-Z-Editor des Tools verwendet werden, um ein Nullpaar (grün dargestellt) auf genau 50Hz zu schieben.

Der resultierende Frequenzgang ist unten dargestellt, wobei zu erkennen ist, dass bei genau 50 Hz eine Kerbe vorhanden ist und die Gruppenlaufzeit von 9 Proben (in violett dargestellt) über das Frequenzband konstant ist.

FIR  Savitzky-Golay filter frequency response, ECG signal processing, ECG DSP, ECG measurement

FIR Savitzky-Golay-Filter-Frequenzgang

Wir lassen den verunreinigten EKG-Datensatz durch unseren optimierten Savitzky-Golay-Filter laufen und passen ihn an die Gruppenverzögerung an, die wir erhalten:

Contaminated ECG with non-stationary 50Hz powerline interference (FIR filtering), ECG signal processing, ECG digital filter, ECG filter designa

Verschmutztes EKG mit nicht-stationären 50Hz-Stromleitungsstörungen (FIR-Filterung)

Wie man sieht, gibt es keine Anzeichen für ein Klingeln, und die ST-Segmente sind jetzt für die Analyse deutlich sichtbar. Beachten Sie auch, wie der Filter (rot dargestellt) das Messrauschen reduziert hat, was die praktische Anwendbarkeit von Savitzky-Golay-Filtern für die biomedizinische Signalverarbeitung unterstreicht.

Ein Savitzky-Golay kann in ASN FilterScript über die Funktion savgolay() wie folgt entworfen und optimiert werden:

ClearH1;  // clear primary filter from cascade

interface L = {2, 50,2,24};
interface P = {2, 10,1,4};

Main()

Hd=savgolay(L,P,"numeric");  // Design Savitzky-Golay lowpass
Num=getnum(Hd);
Den={1};
Gain=getgain(Hd);

Bereitstellung

Dieser Filter kann nun über den automatischen Code-Generator des Werkzeugs in einer Vielzahl von Domänen eingesetzt werden, was einen schnellen Einsatz in Matlab-, Python- und eingebetteten Arm Cortex-M-Geräten ermöglicht.

Wie wählen wir den besten Typ von Digitalfilter für unsere Anwendung aus? Und was sind die Unterschiede zwischen einem IIR-Filter und einem FIR-Filter?