Digitale Filter Martin Schlup 8. Mai 2012

1. Filterstrukturen Dieser Beitrag ist eine kurz gehaltene Einführung in die Darstellung zeitdiskreter Systeme und soll einige elementare Hinweise geben, wie digitale Filter mittels Matlab synthetisiert und grob analysiert werden können. Auf Quantisierungsaspekte (u. A. Wertquantisierung der Filterparameter) wird hier bewusst nicht eingegangen.

1.1. MA-Systeme oder FIR-Filter Bei dieser Klasse von zeitdiskreten LTI-Systemen1 wird das Ausgangssignal y[k] direkt aus dem Eingangssignal u[k] und seiner vorausgehenden Werten (Vorläufer) u[k − i] gebildet. Weil u. A. die Bildung eines gleitenden Mittelwertes (moving average) diese Struktur aufweist, werden diese Systeme MA-Systeme genannt. Da ausserdem vorangehende Ausgangssignale des Systems nicht für das aktuelle Ausgangssignal benutzt werden, heissen diese Systeme auch nicht-rekursiv. Die Differenzengleichung dieser Systeme besitzt folgende Struktur: y[k] = b0 u[k] + b1 u[k − 1] + · · · + bm u[k − m]

(1)

Die Ordnung des Systems entspricht dabei der Grösse des benötigten Schieberegisters welches benötigt wird, um die Vorgänger des Eingangssignals zur Bildung des Ausgangssignals zu speichern. Das System besitzt demzufolge ein endliches Gedächtnis der „Länge“ m. Ferner besitzen MA-Systeme nur Nullstellen 2 und keine von Null verschiedenen Pole 3 , daher werden sie auch als all zero systems bezeichnet. MA-Systeme besitzen eine endliche Stossantwort (finite impulse response)4 deren Dauer m + 1 Takte beträgt. Diese lautet in Operator-Schreibweise: H(q −1 ) = b0 + b1 q −1 + · · · + bm q −m 1 2

3 4

(2)

linear time invariant Als Nullstellen bezeichnet man die Wurzeln des Polynoms b0 + b1 z −1 + · · · + bm z −m , bzw. die m (möglicherweise komplexen) Lösungen zN der Gleichung b0 z m + b1 z m−1 + · · · + bm−1 z + bm = 0. Stichwort: z-Transformation cf. nächster Abschnitt Die Stossantwort eines zeitdiskreten Systems ist die Antwort auf den diskreten Einheitsimpuls: u[k] = δ[k] → y[k] = h[k]. In Operatorschreibweise lautet sie hier: H(q −1 ) = h[0] + h[1]q −1 + · · · + h[m]q −m . Das Symbol q steht dabei für eine Verzögerung um das Abtastintervall Ts : q · x[k] = x[k + 1]. Entsprechend entspricht der Operator q −1 einer zeitlichen Vorverlegung um Ts : q −1 · x[k] = x[k − 1].

1

1. Filterstrukturen Offensichtlich entsprechen die Koeffizienten bk eines MA-Systems gerade den einzelnen Werten der Stossantwort h[k]. Daher werden MA-Systeme in der digitalen Datenverarbeitung auch als FIR-Filter bezeichnet. FIR-Filter weisen unter der Bedingung einer symmetrischen Stossantwort einen linearen Phasengang auf, d. h. harmonische Signale werden unabhängig von ihrer Frequenz, um dasselbe Zeitintervall verzögert (konstante Gruppenlaufzeit): Diese Verzögerungszeit beträgt dann m 2 Ts .

1.2. AR-/ARMA-Systeme oder IIR-Filter Eine andere Struktur ergibt sich, wenn das Ausgangsignal ausser vom Eingangssignal noch von n Vorläufern des Ausgangssignals abhängt. Weil das Ausgangssignal dabei rekursiv gebildet wird (auto regressive), werden diese Systeme AR-Systeme oder auch rekursiv genannt. Die Differenzengleichung dieser Systeme besitzt folgende Struktur: y[k] + a1 y[k − 1] + · · · + an y[k − n] = b0 u[k]

(3)

Die MA- und die AR-Struktur können im Allgemeinen gemeinsam auftreten. Dann hat das resultierende System ARMA-Struktur: y[k] + a1 y[k − 1] + · · · + an y[k − n] = b0 u[k] + b1 u[k − 1] + · · · + bm u[k − m]

(4)

Die Ordnung des allgemeinen Systems entspricht dabei der maximalen Verzögerung mTs oder nTs des Eingangs- oder Ausgangssignals, je nachdem welcher der beiden Werte grösser ist. Das System besitzt demzufolge ein Gedächtnis der „Länge“ max(m, n). AR- und ARMA-, bzw. rekursive Systeme besitzen eine unendliche Stossantwort (infinite impulse response). Desswegen werden sie in der digitalen Datenverarbeitung auch als IIRFilter bezeichnet. Im Allgemeinen weisen diese Systeme keinen linearen Phasengang auf, d. h. deren Gruppenlaufzeit ist frequenzabhängig. Dies führt durch die verschiedenen frequenzabhängigen Signallaufzeiten auf die so genannte Phasenverzerrung. Die diskrete Stossantwort von AR-/ARMA-Systemen, bzw. von IIR-Filtern lautet in (normierter) Operator-Schreibweise: H(q −1 ) =

b0 1 + a1 q −1 + · · · + an q −n

(5)

H(q −1 ) =

b0 + b1 q −1 + · · · + bm q −m 1 + a1 q −1 + · · · + an q −n

(6)

Um die einzelnen Werte der Stossantwort h[k] aus den Gleichungen (5) und (6) zu erhalten, muss eine Polynomdivision durchgeführt werden. Dabei wird klar, dass die Stossantwort im Allgemeinen nicht endlich sein kann. AR-Systeme besitzen nur Pole 5 (all pole systems), ARMA-Systeme Pole und Nullstellen. Rekursive Systeme sind nur dann stabil 6 , wenn sich alle Pole innerhalb des Einheitskreises in der komplexen Zahlenebene befinden. MA-Systeme sind dagegen immer stabil. 5

6

Als Pole bezeichnet man die Wurzeln des Polynoms 1 + a1 z −1 + · · · + an z −n , bzw. die n Lösungen zP der Gleichung z n + a1 z n−1 + · · · + an−1 z + an = 0. Als stabil wird ein System dann bezeichnet, wenn für ein beliebiges, begrenztes Eingangssignal die Ausgangsgrösse ebenfalls begrenzt bleibt.

2

2. Eigenschaften

2. Eigenschaften 2.1. Frequenzgang Der Frequenzgang H(ω) eines zeitdiskreten Systems kann direkt aus der Stossantwort H(q −1 ) durch Ersetzen des Verschiebungsoperators q durch eiωTs erhalten werden: H(ω) =

b0 + b1 e−iωTs + · · · + bm e−miωTs 1 + a1 e−iωTs + · · · + an e−niωTs

(7)

Der Amplitudengang des Systems entspricht dem Betrag |H(ω)| und der Phasengang dem Argument 6 H(ω) des Frequenzgangs. Amplituden und Phasengang eines zeitdiskreten Systems sind kontinuierlich im Frequenzverlauf und periodisch 7 in der Abtastfrequenz fs , so dass es genügt den Frequenzgang nur bis zur halben Abtastfrequenz (Nyquistfrequenz) anzugeben8 .

2.2. Gruppenlaufzeit Die Gruppenlaufzeit (group delay) eines linearen Systems wird als negative Ableitung des Phasengangs ϕ(ω) = 6 H(ω) nach der Kreisfrequenz ω definiert: τGr = −

dϕ(ω) dω

(8)

Die Gruppenlaufzeit ist bei AR- und ARMA-Systemen, bzw. IIR-Filtern im Allgemeinen abhängig von der Frequenz. Bei MA-Systemen, bzw. FIR-Filtern ist sie (unter bestimmten Bedingungen) unabhängig von der Frequenz. Sie ist ein Mass für die Verzögerung eines Signals durch das System. Daher auch das negative Vorzeichen in der Formel. Proportionalität zwischen der Phase und der Frequenz (linearer Phasengang) führt zu einer konstanten Gruppenlaufzeit, bzw. zu einer frequenzunabhängigen Verzögerung des Signals. Für frequenzabhängige Gruppenlaufzeiten bei Systemen mit Tiefpasseigenschaften, kann als Mass für die Signalverzögerung die Gruppenlaufzeit bei ω = 0 herangezogen werden: dϕ(ω) τGr ≈ − (9) dω ω=0

2.3. Vergleich FIR/IIR

Tabelle 1: Vergleich der wichtigsten Eigenschaften zwischen MA- und AR-/ARMA-Systemen Merkmale MA-System (FIR) AR- oder ARMA-System (IIR) Stossantwort endliche Dauer unendliche Dauer Phasengang linear in der Frequenz nichtlinear in der Frequenz Stabilität immer stabil, da nur Nullstellen nur stabil, falls Pole in Einheitskreis Ordnung um identische wesentlcih höher als IIR-Filter, kleiner als FIR-Filter, Spezifikation einzuhalten längere Einschwing- und Totzeit kürzere Einschwing- und Totzeit

7 8

i2π t

Dies lässt mit der 2πfs -Periodizität der komplexen Exponentialfunktion eiωt = e Ts erklären. Dabei weist der Amplitudengang eine gerade und der Phasengang eine ungerade Symmetrie auf.

3

A. Matlab-Hilfe: Digitale Filter (Synthese und Anwendung)

A. Matlab-Hilfe: Digitale Filter (Synthese und Anwendung) Für weitergehende Angaben sollte die Dokumentation für die Signal Processing Toolbox konsultiert werden. Für die Erstellung und die Analyse von Filtern gibt es in Matlab ein spezialisiertes Hilfsmittel: fdatool (Filter Design and Analysis Tool ). Für die Beschreibung dieses Werkzeugs wird auf das Handbuch der Signal Processing Toolbox verwiesen.

A.1. IIR-Filter (infinite impulse response) IIR- oder Rekursiv-Filter können steile übergangsflanken aufweisen, haben aber im Allgemeinen einen nichtlinearen Phasengang. Solche Filter verzerren z. B. Impulse durch die frequenzabhängige Gruppen- und Phasenlaufzeit. Mit dieser Filterart können TP, HP, BP, BS, Integratoren und Allpässe realisiert werden. A.1.1. Synthesefunktionen für Tiefpassfilter Bei diskreten Systemen werden in Matlab Frequenzen immer in normierter Form angegeben. Als Bezugsgrösse wird die Nyquist-Frequenz f2s , d. h. die halbe Abtastfrequenz benutzt9 . Listing 1: Matlab-Funktionen zur Bestimmung der Koeffizienten für diverse IIR-Filter [n,wn]=buttord(wnp,wns,rp,rs); [n,wn]=cheb1ord(wnp,wns,rp,rs); [n,wn]=cheb2ord(wnp,wns,rp,rs); [n,wn]=ellipord(wnp,wns,rp,rs);

% % % %

fuer fuer fuer fuer

Butterworth−Filter Tschebyscheff−Filter 1. Art Tschebyscheff−Filter 2. Art elliptisches oder Cauer−Filter

[b,a]=butter(n,wn); % Butterworth−Filter (monoton abklingender Ampli− % tudengang, halbwegs konstante Gruppenlaufzeit) [b,a]=cheby1(n,rp,wn); % Tschebyscheff−Filter 1. Art % (ripple im Durchlassbereich, peak to peak) [b,a]=cheby2(n,rs,wn); % Tschebyscheff−Filter 2. Art % (ripple im Sperrbereich) [b,a]=ellip(n,rp,rs,wn);% elliptisches oder Cauer−Filter (Standardfilter, % kleinste Ordnung, stark verzerrte Gruppenlaufzeit)

Tabelle 2: Legende zu den Matlab-Funktionen aus Listing 1 wnp wns rp rs n wn b a

9

f

dB Spezifikation TP-Durchlassbereich (auf fs /2 normierte Bandbreite: wnp = −3 fs /2 ) Spezifikation Sperrbereich, normiert, reicht von wns bis 1 (≡ fs /2), wns > wnp spezifizierte maximale Welligkeit (ripple) im Durchlassbereich (Passband) in dB spezifizierte minimale Dämpfung (attenuation) im Sperrbereich (stopband ) in dB Filterordnung fg auf die Nyquist-Frequenz fs /2 normierte Grenzfrequenz: wn = fs /2 Vektor der Filterkoeffizienten bk gemäss Gleichung (5) oder (6) Vektor der Filterkoeffizienten ak gemäss Gleichung (5) oder (6)

Bei normierten Frequenzen spielt es keine Rolle, ob dabei die Frequenz f oder die Kreisfrequenz ω = 2πf benutzt wird.

4

A. Matlab-Hilfe: Digitale Filter (Synthese und Anwendung) A.1.2. Synthesefunktionen für Hochpass-, Bandpass- und Bandsperrfilter (HP, BP, BS) Für die Synthese von HP-, BP- und BS-Filtern werden die selben Funktionen wie für die TPFilter verwendet. Die dazu notwendige Syntax ist an Hand der Funktion butter im Listing 2 angegeben: Listing 2: Matlab-Steuerbefehle um HP-, BP- und BS-Filter an Hand der Funktion butter zu erzeugen [b,a]=butter(n,wn,'high'); % Hochpassfilter (highpass) [b,a]=butter(n,[wn1,wn2]); % Bandpassfilter (bandpass) [b,a]=butter(n,[wn1,wn2],'stop'); % Bandsperrfilter (bandstop)

Analoge Eingabeformate gelten für die anderen Filtertypen und für die Funktionen buttord, cheb1ord, . . . usw., falls HP-, BP- und BS-Filter spezifiziert werden sollen. A.1.3. Synthese eines IIR-Filters mit beliebigem Amplitudengang

Listing 3: Matlab-Funktion um IIR-Flter mit beliebigem Amplitudengang zu erzeugen [b,a]=yulewalk(n,wn,amp);

% Yule−Walker Methode (least square fit)

Tabelle 3: Legende zu Listing 3 wn amp

Vektor mit normierten Frequenzstützwerten Der erste Frequenzwert muss 0, der letzte 1 (≡ fs /2) sein. Vektor mit gewünschten Amplitudengangswerten (Faktor) bei den entsprechenden Frequenzen

A.2. FIR-Filter (finite impulse response) Nichtrekursive, FIR- oder Transversal-Filter (tapped delay line) haben weniger steile übergangsflanken als IIR-Filter gleicher Ordnung. Ihr Phasengang ist aber im Allgemeinen linear, d. h. ihre Gruppenlaufzeit ist konstant. Mit dieser Filterart können TP, HP, BP, BS, Differentiatoren und Hilbert-Transformatoren realisiert werden. A.2.1. Synthesefunktionen für FIR-Tiefpassfilter

Listing 4: Matlab-Funktionen zur Bestimmung der Filterkoeffizienten für diverse FIR-Filter b=fir1(n,wn);

% % b=fir2(n,wn,m); % b=firpm(n,wn,m);% % % %

FIR−Filterentwurf mit Fenstermethode, lineare Phase Parameter koennen mit der Funktion "kaiserord" bestimmt werden punktweise vorgegebener Frequenzgang lineare Phase mit Gruppenlaufzeit: Ts*n/2, equiripple (Parks−McClellan, Standardverfahren) die benoetigte Ordnung von "firpm" kann mit der Funktion "firpmord" bestimmt werden

5

A. Matlab-Hilfe: Digitale Filter (Synthese und Anwendung)

A.3. Anwendung und Eigenschaften A.3.1. Signalfilterung Die Filterung eines Signals u[k] mit einem Digitalfilter mit den Koeffizienten bi (Koeffizienten des Zählerpolynoms) und ai (Nennerpolynom) erfolgt mit den Funktionen filter oder filtfilt. Listing 5: Filter-Aufruf y=filter(b,a,u); y=filtfilt(b,a,u);

% direkte Filterung des Signals u[k] % Filterung des Signals u[k] ohne Phasenverzerrung

Bemerkung: Die Funktion filtfilt filtert das Signal vorwärts und rückwärts, so dass sich die Signalverzögerung und die Phasenverzerrungen auf Grund des Phasengangs aufheben. Durch das zweimalige Filtern verdoppelt sich die Filterordnung. Diese Art der Filterung kann selbstverständlich nur bei einem ganzheitlich vorhandenen Signal angewendet werden und eignet sich nicht zur „Echtzeitverarbeitung“. A.3.2. Frequenzgang und Gruppenlaufzeit Der Frequenzgang (Amplituden- und Phasengang) eines diskreten Systems kann mit der MatlabFunktion freqz, die Gruppenlaufzeit mit grpdelay bestimmt werden, wie dies im folgenden Listing exemplarisch dargestellt wird: Listing 6: Bestimmung von Frequenzgang und Gruppenlaufzeit [H,w]=freqz(b,a,N); Ha=abs(H); % Amplitudengang (Betrag von H) Hp=180*unwrap(angle(H))/pi; % Phasengang in Grad und 2pi−Spruenge entfernen [gd,f]=grpdelay(b,a,N,fs); gd=grpdelay(b,a,f,fs);

Tabelle 4: Legende zu Listing 6 gd N f

Gruppenlaufzeit in Anzahl Abtastschritten der Dauer Ts = 1/fs Anzahl Frequenzwerte zwischen 0 und fs /2, sollte eine Zweierpotenz sein, z. B.: N = 512 = 29 Frequenzwerte

6

A. Matlab-Hilfe: Digitale Filter (Synthese und Anwendung)

A.4. Beispiele A.4.1. Beispiel IIR-TP-Filter

Listing 7: Matlab-Beispiel IIR TP-Filter % Filter−Spezifikationen fs=2e3; % Abtastfrequenz in Hz fg=100; % Grenzfrequenz in Hz ("−3dB"−Frequenz) wp=fg*2/fs; % Durchlassbereich, auf die halbe Abtasfrequenz normierte % Grenzfrequenz rp=1; % max. Welligkeit im Duchlassbereich in dB ws=1.5*wp; % Sperrbereich (normierte Frequenz) rs=40; % min. Daempfung im Sperrbereich in dB % Filtersynthese: Cauer−Filter (elliptisch) [n,wn]=ellipord(wp,ws,rp,rs); [b,a]=ellip(n,rp,rs,wn);% Filterkoeffizienten % Filter−Eigenschaften [H,w]=freqz(b,a,2^10); % Frequenzgang eines zeitdiskreten Systems Ha=abs(H); % Betrag Hp=180*unwrap(angle(H))/pi;% 2pi−Spruenge aus dem Phasenverlauf entfernen f=w*fs/2/pi;% Frequenz aus normierter Kreisfrequenz(!) bestimmen figure(1) subplot(2,1,1) loglog(f,Ha), hold on loglog([0.01 wp wp]*fs/2,10.^(−[rp rp rs]/20),'k',... [ws ws 1]*fs/2,10.^(−[rp rs rs]/20),'k'), hold off grid, axis([fg/10 fs/2 0.001 1]) title('Cauer−Tiefpassfilter') xlabel('\rightarrow Frequenz in Hz') ylabel('\rightarrow Amplitudengang') subplot(2,1,2) semilogx(f,Hp), grid axis([fg/10 fs/2 −400 100]) xlabel('\rightarrow Frequenz in Hz') ylabel('\rightarrow Phasengang in Grad')

7

A. Matlab-Hilfe: Digitale Filter (Synthese und Anwendung) A.4.2. Beispiel FIR-BP-Filter (Fenstermethode)

Listing 8: Matlab-Beispiel FIR BP-Filter % FIR−Filter−Spezifikationen fuer BP−Filter fs=2e3; % Abtastfrequenz in Hz Ts=1/fs; % Abtastintervall in s fSB1o=40; % obere Grenzfrequenz des unteren Stoppbands in Hz fPBu=80; % untere Grenzfrequenz des Passbands in Hz (−3dB−Frequenz) fPBo=126; % obere Grenzfrequenz des Passbands in Hz (−3dB−Frequenz) fSB2u=250; % untere Grenzfrequenz des oberen Stoppbands in Hz rSB1=0.01; % zugelassene Welligkeit im Stoppband1 rPB=0.005; % zugelassene Welligkeit im Passband rSB2=0.01; % zugelassene Welligkeit im Stoppband2 fg=[fSB1o fPBu fPBo fSB2u]; amp=[0 1 0]; % Spezifizierte Amplitudengangwerte (Stopp−, Pass−, Stoppband) dev=[rSB1 rPB rSB2]; % Bestimmung der Filterparameter filtpar = kaiserord(fg,amp,dev,fs,'cell') % liefert [n,wn,beta,filtype] b = fir1(filtpar{:}); a=1; % Nennerpolynom = 1 % Filter−Eigenschaften N=2^9; [H,w]=freqz(b,a,N); % Frequenzgang Ha=abs(H); % Amplitudengang Hp=180*unwrap(angle(H))/pi;% Phasengang in Grad und 2pi−Spruenge entfernen f=w*fs/2/pi;% Frequenz aus normierter Kreisfrequenz bestimmen figure(1) subplot(2,1,1) loglog(f,Ha), grid axis([fPBu/10 fs/2 0.0001 1]) title(['FIR−Bandpassfilter, f_{gu} = ',num2str(fPBu),' Hz, f_{go} = ',... num2str(fPBo),' Hz, f_s = ',num2str(1e−3*fs),' kHz']) xlabel('\rightarrow Frequenz in Hz') ylabel('\rightarrow Amplitudengang') subplot(2,1,2) semilogx(f,Hp), grid axis([fPBu/10 fs/2 −1600 400]) xlabel('\rightarrow Frequenz in Hz') ylabel('\rightarrow Phasengang in Grad') figure(2) plot(f,Hp) grid title(['FIR−Bandpassfilter, f_{gu} = ',num2str(fPBu),' Hz, f_{go} = ',... num2str(fPBo),' Hz, f_s = ',num2str(1e−3*fs),' kHz']) xlabel('\rightarrow Frequenz in Hz') ylabel('\rightarrow Phasengang linear in Grad')

8

A. Matlab-Hilfe: Digitale Filter (Synthese und Anwendung) A.4.3. Beispiel FIR-TP-Filter (Parks-McClellan Algorithmus)

Listing 9: Matlab-Beispiel Parks-McClellan Algorithmus % Spezifikation Bandsperrfilter fs=2e3; % Abtastfrequenz Ts=1/fs; % Abtastintervall n=64; % Filterordnung (hier durch Probieren gesucht) fn=[0 200 250 300 350 1000]; % in Hz wn=fn/(fs/2); % normierte Frequenz, einheitslos amp=[1 1 0.004 0.004 1 1];% Amplitudenwerte % Filtersynthese b=firpm(n,wn,amp);% equiripple (Parks−McClellan, Standardverfahren) a=[1]; % Nennerpolynom = 1 % Frequenzgang und Gruppenlaufzeit N=2^9; % Anzahl Stuetzwerte fuer Frequenz [H,w]=freqz(b,a,N); % Frequenzgang Ha=abs(H); % Amplitudengang Hp=180*unwrap(angle(H))/pi; % Phasengang in Grad und 2pi−Spruenge entfernen f=(fs/2)*w/pi; % Frequenz aus normierter Frequenz bestimmen GD=grpdelay(b,a,N); % Gruppenlaufzeit in Anzahl Abtasintervallen figure(1) subplot(3,1,1) plot(f,Ha,fn,m,'−−k') title(['Parks−McClellan optimal equiripple FIR Filter ',num2str(n),... '. Ordnung, f_s = ',num2str(fs*1e−3),' kHz']) xlabel('\rightarrow Frequenz in Hz') ylabel('Amplitudengang') legend('Filter','Spezifikation',4) grid subplot(3,1,2) plot(f,Hp,[0 fs/2],[0 −n*180/2],'−−k') xlabel('\rightarrow Frequenz in Hz') ylabel('Phasengang in Grad') grid subplot(3,1,3) plot(f,1e3*GD*Ts) axis([0 fs/2 0 1e3*n*Ts]) xlabel('\rightarrow Frequenz in Hz') ylabel('Gruppenlaufzeit in ms') grid

9