Digitale Signalverarbeitung mit MATLAB : Grundkurs mit 16 ausführlichen Versuchen 9783834890825, 3834890820 [PDF]

Zitiervorschau

Martin Werner Digitate Signalverarbeitung mit MATLAB

Aus dem Programm Nachrichentechnik / Kommunikationstechnik Telekommunikation von D. Conrads Kommunikationstechnik von M. Meyer Signalverarbeitung von M. Meyer Digitale Kommunikationssysteme 1 von R. Nocker Digitale Kommunikationssysteme 2 von R. Nocker Nachrichtentechnik von M. Werner Nachrichten-Ubertragungstechnik von M. Werner Signale und Systeme von M. Werner information und Codierung von M. Werner Netze, Protokolle, Schnittstellen und Nachrichtenverkehr von M. Werner

vieweg

Martin Werner

Digitale Signalverarbeitung mit MATLAB Grundkurs mit 16 ausfuhrlichen Versuchen 3., voUstandig iiberarbeitete und aktualisierte Auflage Mit 159 Abbildungen und 67 Tabellen

Studium Technik

SQ

Vieweg

Bibliografische Information Der Deutschen Bibliothek Die Deutsche Bibliothek verzeichnet diese Publikation in der Deutschen Nationalbibliografie; detaillierte bibliografische Daten sind im Internet iiber abrufbar.

Das Buch erschien in der 2. Auflage unter gleichem Titel mit Untertitel Intensivkurs mit 16 Versuchen ebenfalls im Vieweg Verlag

1. Auflage Oktober 2001 2., verbesserte und erweiterte Auflage Marz 2003 3., vollstandig iiberarbeitete und aktualisierte Auflage April 2006 Alle Rechte vorbehalten © Friedr. Vieweg & Sohn Verlag | GWV Fachverlage GmbH, Wiesbaden 2006 Lektorat: Reinhard Dapper Der Vieweg Verlag ist ein Unternehmen von Springer Science+Business Media. www.vieweg.de Das Werk einschlieBlich aller seiner Telle ist urheberrechtlich geschiitzt. Jede Verwertung auBerhalb der engen Grenzen des Urheberrechtsgesetzes ist ohne Zustimmung des Verlags unzulassig und strafbar. Das gilt insbesondere fiir Vervielfaltigungen, Ubersetzungen, Mikroverfilmungen und die Einspeicherung und Verarbeitung in elektronischen Systemen. Umschlaggestaltung: Ulrike Weigel, www.CorporateDesignGroup.de Technische Redaktion: Andreas MeiBner, Wiesbaden Druck und buchbinderische Verarbeitung: Wilhelm & Adam, Heusenstamm Gedruckt auf saurefreiem und chlorfrei gebleichtem Papier. Printed in Germany ISBN-10 3-8348-0043-0 ISBN-13 978-3-8348-0043-5

Vorwort Die digitale Signalverarbeitung ist eine allgegenwartige Schlusseltechnologie des Informationszeitalters. Wir tragen sie mit unseren Mobiltelefonen und MP3-Playem herum. Im Auto hilft sie als Anti-Blockiersystem Unfalle vermeiden. In der Klinik ermoglicht sie als Computertomographie Einblicke ohne Operationen. Taglich und uberall macht sie unser Leben bequemer, produktiver und sicherer. Digitale Signalverarbeitung mit MATLAB®^ stellt eine Auswahl wichtiger Grundlagen und einfacher praktischer Anwendungen vor. Anhand von 16 Versuchen werden die Themen behandelt: •

Zeitdiskrete Signale (Versuche 1 und 2)

•

Signalverarbeitung im Frequenzbereich / FFT (Versuche 3 bis 6)

•

Signalverarbeitung im Zeitbereich / LTI-Systeme (Versuch 7 bis 9)

•

Filterentwurf (Versuch 10 und 11)

•

Stochastische Signale (Versuche 12 und 13)

•

Reale Systeme, A/D-Umsetzung, Wortlangeneffekte (Versuche 14 bis 16)

Zu jedem Versuchen gibt es definierte Lemziele und eine kompakte Einfiihrung. Die Versuchsvorbereitung ist ein wesentiiches Element des Praktikums und kann je nach Vorkenntnissen verschieden lange dauem. Die Versuche sind so gestaltet, dass bei guter Vorbereitung, die Bearbeitungszeit am PC ca. 180 Minuten nicht iibersteigen sollte. Es war leider nicht moglich, alle Versuche gleich aufwandig zu gestalten. Viele Aufgaben sind mit Losungshinweisen versehen. Der ausfuhrliche Losungsteil am Ende des Buches ermoglicht ein effektives Erarbeiten des Lemstoffes und Erfolgskontrolle. Zahlreiche Beispielprogramme und Messbeispiele erleichtem die Durchfiihrung der Versuche. Bei der Auswahl der Themen und der Gestaltung der Versuche wurde darauf geachtet, dass keine tiefergehenden Kenntnisse der digitalen Signalverarbeitung erforderlich sind. Allgemeine Kenntnisse iiber Signale und Systeme aus einer einfiihrenden Lehrveranstaltung oder einem einfiihrenden Lehrbuch werden jedoch vorausgesetzt. Die Versuche werden mit der Simulationssprache MATLAB am PC durchgefiihrt. Fiir die Wahl von MATLAB waren folgende Griinde ausschlaggebend: •

MATLAB ist ein haufig benutztes Werkzeug fiir die digitale Signalverarbeitung und wird weltweit auf PC und Arbeitsplatzrechnem mit unterschiedlichen Betriebssystemen eingesetzt.

•

Die einfache Bedienung sowie die guten graphischen Eigenschaften von MATLAB ermoghchen es, von Beginn an die digitale Signalverarbeitung in den Mittelpunkt zu stellen.

•

Die Kombination aus PC und MATLAB bietet eine preiswerte Schnittstelle zu realen Audiosignalen.

MATLAB® ist ein eingetragenes Warenzeichen der Firma The MathWorks, Inc., U.S.A., w w w. math works .com

VI

Vorwort

•

Alle fur das Praktikum erstellten MATLAB-Programme konnen abgeandert und erweitert werden. Die Experimente lassen sich nach personlichen Bediirfnissen und Interessen modifizieren und erweitem.

•

MATLAB wird in vielen Hochschulen eingesetzt und ist auch als preiswerte Studentenversion (Student Edition) erhaltlich.

Zum Buch sind liber 80 Programme und Datensatze entstanden, die auf der Web-Seite des Vieweg Verlages, www.vieweg.de, kostenlos erhaldich sind. Technische Hinweise und Softwarekompatibilitat Um Verwechslungen vorzubeugen werden Schreibweisen und Begriffe verwendet, wie sie in MATLAB gebraucht werden. Dies betrifft besonders den Dezimalpunkt, 0.1 statt 0,1, und die Bereichsangabe fur Laufindizes, z. B. 0:10 statt 0,1,2,...,10. MATLAB-Befehle, Programmvariablen, und Ahnliches werden durch die Schriftart besonders kenntlich gemacht, wie z. B. h e l p oder p l o t . Alle Programme wurden mit MATLAB Version 7.1 (R14) Service Pack 3 getestet. Da aus didaktischen Griinden nur einfache Befehle verwendet wurden, sollten die Programme auch mit der MATLAB Version 6 laufen. MATLAB verfiigt iiber eine umfassende und gut organisierte Online-Hilfe, die Handbiicher vollwertig ersetzt. Aus diesem Grund konnte auch der Einfiihrungsteil des Buches kurz gehalten werden. Flir das Praktikum ist die MATLAB Signal Processing Toolbox erforderlich. Weitere Informationen zu MATLAB sind auf der Homepage der Firma The MathWorks, U.S.A., www.mathworks.com, oder bei der deutschen Niederlassung The MathWorks GmbH, www.mathworks.de, zu finden. Danksagung Geme bedanke ich mich bei den Studierenden am Fachbereich Elektrotechnik und Informationstechnik der Fachhochschule Fulda, die dieses Praktikum mit viel Engagement, hilfreichen Anregungen und konstruktiver Kritik bereichert haben. Mein besonderer Dank gehort Herm Dipl.-Ing. (FH) Bemd Heil, ohne dessen tatkraftige Untersttitzung im NachrichtentechnikLabor dieses Buch so nicht hatte entstehen konnen. Mit seiner Gitarre (guitar.wav) hat er das Praktikum auf angenehmste Weise bereichert. Fiir die freundhche Untersttitzung bedanke ich mich herzlich bei Herm Courtney Esposito von der Firma The MathWorks.

Fulda, im Februar 2006

Martin Werner

VII

Inhaltsverzeichnis 1

Erste Schritte in MATLAB

1

1.1

Einfiihrung

1

1.2

Programmstart und einfache Befehle

1

1.2.1 1.2.2 1.2.3 1.2.4 1.2.5 1.2.6 1.2.7 1.2.8 1.2.9

MATLAB-Bedienoberflache Einfache arithmetische Operationen Konstanten und Variablen MATLAB heIp-Kommando Vektoren und Matrizen Vordefinierte MATLAB-Funktionen und einfache Graphiken Schreiben eines MATLAB-Programms im E d i t o r / Debugger-Window Verkettete Programme und Unterprogramme MATLAB HELP und MATLAB DEMO

1 1 3 3 4 7 9 11 14

2

Zeitdiskrete Signale

16

2.1

Einfiihrung

16

2.2

Elementare zeitdiskrete Signale

16

2.2.1 2.2.2 2.2.3 2.3 2.3.1 2.3.2 2.3.2 2.3.3

3 3.1 3.2

Einfiihrung Vorbereitende Aufgaben Versuchsdurchfiihrung Audiosignale Einfuhrung Beispiel: Synthese eines Audiosignals Vorbereitende Aufgaben Versuchsdurchfiihrung

Diskrete Fourier-Transformation (DFT) Einfiihrung Grundlagen

16 19 20 20 20 21 23 23

26 26 26

3.2.1

Diskrete Fourier-Transformation

26

3.2.2

Eigenschaften der diskreten Fourier-Transformation

30

3.3

Vorbereitende Aufgaben

31

3.4

Versuchsdurchfiihrung

33

4 4.1 4.2 4.3 4.4

Schnelle Fourier-Transformation (FFT) Einfiihrung Komplexitat Radix-2-FFT-Algorithmus Programmierung der DIT-Radix-2-FFT

36 36 36 37 41

4.4.1

Ordnen der Eingangsfolge (Bit-reversal-Adressierung)

41

4.4.2

Signalverarbeitung im Signalflussgraphen

44

4.5

Vorbereitende Aufgaben

47

4.6

Versuchsdurchfiihrung

48

VIII

Inhaltsverzeichnis

5

Kurzzeit-Spektraianalyse: Grundlagen

50

5.1

Einfiihrung

50

5.2

Grundlagen

50

5.2.1 5.2.2 5.2.3 5.2.4 5.2.5 5.2.6 5.2.7

Abtastung Spektrum des zeitdiskreten Signals Fensterung Diskrete Fourier-Transformation Zero-padding Leakage-Effekt Fensterfolgen

51 52 53 55 58 59 60

5.3

Vorbereitende Aufgaben

62

5.4

Versuchsdurchfiihrung

63

6

Kurzzeit-Spektraianalyse: Beispiele

66

6.1

Einfiihrung

66

6.2

Mehrtonsignal

66

6.2.1 6.2.2 6.2.3 6.3 6.3.1 6.3.2 6.3.3

7 7.1 7.2 7.2.1 7.2.2 7.2.3 7.3 7.3.1 7.3.2 7.4 7.4.1 7.4.2 7.4.3 7.4.4 7.4.5

8 8.1 8.2 8.3 8.4

Mehrfrequenzwahlverfahren Vorbereitende Aufgaben Versuchsdurchfiihrung (DTMF-Signal) Audiosignal Einfiihrung Vorbereitende Aufgaben Versuchsdurchfiihrung (Spektrogramm)

Faltung, Differenzengleichung und LTI-Systeme Einfiihrung Faltung

66 67 68 69 69 71 73

75 75 75

Grundlagen Vorbereitende Aufgaben Versuchsdurchfiihrung (Faltung)

75 77 79

Differenzengleichung 1. Ordnung

79

Goertzel-Algorithmus 1. Ordnung Versuchsdurchfiihrung (Goertzel-Algorithmus 1. Ordnung) Lineare zeitinvariante Systeme (LTI-Systeme) Impulsantwort und Frequenzgang von LTI-Systemen Lineare Differenzengleichung und z-Transformation Differenzengleichung 2. Ordnung - Goertzel-Algorithmus 2. Ordnung Vorbereitende Aufgaben Versuchsdurchfiihrung (Goertzel-Algorithmus 2. Ordnung)

Finite-Impulse-Response-Systeme (FIR-Systeme) Einfuhrung Eigenschaften von FIR-Systemen Vorbereitende Aufgaben Versuchsdurchfiihrung

79 81 82 82 84 87 88 88

89 89 89 92 98

Inhaltsverzeichnis 9

D^

Infinite-Impulse-Response-Systeme (IIR-Systeme)

100

9.1

Einffihrung

100

9.2

Einfluss der Pole auf den Frequenzgang

100

9.3

Blockdiagramm

101

9.4 9.5 9.6 9.7

Impulsantwort Partialbruchzerlegung mit MATLAB Vorbereitende Aufgaben Versuchsdurchfiihrung

102 103 105 Ill

10

EntwurfdigitalerFIR-Filter

113

10.1 10.2 10.3 10.3.1 10.3.2 10.4 10.4.1 10.4.2 10.4.3 10.5 10.5.1 10.5.2 10.5.3 10.6 10.6.1 10.6.2

Einfuhrung FIR-Filterstruktur Toleranzschema Entwurfsvorschrift im Frequenzbereich Vorbereitende Aufgaben Fourier-Approximation Fourier-Reihe des Frequenzganges Vorbereitende Aufgaben Versuchsdurchfiihrung Fourier-Approximation mit Fensterung Einfuhrung Vorbereitende Aufgaben Versuchsdurchfiihrung Chebyshev-Approximation Equiripple-Methode Versuchsdurchfiihrung

113 113 114 114 116 117 117 117 118 118 118 119 119 119 119 121

11

EntwurfdigitalerllR-Filter

124

11.1 11.2 11.3 11.3.1 11.3.2 11.3.3 11.3.4 11.3.5 11.3.6 11.4 11.4.1 11.4.2

Einfuhrung IIR-Filter Entwurf eines Butterworth-Tiefpasses Toleranzschema und Filtertyp Zeitkontinuierlicher Butterworth-Tiefpass Dimensionierung des zeitkontinuierlichen Butterworth-Tiefpasses Vorbereitende Aufgaben (Butterworth-Tiefpass) Bilineare Transformation Vorbereitende Aufgaben (Bilineare Transformation) Entwurf digitaler Tiefpasse nach Standardapproximationen analoger Tiefpasse Einfuhrung Versuchsdurchfiihrung

124 124 127 127 127 128 129 131 132 133 133 134

12

KenngroBen stochastischer Signale

137

12.1 12.2

Einfuhrung Grundlagen

137 137

X 12.2.1 12.2.2 12.3 12.3.1 12.3.2 12.4 12.4.1 12.4.2 11.4.3 11.4.4 11.4.5 11.4.6

Inhaltsverzeichnis Experiment und stochastischer Prozess Zufallszahlen am Digitalrechner Stochastische Signale Vorbereitende Aufgaben Versuchsdurchfiihrung Korrelation stochastischer Prozesse Korrelation, Korrelationsfunktion und Leistungsdichtespektrum WeiBes Rauschen Schatzung der Autokorrelationsfunktion Schatzung des Leistungsdichtespektrums Vorbereitende Aufgaben Versuchsdurchfiihrung

13 Stochastische Signale und LTI-Systeme 13.1 Einfuhrung 13.2 Lineare Abbildung stochastischer Signale 13.2.1 13.2.2 13.2.3 13.3 13.3.1 13.3.2 13.3.3

Grundlagen Vorbereitende Aufgaben Versuchsdurchfiihrung Abbildung stochastischer Signale an LTI-Systemen Grundlagen Vorbereitende Aufgaben Versuchsdurchfiihrung

137 140 140 140 142 143 143 144 146 146 147 148

152 152 152 152 154 154 155 155 155 159

14

Analog-Digital-Umsetzung

161

14.1

Einfiihrung

161

14.2

Digitalisierung

161

14.3 Abtastung 14.3.1 Abtasttheorem 14.3.2 Vorbereitende Aufgaben 14.3.3 Versuchsdurchfiihrung

162 162 163 163

14.4

164

14.4.1 14.4.3 14.4.3 14.4.4 14.4.5

15

Quantisierung Quantisierungskennlinie Maschinenzahlen Quantisierungsfehler Vorbereitende Aufgaben Versuchsdurchfiihrung

164 164 168 170 173

Reale digitate Filter: Koeffizientenquantisierung

175

15.1

Einfiihrung

175

15.2

Wortlangeneffekte

175

15.3

FIR-Filter mit quantisierten Koeffizienten

176

Fehlermodell und Fehlerfrequenzgang Vorbereitende Aufgaben Versuchsdurchfiihrung

176 177 178

15.3.1 15.3.2 15.3.3

Inhaltsverzeichnis 15.4 15.4.1 15.4.2 15.4.3 15.4.4

IIR-Filter mit quantisierte Koeffizienten Kaskadenform Koeffizientenquantisierung und Polausdiinnung Vorbereitende Aufgaben Versuchsdurchfiihrung

X£ 179 179 181 183 187

16

Reale digitale Filter: Quantisierte Arithmetik

188

16.1 16.2 16.2.1 16.2.2 16.2.3 16.3 16.4

Einfuhrung Quantisierte Arithmetik Addition: Uberlauf und groBe Grenzzyklen Multiplikation: inneres Gerausch und kleine Grenzzyklen Ersatzschaltbild flir einen Block 2. Ordnung Vorbereitende Aufgaben Versuchsdurchfiihrung

188 188 188 189 193 194 198

17

Losungen zu den Versuchen

205

17.1 17.2 17.4 17.5 17.6 17.7 17.3 17.8 17.9 17.10 17.11 17.12 17.13 17.14 17.15 17.16

Vorbemerkungen Losungen zu Versuch 2: Zeitdiskrete Signale Losungen zu Versuch 3: Diskrete Fourier-Transformation (DFT) Losungen zu Versuch 4: Schnelle Fourier-Transformation (FFT) Losungen zu Versuch 5: Kurzzeit-Spektralanalyse: Grundlagen Losungen zu Versuch 6: Kurzzeit-Spektralanalyse: Beispiele Losungen zu Versuch 7: Faltung, Differenzengleichung und LTI-Systeme Losungen zu Versuch 8: Finite-Impulse-Response-Systeme (FIR-Systeme) Losungen zu Versuch 9: Infinite-Impulse-Response-Systeme (IIR-Systeme) Losungen zu Versuch 10: Entwurf digitaler FIR-Filter Losungen zu Versuch 11: Entwurf digitaler IIR-Filter Losungen zu Versuch 12: KenngroBen stochastischer Signale Losungen zu Versuch 13: Stochastische Signale und LTI-Systeme Losungen zu Versuch 14: Analog-Digital-Umsetzung Losungen zu Versuch 15: Reale digitale Filter: Koeffizientenquantisierung Losungen zu Versuch 16: Reale digitale Filter: Quantisierte Arithmetik

205 205 206 209 212 213 216 219 222 225 229 233 238 245 247 250

Formelzeichen und Abkiirzungen

254

Literaturverzeichnis

257

Sachwortverzeichnis

260

Erste Schritte in MATLAB 1.1

Einfuhrung

Es liegt in der Natur der Sache, dass ein so machtiges Werkzeug wie MATLAB weder auf wenigen Seiten beschrieben noch in den wichtigsten Funktionen schnell beherrscht werden kann. Dieser Abschnitt soil Sie deshalb bei den ersten Schritten in MATLAB unterstiitzen. Anhand einfacher Beispiele wird gezeigt, wie Sie arithmetische Ausdriicke verarbeiten, einfache Graphiken erzeugen, eigene Programme mit Unterprogrammen und Funktionen erstellen und beniitzen konnen. Im Laufe des Praktikums werden sich Ihnen mit zunehmender Ubung die Moglichkeiten von MATLAB erschlieBen. Die konmientierten Programmbeispiele zu den einzelnen Versuchen und die ausfuhrliche Dokumentation von MATLAB werden Ihnen dabei helfen. Anmerkungen: (i) Ausfiihrlichere Einfiihrungen in MATLAB bieten z. B. [ABRW02][GrGr02][GrGr05] [Pra02]. (ii) Uber Biicher zur Signalverarbeitung mit MATLAB informiert die Webseite von The MathWorks mathworks.com/support/books/index.html. Eine Aufstellung deutschsprachiger Bucher kann ausgewahlt werden. (iii) Im Folgenden wird MATLAB 7.1 (14) mit Service Pack 3 unter Windows® verwendet. Da nur grundlegende MATLAB-Befehle eingesetzt werden, soUten die Beispielprogramme auch mit der Version 6 lauffahig sein.

1.2

Programmstart und einfache Befehle

1.2.1

MATLAB-Bedienoberflache

Nach dem Start von MATLAB erscheint wie in Bild 1-1 der voreingestellte MATLABD e s k t o p mit Mentileiste, Schaltknopfen und drei geoffneten Fenstem: C u r r e n t D i r e c t o r y ©, Command H i s t o r y (D und Command Window (D. Das Fenster Workspace ® verbirgt sich unter dem Fenster C u r r e n t D i r e c t o r y . Abhangig von Ihrer Installation (Versionsnummer, Studentenversion und Zusatzpakete wie SIMULINK) konnen Ihre Bildschirmanzeigen von den hier gezeigten Beispielen abweichen. Anmerkungen: (i) Die Voreinstellung erreicht man uber die Mentileiste mit Desktop ^ Desktop Layout ^ Default, (ii) Die Bedienung der Fenster mit ihren Meniileisten, Schaltelementen usw. geschieht in der ubiichen Art. 1.2.2

Einfache arithmetische Operationen

Wir beginnen mit der Voreinstellung in Bild 1-1. Fiir die ersten Schritte im interaktiven Modus benotigen wir nur das Coinmand Window. Der LFbersichtlichkeit halber schlieBen Sie die anderen Fenster links durch klicken mit der Maus auf das Symbol x am jeweiligen oberen Fensterrahmen. Im Command Window erscheint der MATLAB Prompt >>. Dort konnen Befehle direkt eingegeben werden. Im interaktiven Modus fasst MATLAB alle Eingaben als Befehle auf und flihrt sie nach Moglichkeit direkt aus.

1 Erste Schritte in MATLAB

[File Edft D e ^

Qeskfcop ^ N o w

Help Current t^ectcwit; j CiV'rogramme^ATLABZIVwork

v iQ

{?]

I Shortcuts ^ How to Add \B What's New Comma«d Window

1 AS Files-

[ s |1

IfBeType

©

To g e t s t a r t e d ^ s e l e c t HATLAB Help o r Deittos froits t h e Help menu, i

1» ®

1 Current Directory jWbrkspacej

®

@ l^^MtJ

i

::

.

:

. ;

:

:

1

^

: ; ; :

1 OVFf ,;j

Bild 1-1 MATLAB-Bedienoberflache (MATLAB Desktop) Tippen Sie einfach 2+3 ein und driicken Sie dann die Eingabetaste J zur Ubemahme des Befehls durch den Rechner. »

2 +3 J

Sie erhalten als Ergebnis die Antwort (answer) ans = 5 MATLAB verfiigt iiber die ublichen arithmetischen Operatoren und dariiber hinaus iiber weitere arithmetische und logische Operatoren, Vergleichsoperatoren und spezielle Zeichen, die den Umgang mit Vektoren und Matrizen erleichtem. Ein Beispiel ist der Potenzoperator "". Geben Sie »

3^3 J

ein und Sie erhalten ans = 2 7 Kompliziertere arithmetische Ausdriicke konnen mit Hilfe von Klammem definiert werden. Auf die Eingabe »

4*(4-3+2/4)

J

antwortet MATLAB mit ans = 6 Zum Andem einer Eingabe stehen Ihnen die Cursor-Tasten zur Verfiigung. Mit T konnen Sie friihere Eingaben wieder in die Kommandozeile laden und danach bearbeiten. SchHeBen Sie eine Kommandozeile mit dem Semikolon; ab, so wird die Anzeige des Ergebnisses unterdriickt. Als Beispiel betrachten wir die letzte Eingabe. Mit der Pfeiltaste T holen Sie sie wieder hervor. Nun schlieBen Sie das Konmiando mit einem Semikolon ab.

1.2 Programmstart und einfache Befehle

»

4^(4-3+2/4);

J

MATLAB antwortet jetzt nur mit dem Prompt. 1.2.3

Konstanten und Variablen

Namen von Konstanten und Variablen beginnen in MATLAB stets mit einem Buchstaben. Sie konnen bis zu 63 Zeichen (Buchstaben, Ziffem, Unterstrich) enthalten. MATLAB unterscheidet zwischen GroB- und Kleinschreibung. Anmerkung: Altere Versionen von MATLAB unterscheiden nur die ersten 31 Zeichen, s. Getting Started with MATLAB, Version 7, 2004. Das nachfolgende Beispiel verdeutlicht die Anwendung von Variablen »

A = 2;

J

»

a = 4;

J

»

B = A/a;

J

Der Wert einer Variablen wird durch Eingabe des Namens am Bildschirm angezeigt. » B J B = 0.5 MATLAB verfiigt iiber einige vordefinierte Konstanten. Von besonderer Bedeutung ist die (Kreis-) Zahl ;rund die imaginare Einheit / bzw.j. So lassen sich komplexe Zahlen einfach eingeben. »

z = l+2j

J

z = 1.0000 + 2 . 0 0 0 0 i und »

z = 2+3i J

z = 2.0000 + 3.00001 Die Zahl TT- genauer die in MATLAB verwendete Naherung - erhalt man mit >> p i J ans = 3.1416 Die Anzeige am Bildschirm hangt von der Einstellung des Ausgabeformates ab. Der tatsachlich intern verwendete Wert ist davon unabhangig. Das jeweils letzte Ergebnis wird in der Variablen a n s gespeichert. Beachten Sie auch, dass die vordefinierten Konstanten durch Befehlseingabe uberschrieben werden konnen. Anmerkungen: (i) MATLAB benutzt das Long-Format des lEEE-Floating-point-Standard mit einer Genauigkeit von ca. 16 Dezimalstellen und einem Wertebereich von etwa 10'^^^ bis 10"^-^^^. (ii) Fragen der numerischen (Un-) Genauigkeit sind wichtiger Gegenstand der digitalen Signalverarbeitung und werden in den Versuchen 14 bis 16 genauer behandelt.

1.2.4

MATLAB help-Kommando

Tippen Sie einfach h e l p p i ein und drucken Sie dann die Eingabetaste J zur LFbemahme des Befehls durch den Rechner.

1 Erste Schritte in MATLAB

>> h e l p p i

J

Sie erhalten die kurze Erklarung PI

3.1415926535897

PI = 4'^atan(l) = imag (log(-l) ) = 3.1415926535897 In der zweiten Zeile ist der Zusammenhang zwischen der Zahl K und der ArkustangensFunktion a t a n und dem Imaginarteil imag der natiirlichen Logarithmus-Funktion l o g angegeben. Funktionen in MATLAB werden spater noch genauer vorgestellt. Wenn Sie mehr iiber die Formateinstellung wissen wollen, so konnen Sie sich die Optionen des Befehls f o r m a t mit h e l p f o r m a t J anzeigen lassen. 1.2.5

Vektoren und Matrizen

Vektoren und Matrizen sind als geordnete Abfolgen von Zahlen in natiirlicher Weise als zeitdiskrete Signale aufzufassen und spielen in der digitalen Signalverarbeitung eine herausragende Rolle. MATLAB erleichtert den Umgang mit Vektoren und Matrizen durch spezielle Befehle zu ihrer Erzeugung und Verkniipfung. Vektoren lassen sich mit eckigen Klammem „ [ ] " und dem Semikolon „ ; " durch Angabe der Zahlenwerte einfach erzeugen. Geben Sie dazu folgendes Beispiel ein »

X = [1 2 3] J

Sie erhalten von MATLAB den Zeilenvektor X = 1

2

3

Mit Semikolon zwischen den Zahlen ergibt sich der Spaltenvektor »

y = [ 1 ; 2 ; 3]

y =

J

1 2 3

Mit den letzten beiden Befehlen haben Sie in MATLAB Datenfelder erzeugt, sogenannte A r r a y s , deren Elemente im MATLAB Arbeitsspeicher, dem Workspace, abgelegt sind. tJber den Inhalt und die Organisation des Workspace informieren Sie sich im Command Window durch Eingabe des Befehls whos. » whos J Name

Size

A

1x1

8

double array

B

1x1

8

double array

a

1x1

8

double array

X

1x3

24

double array

Y z

3x1

24

double array

1x1

16

double array (complex)

Bytes

Class

Grand total is 10 elements using 88 bytes

1.2 Programmstart und einfache Befehle

Der Befehl whos listet alle im Speicher abgelegten Variablen mit ihren Dimensionen, s i z e genannt, in der Form „Zahl der Zeilen x Zahl der Spalten" auf. Zusatzlich werden Anzahl der verwendeten Bytes und Typ der Variablen angegeben. Anmerkungen: (i) Die Form der Anzeige kann im Workspace-Fenster im Menii View eingestellt werden. (ii) Falls die Variable ans verwendet wurde, taucht sie ebenfalls hier auf. Altemativ konnen Sie in der Meniileiste unter D e s k t o p die Option Workspace auswahlen. Dann erhalten sie das Workspace-Fenster mit der Ubersicht und zusatzliche Bearbeitungsmoglichkeiten, s. Bild 1-2. Die Anzeige konnen sie beispielsweise mit dem Menii-Punkt View und dann Choose Columns anpassen. ^^^^^^HBBtellctJI Vtew

[eie &&

S-apNcs De^pg ^es}> c l e a r

J

>> whos J >> Wiederholen Sie die Eingaben zweier Vektoren x und y mit »

X = [1 2 3] ; J

»

y = [4; 5; 6 ] ; J

Die Dimension der Vektoren kann auch mit dem Befehl s i z e bestimmt werden. >> s i z e ( x ) ans = 1

J 3

Die Angabe ist wieder im Format „Zahl der Zeilen x Zahl der Spalten". Fiir eindimensionale A r r a y s kann auch der Befehl l e n g t h verwendet werden. Da MATLAB eine relativ aufwendige Speicherverwaltung durchfiihrt, konnen bei kompatiblen Dimensionen viele Befehle direkt auf Vektoren oder Matrizen angewandt werden. Wir betrachten einige Beispiele. Zunachst beginnen wir mit einem einzelnen Element eines Vektors. Durch »

x(2) J

1 Erste Schritte in MATLAB ans = 2 erhalten Sie den aktuellen Wert des zweiten Elements des Vektors x. Mit »

x(2)

= 7 J

ans = 1

1

3

wird der zweiten Komponente der Wert 7 zugewiesen. Beachten Sie, dass MATLAB die Komponenten eines Vektors mit dem Index 1 beginnend adressiert! »

x(0)

J

??? Subscript indices must either be real positive integers or logicals. Mit MATLAB konnen auch Vektoren miteinander verkniipft werden. Die Addition der beiden Vektoren >> x+y J liefert eine Fehlermeldung, da die Dimensionen bzgl. der Addition inkompatibel sind. Der Zeilenvektor x kann nicht mit dem Spaltenvektor y addiert werden. ??? Error using ==> + Matrix dimensions must agree. Ein haufiger Programmierfehler ist die Verkniipfung von Datenfeldem mit inkompatiblen Dimensionen - MATLAB wird Sie durch Fehlermeldungen darauf hinweisen. Durch den Operator ' kann eine Transposition (Zeilen- und Spaltenvertauschung) eines Vektors Oder einer Matrix vorgenommen werden. Mit >> z = y"

J

z = 4

5

6

erhalt man einen Zeilenvektor der nun elementweise zu x addiert werden kann. >> x+z J ans = 5

12

9

Die Addition von zwei Spaltenvektoren ist ebenso moglich. Anmerkung: Mit ' werden Vektoren und Matrizen mit komplexen Zahlen zusatzlich konjugiert. Analog kann eine elementweise Multiplikation durchgefiihrt werden. Dazu wird der Multiplikationsoperator mit einem vorangestellten Punkt erweitert . * . »

X.*z J

ans = 4

35

18

Werden hingegen die beiden Vektoren nur mit dem Multiplikationszeichen * verkniipft, so wird - bei kompatiblen Dimensionen - das Skalarprodukt (14 + 7-5 + 3-6) ausgefiihrt. »

x*y J

a n s = 57 Die explizite Definition von Vektoren und Matrizen durch Eingabe der Elemente kann beschwerlich sein. Um dem abzuhelfen, bietet MATLAB spezielle Befehle an. Mit den folgenden Befehlen werden haufig benotigte Matrizen einfach erzeugt.

1.2 Programmstart und einfache Befehle

>> X = ones(2,3) J

X =

1

1

1

1

1

1

>> X = zeros (3,2) J

X =

0

0

0

0

0

0

>> X = repmat(7,2,3) J

X =

7

7

7

7

7

7

Eine weitere wichtige Moglichkeit Datenfelder zu erzeugen, ist die Anwendung des Doppelpunkt-Operators: . Er erzeugt Datenfelder mit Elementen gleichen Abstandes. Durch »

t = 0:10

J

wird ein eindimensionales Datenfeld (Zeilenvektor) mit der Bezeichnung t erstellt und angezeigt, das die Werte von 0 bis 10 in Einser-Schritten enthalt. t = Columns 1 through 7 0

1

2

3

4

5

6

Columns 8 through 11 7

8

9

10

Die vorgestellten Beispiele geben einen kleinen Einblick in die Moglichkeiten mit MATLAB Vektoren und Matrizen zu erzeugen und zu bearbeiten. Im Laufe der Versuche werden sich weitere dazugesellen. 1.2.6

Vordefinierte MATLAB-Funktionen und einfache Graphiken

Eine Starke von MATLAB ist die umfangreiche Sammlung von vordefinierten Funktionen. Dies trifft sowohl auf die MATLAB-Grundausstattung als auch auf die Erweiterungspakete, Toolboxes genannt, zu. Die in den Versuchen benotigten Funktionen werden bei Bedarf eingefiihrt. Exemplarisch soil nun die Sinusfunktion graphisch dargestellt werden. Dazu wahlen Sie zunachst eine Anzahl von aquidistanten Stiitzstellen. »

t = 0:.1:10;

J

Mit >> y = s i n ( t ) ;

J

erstellen Sie ein Datenfeld, das die Ergebnisse der Sinusfunktion angewandt auf jedes Element von t enthalt. Eine graphische Darstellung der Sinusfunktion in einem eigenen Fenster erzeugen Sie durch den Befehl

1 Erste Schritte in MATLAB

>> p l o t { t , y )

J

Die graphische Darstellung kann mit weiteren Befehlen erganzt werden, s. Bild 1-3. >> g r i d

J

>> xlabel('t \rightarrowM J >> ylabel('y(t) \rightarrowM J >> title('Sinus function') J Mit dem Befehl h e l p p l o t erhalten Sie eine Zusammenfassung von Optionen und eine Liste von Befehlen, mit denen Sie spater graphische Ausgaben ihren Bediirfnissen anpassen konnen. Daruber hinaus bietet Ihnen MATLAB im Graphikfenster nachtraglich Optionen zur Bearbeitung von Graphiken an.

^f^^

"^m;l Insert• Ms^ M^i^^^^^^^

am mm

Bild 1-3 Graphische Darstellung der Sinusfunktion mit MATLAB Die Darstellung am Bildschirm geschieht in der Regel nach linearer Interpolation der Funktionswerte zwischen den Stiitzstellen. Fiir die digitale Signalverarbeitung - wenn nur wenige Funktionswerte dargestellt werden sollen - ist die Interpolation irrefiihrend. Deshalb enthalt MATLAB den Graphikbefehl stem. Machen Sie sich den Unterschied deutUch, indem Sie die graphische Darstellung fiir die Sinusfunktion bei einer geringen Stiitzstellenzahl pro Periode mit dem Plotbefehlen p l o t und s t e m wiederholen. » >> >> >>

t = 0:10; J y = sin(t); J plot(t,y) J grid J

>> xlabel('t \rightarrowM J >> ylabel('y(t) \rightarrowM J >> title('Sinus function with few samples') J Das resultierende Bild 1-4 lasst deudich die interpolierten Geradenstiicke zwischen den Stiitzwerten erkennen.

1.2 Programmstart und einfache Befehle

- tap "M 1

D I* e m k • Sl ^ f) ® «' D 13: a H 1

Sinus function wth few samples

1

.

4

1

1

.

.

.

.

J

K

• /'^

0.5

t 1 ^

1

i ""^-'

1 )

1

2

3

4

6

I 6

7

8

9

10

t-»

1

1

Bild 1-4 Sinusfunktion mit dem Graphikbefehl p l o t bei wenigen Stiitzstellen Fiihren Sie nun die graphische Darstellung mit dem Befehl s t e m durch. >> s t e m ( t , y ) J >> g r i d J >> xlabel('t \rightarrowM J >> ylabel('y(t) \rightarrow') J >> title('Sampled sinus function in discrete-time representation') J Das jetzt dargestellte Bild 1-5 hebt den diskreten Charakter der Stiitzstellen hervor.

Rb i i t

View In^rt Tods iesWtcp m^km jHte^j

Di*asul^Af)®|«|i

Bild 1-5 Sinusfunktion in diskreter Darstellung mit dem Graphikbefehl stem

10

1.2.7

1 Erste Schritte in MATLAB

Schreiben eines MATLAB-Programms im E d i t o r / Window

Debugger-

Nachdem Sie das Command Window kennen gelemt haben, sollen Sie sich nun mit einfachen MATLAB-Programmen vertraut machen. Da MATLAB als Interpreter-Sprache die Befehlszeile sequentiell bearbeitet, liegt es nahe, mehrere Eingaben in einer Text-Datei, dem MATLAB S c r i p t F i l e , zusammenzufassen. Damit MATLAB derartige Dateien als solche erkennen kann, sind sie mit der Endung .m zu versehen, wie beispielsweise s i n . m , oder p i .m. Sie werden deshalb kurz als M - F i l e bezeichnet. Um ein M - F i l e zu erstellen, konnen Sie in der Meniileiste des MATLAB D e s k t o p den Mentipunkt F i l e ^ New '^ M - F i l e anwahlen. Danach erscheint das E d i t o r / Debugger-Window am Bildschirm. Sie konnen nun die in Bild 1-6 angegebenen Programmzeilen eingeben: 9 Editor - C:\Dokumente and EinsteMun&en\werner\Ei&ene Dateien\v1ewe. File

gik

lext

Cell Imh

Debug Desktop

WraJow

hjelp

D 0* •':'x>III'e or-/'a '#i'V,''¥i i« * pji'i 1 2 3 4 ~ 5 6

% plot sirms function % rrs^prograrc's.m t « 0:pi/100:5; Y ~ sin{t); piot(t,7)| % end {scrijst

(In S

Col 10

lOVR

Bild 1-6 Programmerstellung im Editor/Debugger-Window Zum Speichem des Programms gehen Sie nun zur Meniileiste, klicken dort F i l e "^ Save an und geben im File-Save-Fenster den Namen des Programms myprogram ein, der dann automatisch mit der Endung . m versehen wird. Giinstigerweise legen Sie zum Speichem von MATLAB-Programmen, Daten usw. ein eigenes Verzeichnis an. Stellen Sie das Arbeitsverzeichnis C u r r e n t D i r e c t o r y von MATLAB im Command Window auf Ihr Verzeichnis um. Uber den Menii-Punkt D e s k t o p "^ C u r r e n t D i r e c t o r y kann das Arbeitsverzeichnisfenster eingeblendet werden. Beachten Sie auch, mit dem Prozentzeichen % definieren Sie den Rest der Zeile als Kommentar. Nutzen Sie die Moglichkeit der Kommentierung ausgiebig, um Ihre Programme verstandlich zu machen. Mit dem help-Kommando und dem Namen eines M-Files, z. B. h e l p myprogram, werden alle dort direkt am Dateianfang stehenden Kommentarzeilen angezeigt. Damit lassen sich auf einfache Weise Hilfstexte zu Ihren Programmen erstellen. Das Programm kann nun uber verschiede Wege gestartet werden. Im Command Window. miissen Sie dazu hinter dem MATLAB-Prompt den Namen des Programms eingeben - aber ohne die Endung . m. Nun wird das Programm abgearbeitet und die Sinusfunktion graphisch dargestellt. Im E d i t o r / D e b u g g e r - W i n d o w wird das Programm iiber das Menii Debug kurz mit der Steuertaste F5, zur Ausfiihrung gebracht.

Run, Oder

1.2 Programmstart und einfache Befehle

1J_

Anmerkung: In MATLAB sind niitzliche Debug-Werkzeuge eingebaut, mit denen Sie sich spater vertraut machen konnen. Eine weitere Moglichkeit ein M-File zu starten ist die Auswahl des Programms mit der rechten Maustaste im Fenster C u r r e n t D i r e c t o r y und Wahl der Option Run. Bei der Interpretation der Befehle priift MATLAB zunachst, ob die eingegebene Zeichenkette eine GroBe im Workspace benennt. Ist das nicht der Fall, wird in den im MATLAB P a t h angegebenen Verzeichnissen nach einem M - F i l e mit entsprechendem Namen gesucht. MATLAB beginnt dabei immer im aktuellen Arbeitsverzeichnis. Das erste vom Namen passende M - F i l e wird ausgefiihrt. 1.2.8

Verkettete Programme und Unterprogramme

Das kleine Programmbeispiel in Bild 1-6 ist vom Prinzip her eine Verkettung von MATLABProgrammen. Taucht in einem M - F i l e der Name eines weiteren M - F i l e auf, so wird es geoffnet und sein Inhalt wie Tastatureingaben im Command Window interpretiert. Alle neu definierten Datenfelder werden im Workspace abgelegt und stehen auch nach der Bearbeitung des M - F i l e als globale Variablen zur Verfugung. Dies ist fiir kleine Programmbeispiele sehr passend. Bei umfangreichen Projekten wiirden sich jedoch sehr schnell folgende Probleme einstellen: •

iJberlastung des Arbeitsspeichers

•

unstrukturierte Programme

•

unzumutbare Programmlaufzeiten

•

haufige Programmierfehler und unbeabsichtigtes Uberschreiben von Daten

MATLAB unterstiitzt deshalb MATLAB-eigene und anwenderdefinierte Funktionen. Bei den MATLAB-eigenen Funktionen handelt es sich um M - F i l e s oder in sich geschlossene speicher- und laufzeitoptimierte Programm-Module. Anwenderdefinierte Funktionen werden zwar prinzipiell wie Tastatureingaben im Command Window interpretiert, besitzen aber eine definierte Schnittstelle zum aufrufenden Programm und verwalten ihre Daten lokal. Das folgende Beispiel soil Ihnen das Arbeiten mit selbstdefinierten Funktionen erlautem. Wir wahlen ein etwas anspruchsvolleres Beispiel aus der Signalverarbeitung, die Approximation eines periodischen Rechteckimpulszuges durch den Gleichanteil und die Harmonischen der Fourier-Reihe. Bild 1-7 zeigt den periodischen Rechteckimpulszug in normierter Darstellung. Die Periode betragt 1. Die Amplitude altemiert im Abstand der halben Periode zwischen -i-l und 0. x(t)

1 0

1/2

1

2

3

Bild 1-7 Periodischer Rechteckimpulszug in normierter Darstellung

22

1 Erste Schritte in MATLAB

Wir erstellen zunachst eine MATLAB-Funktion die die Funktionswerte des periodischen Rechteckimpulses fiir eine graphische Darstellung liefert. Bevor Sie beginnen loschen Sie zuerst den Workspace mit dem Kommando c l e a r im Command Window. Dann offnen Sie ein neues E d i t o r / D e b u g g e r - W i n d o w und geben die folgenden Zeilen ein. Programmbeispiel 1-1 Rechteckimpulszug function y = rectangular(t,p,w) % y = rectangular(t,p,w) % t : time samples (t>=0) % p : period % w : impulse width % rectangular.m * mw * 10/05/2005

y = zeros(size(t));

% default

amplitude

0

for n = 1:length(t)

X = mod(t(n),p)); % mapping of t into

fundamental

period

if x>=0 & x= und for haben Sie zugriff auf die jeweilige MATLAB-Dokumentation. Die Approximation des periodischen Rechteckimpulszuges geschieht mit seiner aus der Mathematik bekannten Fourier-Reihe 12°° x(0 = - + — y

1

sin((2n + l)-2;^-/)

Dazu erstellen Sie folgendes Progranmi: Programmbeispiel 1-2 Fourier-Reihe function yF = fourier(t,N) % fourier series expansion of periodic rectangular impulse train % yF = fourier(t,N) % t : time samples % N : number of harmonics used (N>=1) % fourier.m * mw * 10/5/2005 yF = 0.5*ones(size(t)); for n=0:N-l yF = yF + {(2/pi)/(2*n+l))*sin((2*n+l)*2*pi*t); end

Speichem Sie die Funktion als M - F i l e mit dem Namen f o u r i e r . m ab.

1.2 Programmstart und einfache Befehle

13

AbschlieBend wird das Hauptprogramm zur Festlegung der Parameter, zum Aufrufen der Funktionen und der graphischen Ausgabe erstellt. Speichem Sie das folgende Hauptprogramm als M - F i l e mit dem Namen f o u r i e r s y n . m ab. Prograimnbeispiel 1-3 Hauptprogramm % Fourier synthesis % fouriersyn.m * mw * 10/05/2005 t = 0:.01:3; % time samples

y = rectangular(t,1,.5); % rectangular synthesis yF = fourier(t,10); % fourier % graphics

impulse with 10

train harmonics

plot(t,y,t,yF) axis([0 3 -0.2 1.2]); grid xlabel('t \rightarrow') ylabeK'y(t), y_{F}(t) \rightarrow') title('Fourier synthesis of rectangular impulse train')

Nach dem Aufruf des Programms f o u r i e r s y n sollten Sie die Graphik in Bild 1-8 erhalten. Sie zeigt das aus der Mathematik bekannte Ergebnis mit den tJberschwingem der abgebrochenen Fourier-Reihe, dem gibbschen Phanomen. Letzteres spielt in der digitalen Signalverarbeitung eine groBe Rolle und wird in einem spateren Versuch noch genauer behandelt.

Bte m

§§w Insert - 1 ^ ^

WM^^

W^'

D * O a \k I ^ ^ f):® 1 € 1:1 :@1 •

Bild 1-8

Periodischer Rechteckimpulszug in normierter Darstellung y(t) und seine Approximation yp(t) durch den Gleichanteil und die ersten zehn Harmonischen

In der MATLAB-Dokumentation wird zwischen einem S c r i p t und einer F u n c t i o n unterschieden. Ersteres erzeugt als Folge von MATLAB-Befehlen globale Variablen, wahrend eine Funktion mit lokalen Variablen arbeitet und Daten iiber ihre Schnittstelle, dem Funktionsaufruf, austauscht. Mit dem Befehl whos oder im Fenster Workspace konnen Sie das an den Programmbeispielen liberpriifen. Die in den Funktionen verwendeten Variablen x und n sind als lokale Variablen nicht im Arbeitsspeicher enthalten.

14

1.2.9

1 Erste Schritte in MATLAB

MATLAB H e l p und MATLAB Demo

Eine wichtige Moglichkeit sich mit MATLAB vertrauter zu machen bietet das MATLABFenster Help, s. Bild 1-9. Dort finden Sie unter dem Link G e t t i n g S t a r t e d eine Einfuhrung in MATLAB. Das Fenster ist u. a. liber das Menii Help, das Fragezeichen ? oder durch anklicken von MATLAB h e l p im Command Window erreichbar. Falls das fiir Sie der erste Kontakt mit MATLAB sein sollte, bedenken Sie, dass ein Programmsystem wie MATLAB nicht in einem Trockenkurs anhand einer Bedienungsanleitung erlemt werden kann. Hier ist „Lemen durch gezieltes Probieren" der sinnvollste Weg. Die folgenden Versuche werden Ihnen dabei Schritt fiir Schritt helfen. MATLAB bietet mit dem Kommando demo eine Auswahl von animierten Einfiihrungen und Anwendungsbeispielen an, s. Bild 1-10. Als Vorbereitung auf die weiteren Versuch sollten Sie sich mit den einfUhrenden Video-Beispielen zur Entwicklungsumgebung D e s k t o p T o o l s and Development E n v i r o n m e n t vertraut machen, s. Bild 1-11.

File idffe view go Favorites Debug Deskltop S^lndow Hefp

Qm Shortcuts IE How to Add [S Whal'sNew ^^ ^ ' K/

Conterifis j index:- Se^ch =i Demos ^ Q Begin Here *i 0 Release Notes 3t^^lnstallatJon •B # MATLAB 9 -0 Communications Toolbox m '0 image Acquisition Toolbox II 0 Image Processing Toolbox S 0 Signal Processing Toolbox ^ Support and Web Services

Sr

V%

Tiki: \ Release 14 with Service Pack 3: Begin Here

Begin Here Reltase 14 with Service Pack 3 If You Are Upgrading from a Previous Release... m Release Notes Highlights nev^ features, installation notes, bug fixes, and compatibility issues.

If You Are Using MATLAB for the First Time... At the heart of MATLAB is a new language that you must learn before you can fully exploit its power. This isnt as hard as it might sound; you can learn the basics of MATLAB very quickly. You v^ill be rev^arded with high-productivity, high-creativity computing power that will change the way you work. If you are a first-time user, the best way to get started is to read thoroughly the Getting Started with MATLAB tutorial with MATLAB open so you can follow along. The tutorial book comes with MATLAB and is available in PDF and for purchase on the MathWorks Web site. If vnit Hnnl w»r>t tn t»kp thp timp tn rpaH it thnrnttrihlw hprp sirp linkXf+ 2 ^W^^A^

n=0,2,...

(4.4)

n=lX--

Die Substitutionen n = 2m

fiir n gerade

n = 2m + l fiir n ungerade und M = N/2 liefem dann

M 5^

38

4 Schnelle Fourier-Transformation

M-l

M-l

X[k] = 2 ] x[2m]wjj^^ + 2 ] 42m + l]w^^'""^^^^ m=0

(4.6)

m=0

Beriicksichtigt man noch die Umformungen fiir die komplexen Faktoren

4'"'=
#r-V%ictow Mle?s^fer-Me!::ctoiyini:f

Frecpency doff^to

:E38:/::: :; « : ; l^eslcige Factor:; 0 i 3 : % ; ; ;

B^m:m^£&WQm^m^:-MS.£di:;:;•; :: ::: Myiob|#ilii|«3dli|::i0.039063

- Window IM-" . Seleci ;v%?ihdQV o o o o o o o o o ! o o o o o o o 6 o A -10

8

10

12

14

16

18

m

Bild 7-4 Faltung der Barker-Codefolge b[n] mit sich selbst nach Zeitumkehr (tr, time reversal)

7.3 Differenzengleichung 1. Ordnung

79

7.2.3

Versuchsdurchfiihrung (Faltung)

D

Kontrollieren Sie das Ergebnis in Aufgabe 7.2-1 mit dem MATLAB-Befehl conv.

n

In der Informationstechnik wird haufig die Faltung von zwei Rechteckimpulsen betrachtet. Fallen Sie mit MATLAB den Rechteckimpulse der Lange 3, x[n] = {1,1,1}, mit sich selbst. Skizzieren Sie das Ergebnis und kontrollieren Sie es durch eigene Rechnung. Fatten Sie mit MATLAB die Rechteckimpulse der Lange 3 und 5, x[n] = {1, 1, 1} und y[n] = {1, 1, 1, 1, 1}. Skizzieren Sie das Ergebnis und Kontrollieren Sie es durch eigene Rechnung.

n

Fuhren sie die Faltung im Programmbeispiel 7-1 auch mit der Barker-Codefolge der Lange 13 Z7M = {1,1, 1,1,1,-1,-1, 1,1,-1,1,-1,1} durch und stellen Sie das Ergebnis graphisch dar. Beachten Sie die GroBenverhaltnisse zwischen dem Maximum und den anderen Werten. Weshalb sind Barker-Codefolgen ftir Synchronisationsaufgaben interessant? Erklaren Sie anhand der Resultate die Besonderheit von Barker-Codefolgen.

7.3

Differenzengleichung 1. Ordnung

7.3.1

Goertzel-Algorithmus 1. Ordnung

Als einftihrendes Beispiel fiir die Anwendung einer Differenzengleichung betrachten wir den Goertzel-Algorithmus, der z. B. in Telefonie zur Erkennung der Tone beim Mehrfrequenzwahlverfahren eingesetzt wird, vgl. Versuch 6. Die Aufgabe, das gesendete Tonpaar beim Mehrfrequenzwahlverfahren zu erkennen, besteht im Grunde darin, die Signalenergien bzgl. alter acht zugelassenen Frequenzen zu bestimmen und zu vergleichen. Eine einfache praktische Naherung liefert der Goertzel-Algorithmus mit der effizienten Berechnung der jeweils spektral nachstliegenden DFT-Koeffizienten (3.1). X[k]=^x[n]'W^N^

(7.7)

n=0

Die Herleitung des Goertzl-Algorithmus geschieht in zwei geschickt gewahlten Schritten. Im ersten wird die Berechnung des DFT-Koeffizienten als Faltungssumme dargestellt. Weil Wj^ = 1, konnen wir dazu (7.7) umformen

X[k]=^x[n]'w/^'^ ""^

(7.8)

n=0

Der Vergleich mit der Faltungssumme (7.1) zeigt prinzipielle Ubereinstimmung dann, wenn im Argument der Exponentialfunktion der Term N-n ^LIS Differenz der Laufindizes des Faltungsergebnisses m und der Faltungssumme n verstanden werden kann.

80

7 Faltung, Differenzengleichung und LTI-Systeme

Als Ergebnis der Faltung des urspriinglichen Signals x[n] mit dem Hilfssignal h[n] ergibt sich das Hilfssignal y[n] n

y[n] = 2 4 / ] • w^^^*^""^^ = x[n] * h[n]

(7.9)

mit M«] = [ w ^ ^ 7

(7.10)

Und mit x[N] = 0 gilt die Ubereinstimmung mit dem gesuchten DFT-Koeffizienten X[k] = y[N]

(7.11)

Im zweiten Schritt wird die Berechnung der Faltungssumme als eine Realisierung des Prinzips der Rtickkopplung (Rekursion) eingeftihrt, wie es auch in den Differenzengleichungen zum Ausdruck kommt. Anmerkungen: Ohne das fundamentale Prinzip der Rtickkopplung ware unsere Welt nicht denkbar. Als fachubergreifende Wissenschaft beschaftigt sich die Kybemetik (N. Wiener, 1948) damit. Andersherum gesehen, ist ebenso zu erwarten, dass sich fur Differenzengleichungen (und Differentialgleichungen) viele Anwendungsfelder ergeben. Vielfaltige technische Anwendungen werden in der Regelungstechnik mit dem Begriff Regelkreis zusammengefasst. Schreibt man die Faltungssumme (7.9) fiir die Werte von y[n] explizit an y[0] = MO]-40] = 40] y[l] = h[0] • 41] + hm' 40] = 41] + h[l] • y[0] y[2] = h[0] • 42] + h[l] • 41] + h[2] • 40] = x[2] + h[l] • (41] + h[l] • 40]) =

(7.12)

= x[2] + h[l]'ym so ergibt sich aus einem Koeffizientenvergleich wegen der Produktdarstellung (7.10) der rekursive Zusammenhang y[n] = x[n]-\-h[l]'y[n-l]

furn>0

und3;[-l] = 0

(7.13)

Es liegt eine lineare Differenzengleichung mit konstanten Koeffizienten (DGL) 1. Ordnung vor y[n] + ai-y[n-l]

= bo'x[n]

(7.14)

mit den Koeffizienten ^1 ~ "^N

u^^ ^0 ~ 1

(7.15)

und den Nebenbedingungen x[n] = 0 yl" 1] = 0

fiir

n^{OX...,N-l} ^ ^ Anfangswert

(7.16)

Die rekursive Struktur der DGL 1. Ordnung wird nochmals im Programmbeispiel 7-2 und im Blockdiagramm des zugehorigen Systems in der Direktform I deutlich sichtbar.

7.3 Differenzengleichung 1. Ordnung

81_

Das Blockdiagramm in Bild 7-5 zeigt die Berechnung des Ausgangswertes y[n] in Abhangigkeit vom aktuellen Eingangswert x[n] und dem vorherigen Ausgangswert j[n-l]. Der Block @ symbolisieren eine Verzogerung (Delay) des eingehenden Signals um einen Takt. Die kreisformigen Elemente stehen flir die Rechenoperationen „Multiplikation des Signals mit der angegebenen Konstanten" bzw. „Addition aller ankommenden Signale". Verzweigen sich Kanten, so werden die Signale entsprechend dupliziert. Werden die Nebenbedingungen (7.16) eingehalten, so kann im Zeitschritt A^ der gesuchte DFTKoeffizient am Systemausgang abgegriffen werden. Zur Berechnung des DFT-Koeffizienten sind mit dem vorgestellten Goertzel-Algorithmus A^ komplexe Multiplikationen und bei den iiblichen reellen Eingangssignalen N reelle Additionen erforderlich. Pro DFT-Koeffizient ergibt sich dadurch keine Erspamis zur direkten Berechnung der DFT-Koeffizienten. Sind jedoch nur wenige Koeffizienten zu bestimmen, z. B. beim Mehrfrequenzwahlverfahren acht von 256, ist die Einsparung insgesamt erheblich. Hinzu kommt, dass nun anders als bei der Radix-2-FFT, die DFT-Lange prinzipiell nicht gebunden ist. FUr den praktischen Einsatz kann der Algorithmus noch effizienter gestaltet werden, so dass nur reelle Operationen durchgefiihrt werden miissen, was im Abschnitt 7.4.3 noch genauer erlautert wird. Eingang

x[n] o

Ausgang

\-r)

1

^ ^^^^

System Bild 7-5 Blockdiagramm eines zeitdiskreten Systems 1. Ordnung in Direktform I fiir den GoertzelAlgorithmus 1. Ordnung Programmbeispiel 7-2 DGL 1. Ordnung ftir den Goertzel-Algorithmus (Programmausschnitt) al = -exp(j*2*pi*k/N) ; Y = x(l); for n = 2:N; y = x(n)-al*y; end

% h[l] % y[0] = x[0] % recursion

% DFT

Xk = -al*y;

coeffizient

7.3.2

Versuchsdurchfiihrung (Goertzel-Algorithmus 1. Ordnung)

D

tiberpriifen Sie den Goertzel-Algorithmus. Erzeugen Sie dazu ein Sinus/Kosinussignal dessen /:-fache Periode 256 ergibt. Wenden Sie den Goertzel-Algorithmus fiir verschiedene DFT-Koeffizienten an und iiberpnifen Sie die Ergebnisse. Schreiben Sie dazu eine MATLAB-Funktion fiir den Goertzel-Algorithmus als System 1. Ordnung function y = goertzell(x,k,N)

82

7 Faltung, Differenzengleichung und LTI-Systeme

D

Erzeugen Sie ein DMFT-Signal der Lange 205 und bestimmen Sie mit dem Goertzel-Algorithmus die den acht DMFT-Tonen am nachsten liegenden DFT-Koeffizienten. Bild 7-6 zeigt das Ergebnis fiir den Wahlton 6. Anmerkung: Die DFT-Lange 205 liefert die besten Resultate fiir die DMFT-Signalerkennung [Mar92][Mit98].

dsplab?_4 : DMFT analysis with Goertzel algorithm Fte

edft

*^sm Insert

lools

Desktop

Window

P]M)B3

Help

1

j

i

m

X ^ ^

40

o

20

>??! IS

M 9 9 < ^QQQQi 20

26

[>QOQ9^

m

>?I

3S

U 40

Bild 7-6 Mit dem Goertzel-Algorithmus 1. Ordnung berechnete DFT-Koeffizienten {N = 205) fiir das DMFT-Signal zum Wahlton 6

7.4

Lineare zeitinvariante Systeme

7.4.1

Impulsantwort und Frequenzgang von LTI-Systemen

Die Systemtheorie definiert Signale als mathematische Funktionen. Systeme sind demgemaB Abbildungen oder Transformationen von Funktionen, s. Bild 7-7 links. Der mathematische Ansatz wird in der digitalen Signalverarbeitung erganzt durch Uberlegungen zur praktischen Realisierung in Hard- und Software. Ihrer Bedeutung entsprechend, werden den digitalen Systemen insgesamt sechs Versuche gewidmet. Zeitdiskrete linearen zeitinvarianten Systeme, kurz LTI-Systeme genannt (Linear Time-Invariant), werden insbesondere durch ihre Reaktionen auf die Erregung mit der Impulsfolge ^n\ und der komplex Exponentiellen e^ ^" . Erregt man ein zeitdiskretes LTI-System mit einer Impulsfolge wie in Bild 7-7 rechts, so ist am Ausgang die Impulsantworten h[n\ zu beobachten. Hangt die Form der Impulsantwort nicht vom Zeitpunkt der Erregung ab, spricht man von einem zeitinvarianten System. Fiir die Berechnung der Ausgangssignale ist wichtig, dass die Signale der digitalen Signalverarbeitung in natiirlicher Weise als Uberlagerung gewichteter und zeitlich verschobener Impulsfolgen aufgefasst werden konnen. Interpretiert man das Eingangssignal dementsprechend und setzt die Linearitat des Systems voraus, d. h. die Anwendbarkeit des Superpositionsprinzips, dann resultieren die Signale am Systemausgang als Uberlagerung ebenso gewichteter und verzogerter Wiederholungen der Impulsantwort. Fiir die Eingangs-Ausgangsgleichung zeitdiskre-

83

7.4 Lineare zeitinvariante Systeme

ter LTI-Systeme ergibt sich demnach im Zeitbereich die Faltung der Eingangsfolge mit der Impulsantwort. y{n\ = x[n] * h[n]

Eingang x[n\ Erregung

System

T{.}

(7.17)

Ausgang • y[n\ Reaktion

fjl^^j^oy^J^on

Bild 7-7 Zeitdiskretes System als Transformation einer Folge (links) und Reaktion zeitdiskreter LTISysteme (rechts) auf die Impulsfunktion und komplex Exponentielle mit der Impulsantwort bzw. dem Frequenzgang Die zweite wichtige Art der Systemreaktion ist die Reaktion auf komplex Exponentielle. Dahinter verbirgt ist das Konzept der Eigenfunktion. Also, dass sich am Systemausgang bis auf einen multiplikativen Faktor das Signal am Systemeingang ergibt [Wer05]. Die komplex Exponentiellen sind Eigenfunktionen von LTI-Systemen. Wird also ein LTI-System mit einer komplex Exponentiellen beliebiger Frequenz erregt, also mit 4^] = ^^ ^^ ^ so reagiert das LTI-System mit y[n] = H(eJ^^ )'e

jQ.Qn

(7.18)

Der konstante, komplexe Faktor H(e^ ^) gibt das tJbertragungsverhalten des Systems fiir die komplex Exponentielle mit der normierten Kreisfrequenz QQ wider. Anmerkung: Dies ist auch der Grund dafiir, dass am Ausgang von LTI-Systemen nur Frequenzanteile auftreten, die auch am Systemeingang eingespeist wurden. Obige Uberlegungen hangen nicht vom Wert der normierten Kreisfrequenz ab. Es gilt daher allgemein: Liegt fiir das Einganssignal eine harmonische Zerlegung in Frequenzkomponenten durch die Fourier-Transformation vor, so kann das Ausgangssignal im Frequenzbereich einfach bestimmt werden

Y(eJ^) = H(eJ^)'X(eJ^)

(7.19)

J^ ) dem Frequenzgang des Systems. mit H(e^ Da sowohl die Impulsantwort als auch der Frequenzgang die Systemreaktion eindeutig beschreiben, miissen beide in engem Zusammenhang stehen. Tatsachlich bilden die Impulsantwort und der Frequenzgang ein Fourier-Paar und konnen, die Stabilitat des Systems vorausgesetzt, durch Fourier-Transformation ineinander umgerechnet werden. r

h[n]^H(eJ^)

(7.20)

Im Weiteren wird, falls nicht anders erwahnt, von kausalen Systemen ausgegangen. Das sind Systeme deren Reaktion erst mit oder nach der Erregung eintritt. Fiir diese ist die Impulsantwort eine rechtsseitige Folge, also h[n\ = 0 fiir n < 0.

84

7 Faltung, Differenzengleichung und LTI-Systeme

Femer nehmen wir an, dass die Systeme auf Eingangssignale mit beschrankten Amplituden mit Ausgangssignalen mit beschrankten Amplituden antworten. Man spricht dann von BIBOStabilitdt (Bounded Input Bounded Output). 7.4.2

Lineare Differenzengleichung und z-Transformation

In der Regel konnen die in der digitalen Signalverarbeitung verwendeten LTI-Systeme durch lineare Dijferenzengleichungen mit konstanten Koeffizienten (DGL) beschrieben werden. N

2

M

aj,y[n-k]=

k=Qi

2 ^/4«-/]

(7.21)

1=0

Damit ist insbesondere die Realisierbarkeit der Systeme durch Hard- oder Software prinzipiell sichergestellt, vgl. Blockdiagramm in Direktform I in Bild 7-5. Zur effektiven Losung der DGL stellt die Mathematik die z-Transformation der Signale zur Verftigung + 00

X(z)=

2

xM-z'""

(7.22)

«=-oo

Anmerkung: Die z-Transformation ist eng mit der Fourier-Transformation fiir Folgen verwandt wie ein Blick in Tabelle 7-1 zeigt. Fiir zeitkontinuierliche Funktionen gelten entsprechende Zusammenhange zwischen den linearen Differentialgleichungen und der Lapalce-Transformation. Durch z-Transformation leitet sich aus der Differenzengleichung (7.21) die Eingangs-Ausgangsgleichung im Bildbereich ab. Y{z)=

H(z)'X(z)

(7.23)

Darin ist H(z) die Ubertragungsfunktion des Systems eine rationale Funktion

k=0

k=\

mit den Zdhlerkoejfizienten bi und den Nennerkoejfizienten aj^. Das Zahler- und Nennerpolynom kann auch in Produktform mit den Nullstellen zo i und Polen Zook geschrieben werden. Das Ubertragungsverhalten der Systeme wird durch die Pole und Nullstellen bis auf eine multiplikative Konstante vollstandig charakterisiert. In den folgenden Versuchen 8 und 9 wird deshalb der Einfluss der Pole und Nullstellen genauer betrachtet. Bei stabilen Systemen liefert die Ubertragungsfunktion auf dem Einheitskreis z = eJ^ den Frequenzgang

//(^)^) = ^ . - H

"' na-zoo..-^'^) k=\

(7,25)

7.4 Lineare zeitinvariante Systeme

85

Wichtige KenngroBen zeitdiskreter LTI-Systeme sind in Tabelle 7-1 im Sinne einer Formelsammlung zusammengestellt. Die KenngroBen und ihre Zusammenhange sollen Ihnen in diesem und den nachsten Versuchen durch beispielhafte Anwendungen verstandlich werden.

Tabelle 7-1 Eigenschaften zeitdiskreter LTI-Systeme [Wer05] h[n] = T{S[n]}

Impulsantwort Eingangs-Ausgangsgleichung im Zeitbereich (Faltung)

(7.26)

y[n] = x[n] * h[n] = 00

00

= 2

x[k]'h[n-k]=

k=-oo

2

h[k]'x[n-k]

(7.27)

k=-oo

5[w] = T{w[n]}

Sprungantwort Zusammenhang zwischen Impulsantwort und Sprungantwort

(7.28) (7.29)

s[n]= 5 ] h[k] k=-oo

h[n] = s[n] — s[n — l]

BIBO-Stabilitat (Bounded Input Bounded Output)

Eigenfunktion und Eigenwert (z-Transformation)

(7.30)

2 |Mn]| auf den Betragsfrequenzgang

9.3

Blockdiagramm

Wir veranschaulichen die rekursive Struktur der IIR-Systeme anhand des Blockdiagramms eines Systems 2. Ordnung in Bild 9-3. Darin sind die Verzogerungsglieder durch y , wie Delay, die Multiplizierer durch die jeweiligen Faktoren und die Addierer durch „-i-" kenntlich gemacht. Anmerkung: In der Literatur werden unterschiedliche Symbole verwendet. Beispielsweise das Symbol des Verstarkers (Dreieck) fiir die Multiplikationen oder T fiir die Verschiebung um einen Zeittakt. Zusatzlich eingetragen sind die ZustandsgroBen si[n\ und S2[n\. Sie beschreiben den inneren Zustand des Systems und spielen bei der Zustandsraumdarstellung der Systeme eine wichtige Rolle [Wer05]. Bild 9-3 zeigt eine mogliche Realisierung, die transponierte Direktform 11. Sie zeichnet sich dadurch die jeweils minimale Anzahl von Verzogerungsgliedem und Multiplizierem aus, die fiir ein rekursives System bestimmter Ordnung moglich sind.

9 Infinite-Impulse-Response-Systeme

102 x[n] O

Bild 9-3

Blockdiagramm eines rekursives System 2. Ordnung in transponierter Direktform II (aQ= 1) mit den ZustandsgroBen s\ [n] und 52 [n]

Anmerkungen: (i) Je nach Anwendung werden auch alternative Strukturen eingesetzt. Deren Diskussion wiirde jedoch den hier abgesteckten Rahmen sprengen. (ii) Die Signal Processing Toolbox von MATLAB unterstiitzt auch die Lattice-Struktur mit dem Befehl l a t c f i l t . Im Versuch benutzen wir die ZustandsgroBen um ein einfaches Programm fiir die Realisierung des Systems abzuleiten. Aus Bild 9-3 erhalten wir fiir die AusgangsgroBe y[n] = Si[n]-\-bQx[n]

(9.1)

und mit dem Verschiebungsoperator D D{x[n]} = x[n-l]

(9.2)

fiir die ZustandsgroBen Si [n] = D {52 [n] + bix[n] — a^ y[n]\

(9.3)

und S2 [n] = D {^2-^[^] ~ ^2 yM}

(9.4)

Die Gleichungen werden in der Aufgabe A9.6-2 als Progranmiiervorschrift fiir ein rekursives Systems 2. Ordnung verwendet.

9.4

Impulsantwort

Die Ubertragungseigenschaften der LTI-Systeme charakterisieren im Bildbereich die Ubertragungsfunktion H{z) und im Zeitbereich die Impulsantwort h[n\, s. Tabelle 7-1. Impulsantwort und Ubertragungsfunktion bilden ein z-Transformationspaar, so dass beide Systemfunktionen ineinander iiberfiihrt werden konnen. Die dafiir benotigten Rechenschritte werden durch spezielle MATLAB-Befehle unterstiitzt. Wir machen uns zunachst anhand des Systems 2. Ordnung in Bild 9-3 mit dem grundsatzlichen Verfahren vertraut. Die tJbertragungsfunktion eines rekursiven Systems 2. Ordnung ist eine rationale Funktion mit einem Zahler- und einem Nennerpolynom. Zahler- und Nennerpolynom konnen mit ihren Nullstellen, den NuUstellen zo bzw. den Polen Zoo, auch als Produkte geschrieben werden.

9.5 Partialbruchzerlegung mit MATLAB

103^

bQ+biz~^-^b2Z~^

Hiz) =

zi

bQZ^+b^z + b2

b^

(Z-ZOI)(Z-ZQ2)

17 =

2 " 7 Vr ^ (^-^^ dQ+diZ -^ci2Z a^z +aiz + a2 % U ^xiJU"^oo2>' Fiir den wichtigen Sonderfall nur einfacher, von null verschiedener Pole kann die tJbertragungsfunktion durch Partialbruchzerlegung in die fiir die inverse z-Transformation giinstige Form iiberfiihrt werden. H(z) = BQ + _^^ + J^ Z

Zool

Z

fiir

z^oi ^ ^oo2

und Zooi.Zooi ^ 0

(9.6)

Zoo2

Tabellen fiir z-Transformationspaare, z. B. in [Wer05], lief em die Korrespondenz z

^ a'^'u[n]fl

(9.7)

z-a also h[n] = BoS[n] + (B^ • z^i + ^2 * ^^2) * ^M

(9.8)

mit der Impulsfolge S[n] und den rechtsseitigen Exponentiellen Biz^iu[n] undB2Z^2^M • Sind die Betrage der Pole kleiner eins, so klingen die Exponentiellen mit wachsender normierten Zeit n ab und das System ist stabil. Bei stabilen, reellwertigen Systemen mit konjugiert komplexen Polpaaren ergeben sich bedampfte sinusformige Anteile, die Eigenschwingungen des Systems.

9.5

Partialbruchzerlegung mit MATLAB

Die Berechnung der Impulsantwort aus der Ubertragungsfunktion wird vorteilhaft durch eine Partialbruchzerlegung der Ubertragungsfunktion vorbereitet. Die Partialbruchzerlegung rationaler Funktionen unterstiitzt MATLAB mit dem Befehl r e s i d u e z . Anhand zweier Beispiele sehen wir, wie der Befehl eingesetzt wird: Beispiel 9-1 System 2, Ordnung mit konjugiert komplexem Polpaar

l-1.4z"^+0.74z"^ Die Zahlerkoeffizienten und die Nennerkoeffizienten werden im Command Window als Vektoren eingegeben. »

b = [1 2 1] ; J

»

a = [1 - 1 . 4 0 . 7 4 ] ;

J

Die Befehlseingabe >> [ r , p , k ]

= residuez(b,a);

J

liefert die Partialbruchentwicklungskoeffizienten (residues) r = -0.1757 - 3.64591

104

9 Infinite-Impulse-Response-Systeme -0.1757 + 3.64591

zu den Polen (pols) p =

0.7000 + 0.50001 0.7000 -

0.50001

Der direkte Term (Durchgriff) ist 1.3514

k =

Mit r/ = r (1) und /?/ = p (1) sind alle Koeffizienten fiir die Partialbruchentwicklung der Ubertragungsfunktion bekannt. n I-PI'Z

n ' Z h ' z

vo ^

1-P2'Z

^

Z-Pl

(9.10)

Z-P2

Es resultiert .r . . . . . . . (-0.1757-73.6459)z (-0.1757+ 73.6459)z —+— (9.11) H.(z) = l35l4 + ^' ' z - ( 0 . 7 + 70.5) z-(0.7-70.5) ^ ^ Beachten Sie, dass bei reellwertigen Systemen die Partialbruchentwicklungskoeffizienten konjugiert komplexer Pole stets selbst zueinander konjugiert komplex sind. Dadurch lassen sich gegebenenfalls Fehleingaben oder numerische Ungenauigkeiten erkennen. Beispiel 9-2 System 2. Ordnung mit doppeltem reellen Pol

rj (,.-

1+^"'

l-z

^+0.25z ^

Die Zahlerkoeffizienten und die Nennerkoeffizienten sind »

b

- [1 0 1 ] ;

>> a =

[1

J

-1 0.25];

J

Der Befehl [ r , i :) , k ]

»

= reslduez(b,a);

J

liefert wieder die Partialbruchentwicklungskoeffizienten r

= - 8 . 0000 5 . 0000

zu den Polen p =

0. 5000 0. 5000

Der direkte Term ist k =

4

(9.12)

9.6 Vorbereitende Aufgaben

105

Die Ubertragungsfunktion resultiert somit in H2iz) = k+-

1

^-Pi'z-'

,

'2

.=4+^ ? ^ +- ''

(l-p,^z-'f

'-^-^

(9.13)

(^-0-5)'

Man beachte, dass hier pi = /?2- Unter Umstanden konnen sich bei der numerischen Berechnung der Pole merkliche Ungenauigkeiten einstellen. Bei der hier verwendeten Bildschirmausgabe im MATLAB-Format s h o r t werden die Werte gerundet dargestellt.

9.6

Vorbereitende Aufgaben

A9.6-1

Gegeben ist das rekursives System 2. Ordnung

1 + z"

(9.14)

H(z) = l-0.8z"^+0.64z"

Berechnen Sie in Tabelle 9-1 die ersten vier Werte der Impulsantwort des Systems in (9.14) anhand des Blockdiagrams in Bild 9-3. Gehen Sie dabei vom energiefreien Zustand aus, d. h. ^i[0] = 52[0] = 0^ und verfolgen Sie Takt fiir Takt die SignalgroBen im Blockdiagramm.

Tabelle 9-1 Impulsantwort h[n] des Systems in (9.14) n

0

x[n]

S2[n]

s\[n]

0

0

h[n]

1 2 3

A9.6-2

Programmbeispiel 9-1 i i r d f 2 t .m realisiert ein IIR-System 2. Ordnung in Direktform II entsprechend dem Blockdiagramm in Bild 9-3. Machen Sie sich mit dem Funktionsaufruf vertraut. Warum werden im Funktionsaufruf die ZustandsgroBen verwendet? Erganzen Sie die eigentliche Filterroutine. Kontrollieren Sie Ihr Programm anhand Ihrer Ergebnisse in Tabelle 9-1. Hinweis: Sechs Befehlszeilen sind ausreichend.

9 Infinite-Impulse-Response-Systeme

106

Programmbeispiel 9-1 IIR- Filter 2. Ordnung in der transponierten Direktform II function [y,sf] = iirdf2t(b,a,x,si) % implementation of a 2nd order IIR system in % transposed direct form II % function [y,xf] = iirdf2t(b,a,x.,xi) % X : input signal b = [bO bl, b2] b : numerator coefficients % % a : denominator coefficients a = [aO al, a2] si = [si s2] % si : initial conditions % y : output signal % sf : final conditions % iirdf2t.m ^ mw * 01/06/2006 if nargin==4 SI = si; % initial conditions else SI = [0,0]; end a = a/a(l); % normalize coefficients to aO = 1 b = b/a(l); y = zeros(size(x)); % allocate memory for the output signal % filtering

A9.6-3

Bestimmen Sie die Pole und Nullstellen des Systems H(z) aus (9.14) und tragen Sie sie in Bild 9-4 ein. ^01,2 =

^ool,2 =

9.6 Vorbereitende Aufgaben

107

>Re

Bild 9-4 Pol-Nullstellendiagramm zu H{z) aus (9.14)

A9.6-4

Berechnen Sie die Impulsantwort h[n] von H(z) aus (9.14) und kontrollieren Sie Ihr Ergebnis mil Tabelle 9-1. ^ Fiihren Sie bitte die Rechnung (ca. 1 Seite) auf einem eignen Blatt durch und tragen Sie unten nur die Ergebnisse ein.

A9.6-5

Im Programmbeispiel 9-2 soil die Impulsantwort eines rekursiven Systems berechnet werden. Erganzen Sie den Programmausschnitt so, dass die Impulsantwort fiir ein System mit nur einfachen Polen bestimmt wird. Hinweis: Sechs Befehlszeilen sind ausreichend

108

9 Infinite-Impulse-Response-Systeme

Programmbeispiel 9-2 Berechnung der Impulsantwort fiir nur einfache Pole [r,p,k] = residuez(b,a); n = 0:20; h = zeros(size(n));

% compute impulse

A9.6-6

% allocate

% partial fraction % normalized time memory for the impulse

expansion variable response

response

Obiges Programmbeispiel soil nun zur Bestimmung der Sprungantwort verwendet werden. Wie muss das Nennerpolynom a modifiziert werden, damit das Programmbeispiel 9-2 die Sprungantwort liefert? Hinweis: Gehen Sie von (7.38) aus und betrachten Sie die Ubertragungsfunktion als Funktion vonz"!

A9.6-7

Im Progranmibeispiel 9-3 i i r p l o t f d. m wird ein Programm zur graphischen Darstellung der SystemkenngroBen im Frequenzbereich vorgestellt, s. Bild 9-5. Machen Sie sich mit dem Programmbeispiel 9-3 soweit vertraut, dass Sie das Programm in der Versuchsdurchfuhrung verwenden konnen. Anmerkung: Als Beispiel wird das System „Cauer" aus der Versuchsdurchfuhrung verwendet. Bild 9-6 zeigt die zugehorige Impulsantwort und Sprungantwort. Die Besonderheiten des Systems werden ebenso wie das in der Versuchsdurchfuhrung verwendete System „Butterworth" in Versuch 11 naher vorgestellt.

9.6 Vorbereitende Aufgaben

109

Programmbeispiel 9-3 Graphische Darstellung der Systemfunktionen in Bild- und Frequenzbereich (keine mehrfachen Pole) function iirplotfd(b,a,N,txt) % pole-zero plot and plot of magnitude and phase of % frequency response and group delay of IIR systems % with non-multiple poles % function iirplotfd(b,a,N,txt) % b : numerator coefficients % a : denominator coefficients % N : number of samples in the frequency domain (default N=1024) % txt : text for figure and display titles % iirplotfd.m * mw * 01/06/2006 if nargin*»Jl-#-*'***V

t

m /

Ss

zulassiger Frequenzgang des Betrags

" ^ ^ ' • • ^ • - • ^ • ^ ^ ' ' • • • • • • ' ' ^

Qn

Qc

Q

K

Bild 10-2 Toleranzschema fiir den Entwurf eines FIR-Tiefpasses nach Betrag mit dem Durchlassbereich ®, dem Ubergangsbereich 0 und dem Sperrbereich (D Bei FIR-Filtem ist eine „linearphasige" Realisierung moglich und in der Regel gewiinscht. In diesem Fall kann das Toleranzschema direkt auf den Frequenzgang bezogen werden, s. Bild 10-3. Der Frequenzgang ist dann rein reell bzw. rein imaginar. Die nachfolgend behandelten Verfahren liefem FIR-Filter mit linearer Phase im Durchlassbereich. Um ein gutes Sperrverhalten zu erreichen, liegen die Nullstellen der Filter im Sperrbereich auf dem Einheitskreis. Wie in den Abschnitten 8 und 9 festgestellt wurde, fiihrt eine Nullstelle auf dem Einheitskreis zu einem Phasensprung um K, was in Bild 10-3 einem Vorzeichenwechsel des Frequenzganges entspricht. Das FIR-Filter ist deshalb streng genonmien nicht linearphasig. Genauer spricht man darum von einem FIR-Filter mit verallgemeinerter linearer Phase (generalized linear phase FIR filter). Im Weitem wird jedoch wie iiblich vereinfachend von linearer Phase gesprochen. Beim Entwurf linearphasiger FIR-Filter sind gewisse Randbedingungen fiir den Frequenzgang zu beachten. Je nach Lange der Impulsantwort ergeben sich bei gerader und ungerader Symmetric der Impulsantwort die Einschrankungen in Tabelle 10-1. Anmerkung'. Das sieht man jeweils am einfachsten an einer Skizze einer moglichen Impulsantwort.

116

10 Entwurf digitaler FIR-Filter

zulassiger Frequenzgang

Bild 10-3 Toleranzschema fur den Entwurf eines FIR-Tiefpasses mit linearer Phase Tabelle 10-1 Eigenschaften der Frequenzgange von linearphasigen FIR-Filtem Filterordnung A^ Impulsantwort gerade: h[n] = h[N-n] ungerade: h[n] = -h[N-n]

10.3.2

gerade

ungerade

keine Einschrankung

H(1) = 0

H(0) = 0, H(1) = 0

HiO) = 0

Vorbereitende Aufgaben

A.10.3-1 Benutzen Sie die Ergebnissen in Tabelle 10-1, um die jeweils nicht realisierbaren Filtertypen, d. h. Tiefpass, Bandpass, Hochpass oder Bandsperre, in Tabelle 10-2 einzutragen. Tabelle 10-2 Einschr ankungen fiir linearpllasige FIR-Filter Symmetrie der Impulsantwort h[n]

Ordnung der Impulsantwort N

nicht realisierbare Filtertypen

gerade

keine Einschrankung

gerade ungerade gerade ungerade ungerade

A.10.3-2 Geben Sie die (nomiierte) Durchlasskreisfrequenz Q/), die (normierte) Sperrkreisfrequenz Q^-, die Durchlasstoleranz ^ und die Sperrtoleranz Ss in Tabelle 10-3 so an, dass fiir die digitale Verarbeitung eines analogen Signals die angegebenen Anforderungen erfiillt werden. Zur Umrechnung der Frequenzangaben benutzen Sie den Zusammenhang gemaB dem Abtasttheorem mit der Abtastfrequenz/^, s. a. Versuch 14. Q = 2;r-///,

fur

0 Transformation

l-\-as z = :I-as

Bild 11-7 Bilineare Transformation der 5-Ebene in die z-Ebene und umgekehrt Der Zusammenhang zwischen den Frequenzvariablen erschlieBt sich durch Einsetzen von s =jo)bzw. z = eJ^. Man erhalt nach kurzer Zwischenrechnung a) = — 'im— a 2

und Q = 23.rct2in(a co) f i i r a > 0 ^ ^

(11.13)

Es ist offensichtlich, dass die Abbildung der sich von -oo bis +oo erstreckenden imaginaren Achse der ^-Ebene auf den Einheitskreis der z-Ebene nicht linear geschehen kann. Bei der Abbildung von der ^-Ebene in die z-Ebene ergibt sich im Frequenzgang die so genannte Arcustangens-Verzerrung. Sie muss beim Filterentwurf berticksichtigt werden. Anmerkung: Eine Erweiterung der Idee fiihrt auf die Transformationen von Tiefpassen in Bandpasse, Hochpasse und Bandsperren. Diese Transformationen gehoren ebenfalls zu den Standardmethoden fiir den Entwurf digitaler Filter und werden in MATLAB implizit verwendet. 11.3.6

Vorbereitende Aufgaben

All.3-6

Entwerfen Sie durch bilineare Transformation des Butterworth-Tiefpasses nach Al 1.3-3 einen Tiefpass im Zeitdiskreten. Die normierte 3dB-Grenzkreisfrequenz soil dabei QsdB = ^/ 5 sein. Hinweis: Verwenden Sie fiir die Pole direkt die normierten dimensionslosen GroBen s^]^ aus Al 1.3-4. Tabelle 11-5 Pole des zeitdiskreten Butterworth-Tiefpasses Pole des zeitdiskreten Systems

Transformationsparameter Zool =

a=

Zoo2 =

Zoo3 =

11.4 Entwurf digitaler Tiefpasse nach Standardapproximationen analoger Tiefpasse All.3-7

133

Skizzieren Sie das Pol-Nullstellendiagramm des zeitdiskreten Systems aus Al 1.3-6. Beachten Sie dabei, dass bei der bilinearen Transformation, die Nullstellen fiir 5-> oo auf z = -1 abgebildet werden. Erganzen Sie deshalb fiir jeden Pol eine Nullstelle. Im{z}, >kreis 1 y^

L

^^^--.

1y 1

1R

\

1 \\

1 ^^^i

y

^-.. 1 1 ;::

^X''"" \

Bild 11-8 Pol-Nullstellendiagramm des Tiefpasses im Zeitdiskreten All.3-8

All.3-9

Auf welcher Ortskurve liegen die Pole in der z-Ebene? Tragen Sie die Ortskurve in Bild 11-8 ein. Hinweis: s. konforme Abbildungen , z. B. in [BSMM97], S. 627f. Geben Sie die Ubertragungsfunktion Hjp{z) des in Al 1.3-6 bis 8 entworfenen Tiefpasses an. Es sei Hjp{\) = 1.

11.4

Entwurf digitaler Tiefpasse nach Standardapproximationen analoger Tiefpasse

11.4.1

Einfuhrung

Chebyshev- und Cauer-Tiefpdsse werden ahnlich wie digitale Butterworth-Tiefpasse entworfen. In Abschnitt 11.3 wurden die einzelnen Schritte am Beispiel eines digitalen ButterworthTiefpasses erlautert. Entsprechend werden die Entwurfe fiir Chebyshev- und Cauer-Tiefpdsse durchgefiihrt. An die Stelle des Potenzverhaltens in der Leistungstibertragungsfunktion (11.1) treten nun Chebyshev-Polynome bzw. jacobische elliptische Funktionen. Das prinzipielle Losungsverfahren andert sich nicht, jedoch sind die Losungen jeweils charakteristisch fiir die unterlegten Funktionen. Eine mathematische Behandlung des Approximationsproblems wtirde den hier abgesteckten Rahmen sprengen, weshalb auf die weiterfiihrende Literatur verwiesen wird. Im Versuch beschranken wir uns darauf, den Filterentwurf mit dem MATLAB-Programm f da t o o l durchzufiihren und die Losungen fur die verschiedenen Filtertypen zu vergleichen.

134 11.4.2

11 Entwurf digitaler IIR-Filter Versuchsdurchfuhrung

Im Versuch werden Tiefpasse nach den Standardapproximationen entworfen und miteinander verglichen. n

Uberpriifen Sie Ihren Filterentwurf in Al 1.3-9 mit dem MATLAB Filter Viewer fvtool. Hinweis: Mit der rechten Maustaste im Bild des Filter Viewer offnet sich ein Menu mit der Auswahl A n a l y s i s P a r a m e t e r s ...und S a m p l i n g F r e q u e n c y .... Darunterhaben Sie eine Reihe von Einstellmoglichkeiten, um die graphische Darstellung anzupassen. So kann die lineare Darstellung des Betragsfrequenzganges oder die Darstellung der Frequenzachse fiir die normierte Kreisfrequenz gewahlt werden.

D

Entwerfen Sie mit Hilfe des MATLAB-Programmes in Programmbeispiel 11-1 d s p l a b l l _ l einen Butterworth-Tiefpass entsprechend Al 1.3-2. Die Abtastfrequenz sei 20 kHz. Schatzen Sie die benotigte Filterordnung ab. Stellen Sie fiir den digitalen Tiefpass das Pol-NuUstellendiagramm, den Frequenzgang des Betrages, der Dampfung, der Phase und der Gruppenlaufzeit graphisch dar. Tragen Sie die Filterordnung fiir den Tiefpass nach Al 1.3-2 in Tabelle 11-6 ein. Hinweis: Da das Toleranzschema fur IIR-Tiefpasse in Bild 11-2 sich in den Toleranzen etwas von dem Schema fur FIR-Tiefpasse in Bild 10-2 unterscheidet, verwenden Sie aus Grunden der Vergleichbarkeit die Toleranzen Sp = 0.095 und L

(12.14)

Anderenfalls ergeben sich Fehler ahnlich der Fensterung in der Kurzzeit-Spektralanalyse. Dariiber hinaus kann das LDS durch Zero-padding in einem engeren Frequenzraster dargestellt werden. Programmbeispiel 12-1 enthalt eine hierzu geeignete MATLAB-Kommandozeile. Man beachte die spezielle Anordnung der Daten ftir den f f t-Befehl. Den Ausgangspunkt bilden die L Schatzwerte der AKF Rxx fiir 1=0 : L - 1 . Die symmetrisch fortgesetzten Werte der AKF fiir I = - L + l : - 1 werden zeitUch gespiegelt angehangt. Das Zero-padding geschieht durch Einfiigen von NuUen. Man beachte auch, das LDS eines reellen Prozesses ist stets rein reell und nicht negativ da es die mittlere Leistung in jedem Frequenzteilband widergibt. Damit lasst sich die Giite des Messergebnisses beurteilen. Gegebenenfalls sind, wo vertretbar, kleine Abweichungen zu korrigieren. Programmbeispiel 12-1 Berechnung des LDS aus den Schatzwerten der AKF % compute pds from acf estimates Rxx with zero-padding to length 1024 L = length(Rxx]); Sxx = fft([Rxx zeros(l,1024-{2*L-l)) Rxx(L:-1:2)]);

12.4 Korrelation stochastischer Prozesse

147

Altemativ zu obigem Weg, konnen Schatzwerte ftir das LDS mit Hilfe der FFT direkt aus den Musterfolgen bestimmt werden. Ein wegen seiner Einfachheit oft verwendetes Verfahren beruht auf dem Periodogramm [OSW99]. Dabei wird das LDS durch die Betragsquadrate der DFT-Koeffizienten eines Signalblockes geschatzt. Aufgrund der stochastischen Natur des Signals konnen benachbarte Messwerte stark schwanken. Zur Glattung des Ergebnisses wird in der Regel sowohl eine Fensterfunktion, haufig das Hamming-Fenster, als auch eine Mittelung iiber mehrere Signalblocke verwendet. Das Verfahren ist in Bild 12-4 illustriert. Aus der Musterfolge x[n] werden mit einer Fensterfolge w[n\ Blocke der Lange A^' herausgeschnitten. Fiir jeden Block wird die DFT berechnet. Die Betragsquadrate der DFT-Koeffizienten werden liber die Blocke gemittelt.

Signalblocki

Die Signalblocke iiberlagem sich dabei mittig um ^. n-+N-\ bei relativ kurzen Musterfolgen geniigend Blocke fiir die Mittelung zur Verfugung zu haben. Ist die BJI^ I2-4 Zerlegung einer Musterfolge mit Musterfolge sehr lang, so kann auf das LFberlappen der Fensterfolge w[n\ in uberder Blocke verzichtet werden. Da dann keine lappende Blocke der Lange A^ Signalwerte doppelt verwendet werden, nimmt in diesem Fall die statistische Zuverlassigkeit des Ergebnisses zu. Fiir eine effiziente Berechnung wird eine Blocklange verwendet, die eine Radix-2-FFT zulasst.

12.4.5

Vorbereitende Aufgaben

A12.4-1 Eine Messwerterfassung hat zu den Merkmalen X und Y (stochastischen Variablen) die Messreihen (Beobachtungen) x\, xi, .... x^ und y\, y2,-", yN ergeben. Anhand der beiden Messreihen soil die Korrelation zwischen X und Y durch den empirischen Korrelationskoeffizienten beurteilt werden. Geben sie die Rechenvorschrift fiir den empirischen Korrelationskoeffizienten r^y an. ^ A12.4-2 MATLAB stellt zwei Befehle fur die Berechung der empirischen Kovarianz bzw. des empirischen Korrelationskoeffizienten zur Verfugung. Fiir die Vektoren x und y der Lange N liefert der Befehl C =

COV(X,Y);

als Ergebnis eine Matrizen zuruck mit C(l,l) C(l,2) C(2,l) C(2,2)

= = = =

sum( (x-mean(x) ) .-^2) ) / (N-l) ; sum((x-mean(x)).*(y-mean(y)))/(N-1); C(l,2); sum((y-mean(y))."2))/(N-l);

Der Befehl R = corrcoef(x,y)

liefert die Matrix mit R(i,j)

=

C(i,j)/sqrt(C(i,i)*C(j,j));

Beachten Sie, die Befehle cov und c o r r c o e f konnen mit verschiedenen Optionen ausgefuhrt werden.

12 KenngroBen stochastischer Signale

148 12.4.6

Versuchsdurchfuhrung

n

Der Bedeutung der Korrelation wird zunachst anhand von Streudiagrammen veranschaulicht. Erzeugen Sie die vier Signale (Zufallszahlenfolgen) xi[n] bis X4[n] mit xl =: randn(500,1) ;

x2 = randn(500,1) + xl;

x3 = randn(500,1);

x4 = randn(500,1) - xl;

Stellen Sie die vier Streudiagramme graphisch dar, indem Sie x l fur die Abszissenwerte und x l , x2, x3 bzw. x4 fur die Ordinatenwerte wahlen, wie z. B. beim Graphikbefehl p l o t ( x l , x3 , ' . ' ) . Vergleichen Sie die Streudiagramme. Welche Aussagen konnen gemacht werden? Was folgt daraus fiir die Wertepaare der Signale? D

Bestimmen Sie mit dem MATLAB-Befehl c o r r c o e f die empirischen Korrelationskoeffizienten der Signale xi[n] bis X4[n] bzgl. Signal xi[n]. Tragen Sie die Resultate in Tabelle 12-5 ein und diskutieren Sie die Ergebnisse mit Blick auf die Streudiagramme. Tabelle 12-5 Korrelationskoeffizienten zu den Signalen jci[n] bis x^ln] Signal /

Pu Erzeugen Sie eine wesentlich groBere Stichprobe fiir die Signale xi[n] bis X4[n] und stellen Sie die bidimensionalen WDFen entsprechend den obigen Untersuchungen graphisch dar, s. Bild 12-5. Diskutieren Sie die Ergebnisse im Zusammenhang mit obigen Untersuchungen. Welchen Einfluss hat der Korrelationskoeffizient auf die bidimensionale WDF der Normalverteilung? Hinweis: Verwenden Sie der Einfachheit halber das Programm dsplabl2_3 mit der Funktion h i s t 2 d fiir das bidimensionale Histogramm. Programmbeispiel 12-2 Schatzung der bidimensionalen WDF % estimation of the bidimensional probability density function (pdf) % dsplabl2_3.m * mw * 01/20/2006

% number of samples % number of bins for (1-dim.) pdfs 1 or 2 % minimum value of bin centers % maximum value of bin centers

N = le6; Ml=30; M2 = 30; MINl=-4; MIN2=-4; MAX1= 4; MAX2= 4;

h2d = [M1,MIN1,MAX1,M2,MIN2,MAX2]; % parameter

hist2d % 2-dim.

histogram

(absolute

vector

for

function

frequency)

xl = randn(N,1); x2 = randn(N,l) + xl; [cl,c2,f2d] = hist2d(xl,x2,h2d);

% normalization

(pdf)

Dl = (MAXl-MINl)/(Ml-1); D2 = (MAX2-MIN2)/(M2-1);

% width of bins for xl % width of bins for x2

12.4 Korrelation stochastischer Prozesse f2d = f2d/(Dl'^D2*N) ;

%

149 % normalization

(relative

frequency)

graphics

FIG=figure('Name','dsplabl2_3 : bidim. pdf','NumberTitle 'off,'Units','normal','Position',[.4 .4 .45 .45]); grid, surfl(cl,c2,f2d) % pdf xlabel('x \rightarrow'), ylabel('\leftarrow y') zlabel('f(x,y) \rightarrow') FIG=figure('Name','dsplabl2_3 : contour lines of bidim. pdf',. 'NumberTitle', 'off','Units' 'normal','Position',[.4 . 37 .45 .45] V = [.01 .04 .08 .12 .14 0.16 [CS,CH] = contour(cl,c2,f2d,V ; grid % contour lines clabel(CS,CH,V); xlabelC'X \rightarrow'), ylabel('y \rightarrow')

Programmbeispiel 12-3 Bidimensionales Histogramm function [cl,c2,f2d] = ]iist2d(x,y, ]i2d) % bidimensional histogram with equidistant bin center spacing % function [cl,c2,f2d] = hist2d(x,y,h2d) % X : input sequence 1 % y : input sequence 2 % if y is scalar, y is used as shift parameter, i.e. % histogram of pairs x[n] and x[n+y] is computed % h2d : parameters for bin settings % h2d = [M1,MIN1,MAX1,M2,MIN2,MAX2] % Ml, M2 : number of equally sized bins % MINI, MIN2 : minimum value of bin centers MAXl, MAX2 : maximum value of bin centers % bin centers with respect to xl % cl % c2 bin centers with respect to x2 % f2d histogram (absolute frequency) % hist2d.m '"^ mw * 01/20/2006 Lx length(x); Ly = length(y); Dl (h2d(3)-h2d(2))/(h2d(l)-l); % width of bins of x D2 (h2d(6)-h2d(5))/(h2d(4)-l); % width of bins of y h2d(2)+Dl*(0:h2d(l)-l); cl h2d(5)+D2*(0:h2d(4)-l); c2 h2d(2)-Dl/2; Bl B2 h2d(5)-D2/2 if Ly==l % one sequence x and y=x[n+m] m = y; y = x(l+m:Lx); X = x(l:Lx-m); Ly = length(y); end L = min(Lx,Ly); x = x(l:L); y = y(l:L) X = ceil((x-Bl)/Dl) ; X = min(x,h2d(l) ) X = max(X,1) y = ceil({y-B2)/D2) • y = min(y,h2d(4) ) y = max(y,1) f2d = zeros(h2d(4),h2d(l)); % histogram for n=l:L f2d(y(n),x(n)) = f2d(y(n),x(n)) + 1; end

150

12 KenngroUen stochastischer Signale

dspJab12__3 : b i d i m . p d f

\-3^m111

View In-serfe l o o \s l^skfcop Jincksw HMP

t h? on 1 "IS 0,05.

0. 5 4 -^-"-""^ •2 S

i

0

4

dsplabi 2_3 : contour lines of bidim. pdf Ffe

Idte

^ w

Insert

lools

geskfeop Window

Help

4f

Bild 12-5 Graphische Darstellung der bidimensionalen WDF einer Normalverteilung (oben) und ihrer Hohenlinien (unten) mit dem Programm dsplabl2_3 n

Testen Sie, ob der MATLAB-Rauschgenerator r a n d n unkorrelierte Zufallszahlen erzeugt. Fiillen Sie dazu Tabelle 12-6 aus. Ftir die Uberpriifung verwenden Sie den MATLAB-Befehl aus der Signal Processing Toolbox x c o r r beispielsweise in der Form [R,l] = x c o r r ( r a n d n ( l e 4 , 1 ) , 2 0 , ' u n b i a s e d ' ) ;

151

12.4 Korrelation stochastischer Prozesse Tabelle 12-6

Schatzwerte fiir die AKF RxxU] der durch den MATLAB-Befehl randn erzeugten Zufallszahlenfolge fiir den Stichprobenumfang A^

N

max(l/?xx[/]l)fur/=l,2, ...,20

Rxxm

10^ 10^ 10*

n

Horen Sic sich eine Probe des von r a n d n erzeugten Rauschsignals z. B. mit dem Befehl s o u n d s c an. Variieren Sie dabei die Grenzfrequenz (Abtastfrequenz) des Rauschsignals.

D

Bestimmen Sie die AKF und das LDS des Sprachsignals in der Datei speech.wav. Welche Grundfrequenz hebt sich im LDS deutlich hervor? Hat dies Auswirkungen aufdieAKF?

tjber welchen Zeitraum liegt eine starke Korrelation (Koharenz) vor. Als Schatzwert fiir die Koharenzdauer verwenden Sie die Zeitverschiebung bis die AKF 50% ihres Maximalwertes erstmalig unterschreitet.

Vergleiche Sie die Ergebnisse auch mit den Resultaten in Versuch 6, Abschnitt 6.3.

H Hinweise zu MATLAB-Funktionen und M-Files Im Folgenden werden fiir die Versuchsdurchfiihrung ntitzliche MATLAB-Befehle und -Funktionen aufgelistet, zu denen Sie sich mit Hilfe der Help-Funktion Erlauterungen und Beispiele am Bildschirm anzeigen lassen konnen, s. a. vorherige Versuche. Tabelle 12-7 MATLAB-Befehle - benutzte Programme und Dateien, s. a. vorhergehende Versuche Graphikbefehle

clable,

mathematische Funktionen

c e i l , c o r r c o e f , cov, erf, h i s t c , mean, s t d , v a r

Funktionen der Signalverarbeitung

rand,

vorbereitete Programme und Dateien

contour,

randn,

dsp l a b l 2 _ l . m , dsp labl2_4.m,

surfl,

zlabel erfinv,

hist,

xcorr

dsplabl2_2.m, dsplabl2_3.m, histogram.m, hist2d.m, speech.wav

152

13

Stochastische Signale und LTI-Systeme

13.1

Einfiihrung

Dieser Versuch knlipft an den Untersuchungen zu den statistischen KenngroBen stochastischer Signale an. Im Mittelpunkt steht die Frage nach den Veranderungen der KenngroBen bei Filterung der Signale durch LTI-Systeme. Lemziele: Nach Bearbeiten dieses Versuches konnen Sie •

den linearen Mittelwert und die Varianz einer stochastischen Variablen durch eine iineare Abbildung einstellen

•

den linearen Mittelwert und die Varianz am Ausgang eines LTI-Systems berechnen

•

Die Zeit-Autokorrelationsfunktion und die Leistungsfibertragungsfunktion eines LTI-Systems berechnen

•

die Autokorrelationsfunktion und das Leistungsdichtespektrum am Ausgang eines LTI-Systems angeben

•

die Bedeutung des Korrelationskoeffizienten anhand der zweidimensionalen WDF einer Normalverteilung erlautem

13.2

Lineare Abbildung stochastischer Signale

13.2.1

Grundlagen

Stochastische Signale treten in der digitalen Signalverarbeitung an vielen Stellen auf: •

Informationstragende Signale sind von zufalligem Charakter, ansonsten waren sie bekannt und brauchten nicht iibertragen oder verarbeitet zu werden.

•

Informationstragende Signale treten aufgrund der physikalischen Rahmenbedingungen haufig zusammen mit Storsignalen auf, wie z. B. thermisches Rauschen in Verstarkem.

•

Reale digitale Systeme arbeiten mit endlicher Wortlange, so dass nach einer Multiplikation das so genannte Rundungsrauschen auftritt.

Die Verarbeitung stochastischer Signale erzeugt neue stochastische Signale mit neuen Eigenschaften. Man spricht von Abbildungen stochastischer Variablen und Prozesse. In der Informationstechnik ist es wichtig, derartige Abbildungen zu verstehen, um beispielsweise gezielt storunterdriickende MaBnahmen anwenden zu konnen. Um den vorgesehenen Rahmen nicht zu sprengen, beschranken wir uns im Folgenden auf die Basisoperationen von LTI-Systemen: die Multiplikation eines Signals mit einer Konstanten und die Addition von Signalen. >

Lineare Abbildung einer stochastischen Variablen

Wir betrachten zunachst die Multiplikation der stochastischen Variablen X mit der reellen positiven Zahl a. Bin anschauliches Beispiel liefem die WDFen bei der Abbildung einer gleich-

13.2 Lineare Abbildung stochastischer Signale

153

verteilten stochastischen Variablen in Bild 13-1. Wichtig ist, dass die Flachen unter den WDFen als Ma6 fiir die Wahrscheinlichkeiten vor und nach der Abbildung stets eins ergeben miissen. Jx(x)

^fyiy) Flache = 1

-2

-1

0

1

2

3

1/6

S^f^^§^S^^?5S^^

-2

4 X

Flache = 1

-1

0

1

2

3

4

J

Bild 13-1 WDF vor und nach der Abbildung Y=1X Dementsprechend ergibt sich fiir die allgemeine lineare Abbildung (13.1)

Y = a-X+b die neue WDF

\a\

(13.2)

\ a

Fiir den linearen Mittelwert und die Varianz folgt aus der Erwartungswertbildung //y = a-lUx + ^ >

und

(13.3)

Gy = a -Gx

Addition zweier stochastischer Variablen

Die Addition zweier stochastischer Variablen (13.4) fiihrt zu etwas aufwandigeren (Jberlegungen [Wer05]. Im wichtigsten Sonderfall unabhangiger stochastischer Variablen X\ und X2 resultiert die bidimensionale WDF aus der Faltung der WDFen der Summanden.

/F(y)= //xiW/x2(>'~-^)

(13.5)

^^

Fiir den Fall der Addition zweier unabhangiger stochastischer Variablen

Y=

a'X^^b'X2

(13.6)

ergibt sich unabhangig von den jeweiligen Verteilungen fiir den linearen Mittelwert und die Varianz der Summe stets /Uy = ^//j + bjLix^

und

Oy =a -az +b

-ai

(13.7)

Anmerkungen: Fur manche Anwendungen ist es nutzlich zu wissen, dass jede Linearkombination gemeinsam normalverteilter stochastischer Variablen wieder auf eine Normalverteilung fuhrt. Wendet man die FFT auf einen Block normalverteilter stochastischer Variablen an, so sind die DFT-Koeffizienten ebenfalls normalverteilt.

154 >

13 Stochastische Signale und LTI-Systeme Zentraler Grenzwertsatz

Auf der Addition unabhangiger stochastischer Variablen beruht auch der zentrale Grenzwertsatz der Wahrscheinlichkeitsrechnung. Unter relativ allgemeinen Bedingungen kann gezeigt werden, dass die (Jberlagerung einer Vielzahl in ihrer Wirkung verschwindend kleiner, zufalliger Beitrage auf eine normalverteilte Summenwirkung fiihrt. Ein physikalisches Beispiel ist das thermische Rauschen in einem Widerstand, bei dem sich die irregularen Warmebewegungen der Elektronen zu einer am Widerstand messbaren, normalverteilten Spannung iiberlagern. Wir betrachten den Sonderfall unabhangiger, identisch verteilter stochastischer Variablen X/ mit den linearen Mittelwerten ju und den Varianzen a^. Nach dem Grenzwertsatz von Lindeberg-Levy [BSMM99] erhalten wir mit der Abbildung Xi-M L vie

(13.8)

fiir A^ —> oo eine normierte, normalverteilte stochastische Variable Y. 13.2.2

Vorbereitende Aufgaben

A13.2-1 Der MATLAB-Befehl r a n d n erzeugt eine Musterfolge eines normierten GauBprozesses. Geben Sie die notigen MATLAB-Befehle an, um eine Musterfolge mit Varianz 0.5 und linearem Mittelwert 0.3 zu erzeugen.

A13.2-2 Betrachten Sie die Addition zweier unabhangiger, in [0,1] gleichverteilter stochastischen Variablen. Zeichnen Sie in Bild 13-9 die resultierende WDF der Sunmie ein. Hinweis: keine lange Rechnung erforderlich.

fix)

1 A

2

Bild 13-2 WDF der Summe zweier unabhangiger, in [0,1] gleichverteilter stochastischer Variablen 13.2.3

Versuchsdurchfuhrung Den Ausgangspunkt bilden zwei in sich unabhangige, in [0,1] gleichverteilte stochastische Prozesse Xi{n\ undX2{n\. Erzeugen Sie mit MATLAB zwei Musterfolgen x\[n\ und X2{n\. Durch y{n\ = xi[n] + X2[n] erzeugen Sie eine Musterfolge des Summenprozesses Y[n] = Xi[n] + X2[n]. Messen Sie die WDF des Summenprozesses. Welche WDF ergibt sich? Hinweis: Siehe histogram

13.3 Abbildung stochastischer Signale an LTI-Systemen

155

D

Wiederholen Sie den Versuch fiir die Addition mehrerer Prozesse, z. B. 10. Welche WDF ergibt sich jetzt naherungsweise und warum? Kontrollieren Sie auch den linearen Mittelwert und die Varianz des Summenprozesses. Geben Sie Ihr Ergebnis und die approximierende WDF graphisch aus.

n

Kontrollieren Sie mit MATLAB Ihr Ergebnis aus der Vorbereitung A13.2-1.

n

Addieren Sie nun elementweise die Musterfolgen zweier unabhangiger, normalverteilter Prozesse Xi[n] und X2[n] mit linearem Mittelwert 0.3 und Varianz 0.5. Welche WDF ergibt sich jetzt naherungsweise? Geben Sie Ihr Messergebnis und die theoretische WDF gemeinsam graphisch aus. Kontrollieren Sie auch den linearen Mittelwert und die Varianz des Summenprozesses Y[n],

13.3

Abbildung stochastischer Signale an LTI-Systemen

13.3.1

Grundlagen

Durch Filterung eines stochastischen Prozesses mit einem zeitdiskreten kausalen und reellwertigen LTI-System entsteht ein neuer Prozess, dessen KenngroBen vom Eingangsprozess und dem System abhangen. Die wichtigsten Zusammenhange zwischen den ProzesskenngroBen am Eingang und am Ausgang eines LTI-Systems sind in Tabelle 13-1 zusammengestellt. Besonders einfach ist der Fall weiBen Rauschens am Eingang, da dann die Korrelation am Systemausgang bis auf einen konstanten Faktor, der Leistung am Eingang, durch das LTISystem bestimmt wird. Am Ausgang ergibt sich die AKF des Prozesses als die Zeit-AKF der Impulsantwort RYYU] = (^x'^n^RhhU]

= ol'RhhU]

(13.9)

und flir das LDS die Leistungsubertragungsfunktion des Systems. SYY(^) = CTJc'^hh(^) 13.3.2

(13-10)

Vorbereitende Aufgaben

A 13.3-1 Ein unkorreUertes, normiertes und normalverteiltes Rauschsignal wird mit einem System mit der Ubertragungsfunktion (13.11) gefiltert. 0 3z z-0.8 Berechnen Sie die AKF des Signals am Filterausgang RyyU] im Zeitbereich aus der ^ Impulsantwort h[n]. Setzen Sie die Losungen unten ein. Hinweis: Bestimmen Sie dazu erst die Impulsantwort des Systems, s. Versuch 9, und berechnen Sie dann die Zeit-AKF durch die Faltung. Beachten Sie, dass die Zeit-AKF eines reellen LTI-Systems eine reelle gerade Funktion ist. Es ist daher vorteilhaft, die Berechnung zunachst nur fiir den Fall / > 0 durchzufuhren und dann die Losung fur / < 0 symmetrisch fortzusetzen.

156

13 Stochastische Signale und LTI-Systeme

Tabelle 13-1 Stochastische Signale und zeitdiskrete kausale und reelle LTI-Systeme

h[n]eR

reelle und rechtsseitige Impulsantwort

mit h[n] =

0\/n0] =

gesuchte AKF

«

^,,^(z) =

H{z)-H{z-^)

=

RYYU] =

Hinweis: Die Zeit-AKF eines reellwertigen LTI-Systems ist eine reelle und gerade Funktion. Reell bedeutet, dass komplexe Pole nur als konjugiert komplexe Paare auftreten. Gerade bedeutet, dass das Inverse zu jedem Pol ebenfalls einen Pol der z-Transformierten der Zeit-AKF ^hh(z) ergibt. Dabei sind die Pole im Einheitskreis der z-Ebene dem rechtsseitigen und die Pole auBerhalb dem linksseitigen Anteil der AKF zuzuordnen. Fiir die Berechnung im Bildbereich werden nur die Pole innerhalb des Einheitskreises, die kausalen Pole, bei der Bestimmung des rechtsseitigen AKF-Anteils mit der Partialbruchzerlegung als Pole ausgewertet.

Im{z} 1

i

^ -i1 | ir^ i-j

-~---\~--~\-A\----Re{zl

1 ji 1 ^ -—j-—|--r^ i

i

i

Bild 13-3 Pol-Nullstellendiagramm der z-Transformierten der AKF RyrU]

13 Stochastische Signale und LTI-Systeme

158

A13.3-3 Ein stationares, unabhangiges und in [0,1] gleichverteiltes Rauschsignal wird mit einem System mit der (Jbertragungsfunktion (13.11) gefiltert. Berechnen Sie die linearen Mittelwerte, die quadratischen Mittelwerte und die Varianzen am Systemeingang und Systemausgang. Tragen Sie die Ergebnisse in Tabelle 13-2 ein. ^ Hinweis: Beachten Sie, dass der Eingangsprozess in sich unkorreliert, aber nicht mittelwertfrei ist. Die Berechnung des linearen Mittelwertes bzw. der Varianz kann deshalb nicht direkt wie in (13.9) erfolgen, so dass eine etwas langere Rechnung erforderlich ist.

Tabelle 13-2 Abbildung einfacher KenngroBen durch das LTI-System Hi{z)

M

mi

C72

Eingang Ausgang

A 13.3-4 Geben Sie die notwendigen MATLAB-Progranmizeilen an, um bei bekannten Filterkoeffizienten, Zahler- u. Nennerkoeffizienten b bzw. a, die Zeit-AKF des Filters zu berechnen. Hinweis: Verwenden Sie die MATLAB-Befehle impz und conv. Es geniigen zwei Befehlszeilen.

A13.3-5 Geben Sie die notwendige MATLAB-Programmzeile an, um bei bekannten Filterkoeffizienten, die Zahler- u. Nennerkoeffizienten b bzw. a, die Leistungslibertragungsfunktion des Filters zu bestimmen. Hinweis: Verwenden Sie den MATLAB-Befehl f reqz.

13.3 Abbildung stochastischer Signale an LTI-Systemen

159^

13.3.3

Versuchsdurchfuhrung

D

Verwenden Sie das System 1. Ordnung mit der Ubertragungsfunktion (13.11). Das Eingangssignal sei eine stationare, unabhangige und in [0,1] gleichverteilte Zufallszahlenfolge. Bestimmen Sie die WDF am Systemausgang durch Simulation mit MATLAB. Approximieren Sie die geschatzte WDF durch die einer Normalverteilung und stellen Sie das Ergebnis und die Naherung graphisch dar. Schatzen Sie auch den linearen und den quadratischen Mittelwert und die Varianz am Systemausgang. Kontrollieren Sie die Ergebnisse anhand Ihrer Vorbereitung in Tabelle 13-2. Hinweis: MATLAB-Befehle f i l t e r , histogram

n

Schatzen Sie die AKF und das LDS zum System mit der Ubertragungsfunktion (13.11). Erregen Sie das System mit normiertem, gauBschen Rauschen, s. r a n d n . Kontrollieren Sie Ihre Rechnung? Hinweis: Verwendet man den Zufallszahlengenerators rand, ist der Gleichanteil zu beriicksichtigen. Obendrein ergeben sich relativ stark schwankende Simulationsergebnisse, weshalb hier darauf verzichtet wird.

n

Schatzen Sie die AKFen und Leistungsubertragungsfunktionen der Systeme mit den untenstehenden tJbertragungsfunktionen Hiiz) und H^iz). Diskutieren Sie die Ergebnisse. Geben Sie auch die KenngroBen der Systeme wie in Versuch 9 an. Hinweis: Verwenden Sie den Zufallszahlengenerator randn und zur Filteranalyse f v t o o l . ^^ ^ ^ 0.06z^+0.12z + 0.06 H2(z) = 5 z^ -1.3Z + 0.845

n

^ und

^ _ ^ 0.845z^-1.3z + l H^{z) = —z z^ -1.3Z + 0.845

Schatzen Sie die bidimensionale WDF der Wertepaare einer Musterfolge (x[n], x[n+l]) am Systemausgang des Filters H2iz) fiir die Verschiebungen / = 1, 2, 3, 4 und 5. Als Erregung verwenden Sie ein normiertes normalverteiltes Eingangssignal. Stellen Sie die WDFen und ihre Hohenlinien graphisch dar. Diskutieren Sie die Ergebnisse. Uberlegen Sie in welchen Bereichen Amplitudenpaare gehauft auftreten und in welchem Zusammenhang die Haufungen mit dem Korrelationskoeffizienten stehen. Hinweis: s. a. Programmbeispiele 12-2 und 12-3 mit h i s t 2 d in Versuch 12.

D

Stellen Sie zu obiger Versuchsaufgabe auch die Streudiagramme der Paare {x[n\, x{n-\-l\) am Systemausgang des Filters Hiiz) fur die Verschiebungen / =1, 2, 3, 4 und 5 dar. Als Erregung verwenden Sie wieder ein normiertes, normalverteiltes Eingangssignal. Fiir das Streudiagramm interpretieren Sie die Wertepaare als Koordinaten (jc, y) in der (jc, 3;)-Ebene. Tragen Sie dort die Wertepaare als Punkte ein. Zur besseren Ver-

160

13 Stochastische Signale und LTI-Systeme gleichbarkeit normieren Sie die angezeigten Werte einer Musterfolge jeweils auf dem Maximal wert eins. Uberlegen Sie wieder in welchen Bereichen Amplitudenpaare gehauft auftreten und in welchem Zusammenhang die Haufungen mit dem Korrelationskoeffizienten stehen.

H Hinweise zu MATLAB-Funktionen und M-Files Im Folgenden werden fiir die Versuchsdurchfiihrung nlitzliche MATLAB-Befehle und -Funktionen aufgelistet, zu denen Sie sich mit Hilfe der Help-Funktion Erlauterungen und Beispiele am Bildschirm anzeigen lassen konnen, s. a. vorherige Versuche. Tabelle 13-3 MATLAB-Befehle - benutzte Programme und Dateien, s. a. vorhergehende Versuche Graphikbefehle

scatter

Funktionen der Signalverarbeitung

conv, f f t , f i l t e r , impz, r a n d , r a n d n ,

vorbereitete Programme und Dateien

dsplabl3_l.m, dsplabl3_2.m, dsplabl3_3.m, dsplabl3_4.m, histogram.m, hist2d.m

freqz, xcorr

fvtool,

161

14

Analog-Digital-Umsetzung

14.1

Einfuhrung

In diesem Versuch werden die beiden Schritte vom analogen zum digitalen Signal vorgestellt: die Abtastung und die Quantisierung. Die Abtastung erzeugt das zeitiich diskretisierte Signal, wobei in der Regel das Abtasttheorem zu beachten ist. Die Quantisierung ist unvermeidlich, well zur Darstellung der Signalamplituden auf Digitalrechnem nur eine begrenzte Anzahl von Binarstellen zur Verfiigung steht. Lemziele: Nach Bearbeiten dieses Versuches konnen Sie •

die notwendigen Verarbeitungsschritte der Analog-Digital-Umsetzung anhand eines Blockschaltbildes skizzieren

•

das Abtasttheorem angeben und seine Bedeutung erklaren

•

die Quantisierung anhand der Quantisierungskennlinie erlautem

•

ein einfaches Modell fur den Quantisierungsfehler angeben und das SignalQuantisiemngsgerausch-Verhaltnis in Abhangigkeit der Wortlange abschatzen

•

das Gleitkomma-Format und das Festkomma-Format vorstellen und beziiglich der Dynamik und Prazision bewerten

•

Zahlen in Zweierkomplement-Darstellung und im Gleitkomma-Format IEEE 754-1985 angeben

14.2

Digitalisierung

Die prinzipiellen Verarbeitungsschritte zur Digitalisierung eines analogen Signals zeigt Bild 14-1. Die Vorfilterung durch den Tiefpass mit der Grenzfrequenz/^ kann unterbleiben, falls das Eingangssignal bereits entsprechend bandbegrenzt ist. Zunachst wird bei der zeitlichen Diskretisierung, der idealen Abtastung (Sampling Operation), jeweils alle Abtastintervalle T^ ein Abtastwert (Sample) aus dem analogen Signal entnommen. Die zeitdiskrete Abtastfolge x[n^ besitzt noch wertkontinuierliche Amphtuden.

analoges Signal 1

x\t)

zeitdiskretes Signal

digitales Signal

1

Tiefpass

x(t)

\ -^[^1 = ^it - nTs) Abtastung 1

Abtastintervall T^

Quantisierung

UinWQ ^

f Quantisierungskennlinie

Bild 14-1 Digitalisierung: Verarbeitungsschritte vom analogen zum digitalen Signal

162

14 Analog-Digital-Umsetzung

Bei der Quantisierung werden den Amplituden Werte aus einem diskreten Zeichenvorrat, den Maschinenzahlen, zugewiesen. Die Zuweisung wird durch die Quantisierungskennlinie definiert. Das digitale Signal [jc[n]](2 entsteht. Technisch wird die Digitalisierung in Analog-Digital-Umsetzern (A/D-Umsetzer) realisiert [TiSc99]. Dabei kommen meist integrierte mikroelektronische Bausteine zum Einsatz, die gemaB den Anforderungen aus einem reichhaltigen Angebot gewahlt werden konnen. Grundsatzlich wachst die Komplexitat der A/D-Umsetzer mit der Hohe der Abtastrate und Genauigkeit der Zahlendarstellung. Typische Bereiche fiir die Wortlangen sind 8 bis 16 (20) Bits und Abtastfrequenzen bis 100 (1000) MHz.

14.3

Abtastung

14.3.1

Abtasttheorem

Eine sinnvolle zeitliche Diskretisierung liegt vor, wenn die Veranderungen des analogen Signals durch die Abtastfolge gut wiedergegeben werden. Damit das analoge Signal aus der Abtastfolge durch eine Interpolation hinreichend genau wieder gewonnen werden kann, muss ein sich schnell andemdes Signal haufiger als ein dazu relativ langsam veranderliches Signal abgetastet werden. Diese grundsatzliche (Jberlegung wird im Abtasttheorem prazisiert. Abtasttheorem Eine Funktion x{t), deren Spektrum fiir | / | > /^ null ist, wird durch die Abtastwerte x{t = nT^) vollstandig beschrieben, wenn das Abtastintervall T^, bzw. die Abtastfrequenzfs, so gewahlt wird, dass 1 '

fs

(14.1)

Die so abgetastete Funktion kann aus den Abtastwerten durch die si-Interpolation fehlerfrei rekonstruiert werden. 00

x{t)= 2J x{nT^)- si(/,;r[/ -«?;])

(14.2)

«=-00

Beachten Sie die hier gewahlte FormuHerung des Abtasttheorems. Bei der Grenzfrequenz/g ist das Spektrum des analogen Signals bereits zu null angenommen. Denn soil das Abtasttheorem mit /^ = 2fg eingehalten werden, ist ein Spektralanteil bei der halben Abtastfrequenz auszuschUeBen. Die zur Interpolation verwendeten si-Impulse entsprechen im Frequenzbereich einem idealen Tiefpass mit der Grenzfrequenz fg. Eine Interpolation mit einem idealen Tiefpass liefert das urspriingliche zeitkontinuierliche Signal. Die praktische Anwendung der si-Interpolation geschieht naherungsweise in Digital-Analog-Umsetzern: Ein Abtast-Halteglieder erzeugt zunachst ein treppenformiges analoges Signal, das anschlieBend durch einen Tiefpass geglattet wird. Die Auswirkung der Abtastung auf die Spektren der Signale wird in Versuch 5 KurzzeitSpektralanalyse: Grundlagen vorgestellt. Wird das Abtasttheorem verletzt, kommt es zur Uberfaltungen des periodisch fortgesetzten Spektrums, d^m Aliasing.

14.3 Abtastung 14.3.2

163^

Vorbereitende Aufgaben

A14.3-1 Mit dem Programmbeispiel 14-1 wird im Versuch der Einfluss der zur Abtastung notwendigen Bandbegrenzung horbar gemacht. Ein Audiosignal wird bandbegrenzt, unterabgetastet und auf der PC Sound Card ausgegeben. Machen Sie sich mit dem Programm vertraut und erklaren Sie insbesondere wie die Bandbegrenzung durchgeflihrt wird. Hinweis: Da nicht jede PC Sound Card eine weitgehend beliebige Einstellung der Abtastfrequenz erlaubt, wird im Programmbeispiel auf die explizite Unterabtastung verzichtet. Es wird die vor der Unterabtastung notwendige Bandbegrenzung durchgefuhrt, so dass sich ein der Unterabtastung entsprechender Horeindruck ergibt. Programmbeispiel 14-1 Unterabtastung eines Audiosignals (Tiefpass-Filterung) % subsampling of an audio signal using lowpass filtering by fft % dsplabl4_l.m * mw * 01/23/2006 fprintf(•01/23/2006/dsplabl4_l : subsampling \n') filename = input('name of audio file (*.wav): ','s'); [x,fs,bits] = wavread(filename); % load audio signal X = x/max(abs(x)); % scaling X = fft(x); Lx = length(x); K = 0; while K~=7 K = menu([filename,' fs = ',num2str(fs),' Hz'],'original',... 'fs / 2','fs / 4 •,'fs / 8 '^'fs / 16 '^'fs / 32 '^'exit'); SR = 2^(K-1); % subsampling factor

% lowpass

filtering

using

fft

index = ceil(Lx/(2*SR)); XO = [X(l:index); zeros(Lx+l-2*index,1); X(Lx+2-index:Lx)]; x_lp = real(ifft(XO)); if K~=7 soundsc(x_lp,fs,bits); end end

14.3.3

Versuchsdurchfiihrung

n

Machen Sie sich die Bedeutung der Abtastfrequenz mit dem Programmbeispiel 14-1 klar, indem Sie sich das Sprachsignal s p e e c h . w a v bei verschiedenen Abtastfrequenzen anhoren. Bis zu welcher Grenzfrequenz ist das Sprachsignal noch verstandlich?

D

Erzeugen Sie Ausschnitte eines Kosinussignals mit der jeweiligen Dauer von 0.5 s bei der Abtastfrequenz 8 kHz. Variieren Sie die Signalfrequenz von 1 kHz bis 8 kHz und geben Sie die Signale kurz hintereinander liber die PC Sound Card aus. Welcher Effekt stellt sich? Diskutieren Sie das Ergebnis. tjberlegen Sie: Was geschieht bei der Abtastung eines Kosinussignals mit der Frequenz 4 kHz bei einer Abtastfrequenz von 8 kHz? Welche Moglichkeiten ergeben sich fur die Abtastfolge, wenn die Anfangsphase der Abtastzeitpunkte zufallig ist? Welche Information geht dabei verloren?

164

14 Analog-Digital-Umsetzung

14.4

Quantisierung

14.4.1

Quantisierungskennlinie

Das analoge Signal x(t) sei auf das Amplitudenintervall [-1,1 [ begrenzt. Falls nicht, wird das Signal vor der Quantisierung geeignet skaliert. Im Weiteren wird stets von einem Quantisierungsbereich von -1 bis +1 ausgegangen. Zur Darstellung der Abtastwerte sollen je w Bits zur Verfugung stehen. Man spricht von der Wortldnge und schreibt z. B. kurz w = 3. Mit w Bits konnen genau 2^ Quantisierungsintervalle dargestellt werden. Bei der gleichformigen Quantisierung teilt man den Quantisierungsbereich in 2^ gleichgroBe Intervalle mit der Quantisierungsintervallbreite 2 = 2'^^"^)

(14.3)

Anmerkung: Ublich sind auch die Sprechweisen Quantisierungsstufe und Quantisierungsstufenhohe. Die exakte Beschreibung der Quantisierung geschieht durch die Quantisierungskennlinie, wie z. B. in Bild 14-2. Sie definiert die Abbildung der kontinuierlichen Abtastwerte auf die Reprdsentanten. Die Reprasentanten sind in Bild 14-2 als kleine Quadrate kenndich gemacht. Bei der gleichformigen Quantisierung mit Runden liegen die Reprasentanten jeweils in der Mitte der Quantisierungsintervalle, so dass der Abstand zwischen Abtastwert und Reprasentant die halbe Intervallbreite nicht iiberschreitet. Anhand von Bild 14-2 lassen sich die zwei grundsatzlichen Probleme der Quantisierung erkennen: ® Eine Ubersteuerung tritt auf, wenn das Eingangssignal auBerhalb des vorgesehenen Aussteuerungsbereichs liegt. In der Regel tritt dann die Sdttigung ein und es wird der Maximalwert bzw. der Minimalwert der Reprasentanten ausgegeben (Sdttigungskennlinie). ® Eine Untersteuerung Hegt vor, wenn das Eingangssignal (fast) immer viel kleiner als der Aussteuerungsbereich ist. Im Extremfall entsteht so genanntes granulares Rauschen bei dem das quantisierte Signal scheinbar regellos zwischen den beiden Reprasentanten um die Null herum wechselt. Dieses Problem kann umgangen werden, wenn ein Quantisierungsintervall um die Null eingefuhrt wird, wie unten in Bild 14-2. Bei jeder Anwendung ist stets wichtig, den Eingang des A/D-Umsetzers passend auszusteuem. Zu den grundsatzlichen Problemen konnen im praktischen Einsatz durch nicht perfekte A/DUmsetzer weitere statische und dynamische Fehler hinzutreten [TiSc99]. 14.4.2

Maschinenzahlen

Bei der Zahlendarstellung auf Digitalrechnem, den Maschinenzahlen, wird in der Regel das Dualsystem mit der Basis 2 zugrunde gelegt. Es kommen zwei grundsatzliche Formate zum Einsatz: das Gleitkomma-Format und das Festkomma-Format.

14.4 Quantisierang

165

W(2^

Sattigungs kennlinie

k

1

Reprasentant 1/2 i

1

X.''

Beispiele

i^

-f-

1/4 i-;>^«(____]____-

ol/ i -1

'

-1/2 li^'

1-1/2

—i—1

1-3/4 I* -

A

1/2

1 1-1/4

I—i—nf-_ j — I

Y'

^ A' 1

X

MQ

0.7

0.5

0

0

-0.2

-0.25

-1.1

-1

1

Sattigungs kennlinie 3/4 1 1/2 [ 1 — i - «

X

W(2

1.1

0.75

0.7

0.75

1-1/2

0

0

1-3/4

-0.2

0

1/4 1 l-^-T^^

-1

Beispiele

H-kL-i—

--!-—i— 1

-1/2

^ X

* 7 ^ - < 1-1/4 L____U^^ ^

•''' -^-—^-

- ^ — I

__J

£

1 -1

Bild 14-2 Quantisierungskennlinie mit Abschneiden (oben) und Runden (unten) bei gleichformiger Quantisierung mit der Wortlange 3 Bits und Zweierkomplement-Darstellung •

Gleitkomma-Format

Im Gleitkomma-Format werden die Zahlen durch das Vorzeichen, den Exponenten und die Mantisse dargestellt, s. Bild 14-3. Sie werden tiblicherweise so skaliert, dass die hochstwertige Stelle der Mantisse eins ist. Sie kann deshalb weggelassen (gemerkt) werden und so effektiv eine Stelle mehr zur Verfugung steht. Man spricht von der normalisierten Form. Im Beispiel des weit verbreiteten 64-Bit-Formats (Double Precision) nach IEEE 754-1985 bedeutet das, s. a. Bild 14-3,

x = (-\y '2^ '{\ + Mf)

(14.4)

mit dem Exponenten mit Bias 10

£ = 2]^/-^^-1023

(14.5)

/=0

und der Mantisse (Significand) 51

-52

•Ifr i=0

(14.6)

14 Analog-Digital-Umsetzung

166

31 30

24 22

Single Precision P"T^^ 63 162

Double Precision

J 2 51

^i)^>^n^:^--'-

Vorzeichen IBit

/22 ••• /O

/l /o

/51/5O

Mantisse 23 / 52 Bits

Exponent 8/11 Bits

Bild 14-3 Gleitkommadarstellung nach IEEE 754-1985 fiir Single Precision und Double Precision Wir machen uns die Anwendung an Beispielen deutlich. Beispiel 14-1 Gleitkomma-Format einer nattirlichen Zahl Im Beispiel der Zahl 100 ergibt sich zu (14.4) zunachst die Zerlegung 100 = (-1)^ • 2^ • (1 + 0.5625)

(14.7)

und daraus das zugehorige Bitmuster

100^ =oioo'oooo'oioiiooroooo...oooo "

U I

^ = 6+1023

I I

My =0.5625

(14.8)

I

Beispiel 14-2 Gleitkomma-Format einer negativen Dezimalzahl Im Beispiel der Zahl -0.40625 fuhrt die Zerlegung (14.4) auf (14.9)

-0.40625 = (-1)^ -2"^ -(1 + 0.625) und das zugehorige Bitmuster ist -0.40625^ = 1 0 i n i l l ' 1 1 0 1 1010'0000...0000 "

U I

s

^ = -2+1023

I

I

My =0.625

(14.10)

I

MATLAB am PC verwendet das Format IEEE 754-1985, das bestimmte Bitmuster fur Ausnahmen reserviert, wie zum Beispiel I n f (positiver tJberlauf, 00), - i n f (negativer Uberlauf, - 00) und NaN fiir Not a Number. Letzteres resultiert als Ergebnis fiir undefinierte Ausdrticke, wie z. B. bei Division von 0 / 0 . Um die Ausnahmen anzuzeigen ist der maximal moghche Exponent e^^x = 2^-1 = 255 bzw. 211-1 = 2047 reserviert, s. Tabelle 14-1. Man beachte auch den Unterschied zwischen den normahsierten und denormalisierten Zahlen. Die kleinste von null verschiedene positive Zahl in normalisierter Form ist 2^^^. Durch denormalisieren kann diese Zahl auf 2-23 • 2^min bzw. 2-52 . 2^niin verkleinert werden. So wird das Problem des Zahlenunterlaufs wirksam reduziert.

14.4 Quantisierung

167

Tabelle 14-1 Codierung von Gleitkommazahlen (IEEE 754-1985) reelle Zahl

Exponent

Mf

Mantisse

I — e < ^niax

fm^ -"Jo

1+My

(-ly '2^ '{l + Mf)

denormalisiert

0, ...,0

fm^ •••'/o

Mf

(_l)^.2^min .Mj-

01

0, ...,0

0, ...,0

-

0

NaN

^max

1,..., 1

-

-

Inf, - I n f 1

^max

0, ...,0

-

-

Typ normalisiert

E = e - Bias

0 0 und -0 werden ebenso wie Inf und -Inf durch das Vorzeichenbit 5 = 0 bzw. s = -\ unterschieden Formatbeispiele

m

^max

F

Single Precision

22

255

-126

127

127

Double Precision

51

2047

-1022

1023

1023

^max ~ ^max

^

-Bias

Bias

Festkomma-Format Flir die Darstellung von Festkomma-Zahlen gibt es verschiedene Moglichkeiten: Betrag und Vorzeichen, Einerkomplement und Zweierkomplement. Wir beschranken uns im Folgenden auf das in der digitalen Signalverarbeitung haufig verwendete Festkomma-Format im Zweierkomplement (Two's Complement) mit der Wortlange w Bits. Im Zweierkomplement-Format gibt es genau 2^ verschiedene Maschinenzahlen, wobei die Null eindeutig ausgedriickt wird. w—1

jc = -flo2^ + 2 ^ / ^ ' ™^ a,E{0,l}

f u r - l < x < l - 2-w+\

(14.11)

i=\

Die Negation geschieht vorteilhaft durch Komplementbildung und Addition eines Bits mit geringster Wertigkeit, dem LSB {Least Significant Bit). w—\

x = -flo2^ + 2^a,-2"^+2"''"^^ m i t a , G { 0 , l }

f u r - l < x < l - 2 " -w-hl

(14.12)

Anmerkungen: Das Zweierkomplement-Format ermoglicht prinzipiell die Darstellung der Zahl -1. Bei manchen Anwendungen, speziell mit digitalen Signalprozessoren, wird aus Griinden der Kompatibilitat auf die Verwendung der -1 verzichtet. Das Zweierkomplement-Format ermoglicht die einfache Addition von positiven und negativen Zahlen.

14 Analog-Digital-Umsetzung

168 Beispiel 14-3 Festkomma-Format einer negativen Dezimalzahl

Im Beispiel der Zahl -0.328125 ergibt sich aus (14.11) zunachst die Darstellung des Betrages. |-0.328125| = 2"^ + 2"^ + 2"^

(14.13)

Bei einer Wortlange von 16 Bits ist dann das zugehorige Bitmuster mit dem LSB rechts |-0.328125|^ =0010'1010'0000'00002c

(14.14)

Durch Zweierkomplement-Bildung ergibt sich daraus -0.328125^ =110r0110'0000'00002^

(14.15)

Einen pauschalen Vergleich fiir die praktische Anwendung der beiden prinzipiellen Formate bei gleicher Wortlange zeigt Tabelle 14-2. Was die numerischen Eigenschaften betrifft werden in den Anwendungen meist zwei Kriterien herangezogen: •

die Dynamik, das Verhaltnis von groBter und kleinster positiver darstellbarer Zahl

•

die Prdzision, die Feinheit der Zahlendarstellung, d. h. je groBer das Quantisierungsintervall desto kleiner die Prazision Tabelle 14-2 Vergleich des Festkomma- und Gleitkomma-Formats fur Digitalrechner bei gleicher Wortlange Festkomma-Format Aussteuerungsbereich (Dynamik)

kleiner

groBer

Quantisierungsfehler (Prazision)

konstant

exponentenabhangig

meist erforderlich

meist nicht erforderlich

Hardware-Komplexitat

kleiner

groBer

Energieaufnahme

kleiner

groBer

schwieriger

einfacher

niedriger

hoher

Skalierung

Software-Portabilitat Preis

14.4.3

Gleitkomma-Format

Quantisierungsfehler

Aus den quantisierten Werten kann das ursprlingliche Signal bis auf kiinstliche Spezialfalle nicht mehr fehlerfrei rekonstruiert werden. Der Quantisierungsfehler wird durch die Wortlange kontrolliert. Je groBer die Wortlange, umso kleiner der mogliche Quantisierungsfehler, desto groBer die Prazision der Signaldarstellung. Mit wachsender Wortlange nehmen jedoch auch die Komplexitat der arithmetischen Operationen und der Speicherbedarf in der Signalverarbeitung zu. In praktischen Anwendungen ist deshalb je nach Aufgabenstellung zwischen der gewunschten Prazision und dem notwendigen Aufwand abzuwagen. Um im konkreten Fall die Frage nach der Prazision zu beantworten, muss sie quantitativ messbar sein. Dazu verwendet man das Modell einer additiven Stoning, einem Fehlersignal, hier Quantisierungsgerdusch genannt, s. Bild 14-4.

14.4 Quantisierung

169

Anmerkung: Die (mathematische) Quantisierung wirkt auf die Folgenelement bzw. Momentanwerte unabhangig von friiheren oder folgenden Werten. Sie ist deshalb gedachtnislos. Der Einfachheit halber betrachten wir im Folgenden ein zeitkontinuierliches dreieckformiges Eingangssignal. Ein iibersichtliches Beispiel liefert die Quantisierung mit Runden des periodischen, dreieckformigen Signals x(t) in Bild 14-5. Dann wird am Ausgang die Quantisierungskennlinie sichtbar.

Eingangs- Quantisierung quantisierter ^ ^ Wert wert

e

[X\Q

A=[x] -X

Bild 14-4 Ersatzmodell fiir die Quantisierung mit dem Quantisierungsgerausch A Im unteren Bild ist das resultierende Fehlersignal A(t) aufgetragen. Betrachtet man den Zeitpunkt r = 0, so ist x(0) = 0 und [JC(0)](2 ebenfalls 0. Mit wachsender Zeit steigt das Eingangssignal linear an und das Fehlersignal fallt zunachst linear. Bei t = TQ/16 uberschreitet das Eingangssignal die Grenze des Quantisierungsintervalls. Danach ist das Eingangssignal zunachst kleiner als der zugewiesene Reprasentant. Das Fehlersignal springt von -Q/2 auf +Q/2, um iiber das gesamte Quantisierungsintervall linear wieder bis -Q/2 zu fallen. In der Mitte des Quantisierungsintervalls ist der Quantisierungsfehler null. Entsprechendes kann fiir die anderen Signalabschnitte liberlegt werden. Man beachte, dass in Bild 14-5 das Signal nicht vollstandig zwischen -1 und 4-1 ausgesteuert wird. Wegen der exakten Darstellung der Null im Zweierkomplement-Format fehlt die Darstellung der Eins. Im Falle eines Eingangswertes eins ergabe sich der Quantisierungsfehler -Q. Fiir Eingangswerte groBer gleich 1-2-^ wird der groBte darstellbare Wert l-2-(w-l) zugewiesen (Sattigung). Obige einfache Uberlegungen ermoglichen es, die Qualitat der Quantisierung quantitativ zu erfassen. Als QualitatsmaB wird das Verhaltnis der Leistungen des Eingangssignals und des Quantisierungsgerausches, das Signal- Quantisierungsgerausch- Verhdltnis, kurz SNR (Signal-to-Noise Ratio), zugrunde gelegt. Im Beispiel des dreieckformigen Signals ergibt sich bei Vollaussteuerung die mittlere Signalleistung 1

1

A(t) = [x(t)]Q - x(t) QI2

Bild 14-5 Quantisierung eines periodischen dreieckformigen Signals (oben) und das dabei entstehende Quantisierungsgerausch (Fehlersignal) A{t) (unten)

Die mittlere Leistung des Quantisierungsgerausches kann wegen des abschnittsweise linearen Verlaufs ebenso einfach berechnet werden, s. Bild 14-5. Der Sattigungseffekt der Quantisierung fiir Eingangswerte bei Vollaussteuerung nahe 1 wird - geniigend groBe Wortlange vorausgesetzt - fiir die vereinfachten Modelliiberlegungen vemachlassigt.

170

14 Analog-Digital-Umsetzung

Dann sind die Werte des Quantisierungsgerausches auf das Intervall [-2/2, Q/2] beschrankt. Die mittlere Leistung ergibt sich demzufolge wie in (14.16) zu

^ = i . ( e f = 21=31!: 3 [2)

12

(14.17,

3

wobei die Quantisierungsintervallbreite durch die Wortlange ersetzt wurde. Fiir die vereinfachten Modelltiberlegungen resultiert das SNR im logarithmischen MaB = 10-logio2^'^ dB = 20-wlogio2 d B « 6 - w dB

(14.18)

Das SNR verbessert sich demzufolge um 6 dB pro zusatzlichem Bit Wortlange. Anhand der Modellrechnung erkennt man, dass das SNR von der Art des Signals abhangt; genauer von der Verteilung der Signalamplituden. Im Falle eines im gesamten Aussteuerungsbereich gleichverteilten Eingangsprozesses ergibt sich das SNR wie in (14.18). Das gefundene Ergebnis liefert eine brauchbare Naherung fiir viele praktische Anwendungsfalle: Fiir eine symmetrische gleichformige Quantisierung mit hinreichender Wortlange w in Bits und VoUaussteuerung gilt die 6dB-pro-Bit-Regel ^/A^|^B^6-wdB

(14.19)

Fiir die Naherung liegt erfahrungsgemaB eine hinreichende Wortlange vor, wenn das Signal mehrere Quantisierungsintervalle durchlauft. Viele Signale, wie Sprach- und Audiosignale, haben einen hohen Dynamikumfang, wobei die kleinen Amplituden relativ haufig vorkommen. Bei VoUaussteuerung ist die Signalleistung deshalb deutlich kleiner als 1/3, so dass das tatsachliche SNR um einige dB kleiner ist als in (14.19).

14.4.4

Vorbereitende Aufgaben

A14.4-1 Uberlegen Sie, wie Sie mit einem MATLAB-Programm die kleinste Maschinenzahl groBer eins bestinmien konnen. 's. Hinweis: (l+£) mit £> 0, wenige Programmzeilen gentigen. A14.4-2 Geben Sie fiir die Wortlange w in Bits die fehlenden Eintrage zum Zweierkomplement-Format in Tabelle 14-3 an. A14.4-3 Geben Sie fiir die Zahlenwerte in Tabelle 14-4 im Zweierkomplement-Format fiir die Wortlange w = 8 (Bits) in hexadezimaler Form an. Runden Sie gegebenenfalls die letzte Stelle. Tragen Sie die quantisierten Werte auch als iibliche Gleitkommazahlen ein.

14.4 Quantisierung

171

Tabelle 14-3 Zweierkomplement-Format wBits

Wortlange

8 Bits

16 Bits

kleinste positive Zahl (LSB)

groBte positive Zahl (1-LSB)

kleinste negative Zahl

Dynamik

Prazisioni

0 maximaler Quantisierungsfehler bei Runden (aussteuerungsunabhangig)

Tabelle 14-4 Quantisierung im Zweierkomplement-Format w = 8 (Bits) mit Hexadezimaldarstellung (Hex) der Bitmuster JC

^2c

Wfi

X

0.996

7FHex

0.9921875

0

0.125

-0.125

0.004

-1

^2c

Me

A14.4-4 Machen Sie sich mit der Zahlendarstellung nach IEEE 754-1985 in MATLAB vertraut. Bestimmen Sie dazu fiir die normalisierte Form den groBten und den kleinsten positiven Zahlenwert, s. Tabelle 14-5. Geben Sie auch die kleinste Maschinenzahl groBer eins an. Hinweis: Geben Sie zunachst die theoretischen Werte an mit der Wortlange des Exponenten Wg und der Mantisse w^, Spezialisieren Sie dann auf die angegebenen Formate. Verwenden Sie gegebenenfalls geeignete Schatzwerte in der Dezimal-Darstellung.

172

14 Analog-Digital-Umsetzung

Tabelle 14-5 Zahlendarstellung IEEE 754-1985 (normalisiert) theoretischi

Single Precision

Double Precision

groBte positive Maschinenzahl kleinste (normalisierte) Maschinenzahl > 0 kleinste (normalisierte) Maschinenzahl^ > 1 Dynamik

Prazision^

1) Wortlange des Exponenten w^ und Wortlange der Matisse w^ 2) Maschinen-Epsilon 3) maximaler Quantisierungsfehler bei Runden aufgrund der Wortlangenbeschrankung der Mantisse (aussteuerungsunabhangig) A14.4-5 Machen Sie sich mit dem Programmbeispiel 14-2 zur Simulation einer Quantisierung im Zweierkompliment-Format vertraut. Programmbeispiel 14-2 Quantisierung im Zweierkomplement-Format function xq = quant2c(x,w,TMode) % tow's complement quantizer within the range of -1 to 1 (saturation) % function xq = quant2c(x,w,TMode) % X : input signal % w : word length (# of bits) % TMode : truncation mode % 't' - truncation (rounding towards minus infinity) % -J-' - rounding to nearest quantization level xq : quantized signal % % quant2c.m * mw * 01/25/2006

LSB = 2'^(-w+l); xq = min(l-LSB,x) ; xq = max(-1,xq); if TMode=='t' xq = floor(xq/LSB)*LSB; else xq = round(xq/LSB)*LSB; end

% least significant bit % clipping (saturation)

% truncation % rounding

14.4 Quantisierung

173

14.4.5

Versuchsdurchfuhrung

n

Bestimmen Sie die kleinste von eins verschiedene Maschinenzahl mit Ihrem vorbereiteten Programm.

D

Bestimmen Sie die MATLAB-Variablen e p s , r e a l m a x und r e a l m i n und machen Sie sich deren Bedeutung klar. Vergleichen Sie die Werte mit den entsprechenden Werten aus Ihrer Vorbereitung.

eps realmin realmax

n

Nehmen Sie die Quantisierungskennlinien zum Zweierkomplement-Format mit der Wortlange von drei Bits fiir Abschneiden und fiir Runden auf. Stellen Sie auch die Quantisierungsfehler graphisch dar. Hinweis: Verwenden Sie das Programm d s p l a b l 4 _ 2 in Programmbeispiel 14-3 mit der Funktion q u a n t 2 c . m

Programmbeispiel 14-3 Testprogramm fur die Quantisierung im Zweierkomplement-Format mit Runden mit dem Progranun q u a n t 2 c % test of tow's complement quantizer function quant2c % dsplabl4_2.m * mw * 01/25/2006 n=-l:.001:1;x=n; % input signal xq = quant2c(x,3,'r'); % quantization with rounding % graphics FIGl=figure('Name','dsplabl4_2 : quantization ','NumberTitle','off'); subplot(2,1,1),plot(n,xq,n,xq,'.',n,x),grid xlabel('n \rightarrow') ylabel('[x]_{Q} \rightarrow') subplot(2,1,2),plot(n,xq-x),grid xlabel('n \rightarrow') ylabel('\Delta \rightarrow')

D

Quantisieren Sie ein Sinussignal mit Amplitude eins mit den beiden Quantisierungskennlinien aus der vorhergehenden Aufgabe fiir die Wortlange von vier Bits. Stellen Sie das Signal vor und nach der Quantisierung und den zugehorigen Quantisierungsfehler graphisch dar. Hinweis: Modifizieren Sie das Programm dsplabl4_2 . m entsprechend.

n

Messen Sie das Signal-Quantisierungsgerausch-Verhaltnis (SNR) fiir ein lineares und ein sinusformiges Signal in Abhangigkeit von der Wortlange. Verwenden Sie die Wortlangen von 2 ... 16 Bits. Wiederholen Sie die Messung fiir das Sprachsignal s p e e c h . wav. Stellen Sie die Ergebnisse in dB in einem Bild dar. Vergleichen Sie Ihr Ergebnis mit der Abschatzung durch das Fehlermodell, der 6dB-pro-Bit-Regel. Hinweis: Verwenden Sie Programm dsplabl4_3 .m in Programmbeispiel 14-4.

174

14 Analog-Digital-Umsetzung

Programmbeispiel 14-4 Messimg des SNR fiir die Quantisierung im ZweierkomplementFormat mit Runden flir ein lineares und sinusformiges Signal und ein Sprachsignal % tow's complement quantization of linear and sinusoidal signal and % speech signal % dsplabl4_3.m * mw * 01/25/2006 n = 0:.001:1; L = length(n); w = [2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]; % word length % linear input signal SNR = zeros(size(w)); X = n; % linear signal for k=l:length(w) % signal power S = sum(x. ^^2 )/L; xq = quant2c(x,w(k),'r'); % quantized signal (rounding) N = sum((xq-x).^2)/L; % noise power SNR(k) = S/N; % signal-to-noise ratio end % graphics FIGl=figure('Name','dsplabl4_3 : quantization SNR','NumberTitle','off); plot(w,10*loglO(SNR),'o',w,10*loglO(SNR),':',w,6*w), grid xlabel('word length w [bits] \rightarrow') ylabel('SNR [dB] \rightarrow')

n

Untersuchen Sie den Einfluss der Quantisierung auf Audiosignale durch Hortests. Wiederholen Sie die Hortests fiir das Sprachbeispiel s p e e c h . w a v . Bis zu welcher ]V[indestwortlange ist das Sprachsignal noch verstandlich? Probieren Sie auch andere Audiosignale aus. Hinweis: Skalieren Sie den Aussteuerungsbereich des Sprachsignals geeignet.

H Hinweise zu MATLAB-Funktionen und M-Files Im Folgenden werden fiir die Versuchsdurchfuhrung ntitzliche IVIATLAB-Befehle und -Funktionen aufgelistet, zu denen Sie sich mit Hilfe der Help-Funktion Erlauterungen und Beispiele am Bildschirm anzeigen lassen konnen, s. a. vorherige Versuche. Tabelle 14-6 MATLAB-Befehle - benutzte Programme und Dateien, s. a. vorhergehende Versuche spezielle Befehle und Konstanten

eps, realmax, realmin, format

vorbereitete Programme und Dateien

dsplabl4_l.m, dsplabl4_lb.m, dsplabl4_lc.m, dsplabl4_2.m, dsplabl4_2b.m, dsplabl4_2c.m, dsplabl4_3 .m, dsplabl4_4 .m, quant2c.in, speech.wav, guitar.wav, handel.wav

175

15

Reale digitale Filter: Koeffizientenquantisierung

15.1

Einfiihrung

In der Systemtheorie werden digitale Signale und Systeme zunachst unter idealisierten Bedingungen betrachtet. Reale Systeme, wie beispielsweise bei der MATLAB-Simulation am PC, bei einer Implementierung auf einem digitalen Signalprozessor oder einer dedizierten Hardware, arbeiten jedoch stets mit endlicher Wortlange. Es entstehen zwangslaufig Quantisierungsfehler, wenn die Koeffizienten und Variablen auBerhalb des Darstellungsbereiches der Maschinenzahlen liegen. Lemziele: Nach Bearbeiten dieses Versuches konnen Sie •

die Fehlerquellen bei realen digitalen Filtem nennen

•

eine Abschatzung fur die maximale Frequenzgangsabweichung durch Quantisierung der Koeffizienten bei FIR-Filtem angeben

•

mit MATLAB den tatsachlichen Fehler durch Quantisierung der Koeffizienten bei FIR-Filtem bestimmen

•

die Wirkung der Koeffizientenquantisierung von IIR-Filtem auf den Frequenzgang anhand des Pol- NuUstellendiagramms erlautem

•

den Effekt der Polausdtinnung anhand einer Skizze in der z-Ebene erlautem

•

ein IIR-Filter hoherer Ordnung in eine fur die Realisiemng giinstige Kaskadenform umwandeln

•

mit MATLAB den tatsachlichen Fehler durch Quantisiemng der Koeffizienten bei IIR-Filtem bestimmen

•

das MATLAB Filterdesign-und-analyse-Werkzeug f da t o o l gezielt einsetzen

15.2

Wortlangeneffekte

Praktische Systeme zur digitalen Signalverarbeitung werden haufig auf speziellen Mikrokontrollem iplementiert, den digitalen Signalprozessoren (DSP, Digital Signal Processor). Nach der Art der Zahlendarstellung werden zwei Gmppen unterschieden: FestkommaProzessoren und Gleitkomma-Prozessoren. Wichtige Vor- und Nachteile der beiden ProzessorArchitekturen stellt Tabelle 15-1 gegentiber. Beim praktischen Einsatz digitaler Systeme sind verschiedene Fehlerquellen zu beriicksichtigen: • •

Quantisierungsfehler bei der A/D- Umsetzung Quantisierungsfehler bei der Darstellung der Filterkoeffizienten ® mogliche Folgen: Entwurfsspezifikationen werden verletzt. Filter wird instabil © Filter bleibt ein lineares System

176 •

15 Reale digitale Filter: Koeffizientenquantisierung

Arithmetikfehler innerhalb des Filters treten bei Wortlangenverkiirzung nach einer Multiplikation oder einem Uberlauf nach einer Addition auf ® mogliche Folgen: inneres Gerausch, kleine und groBe Grenzzyklen ® Filter ist kein lineares System mehr

In diesem Versuch werden die Auswirkungen der Koejfizientenquantisierung auf das Ubertragungsverhalten digitaler Filter untersucht. Dabei wird nach FIR- und IIR-Filtem unterschieden, da letztere wegen ihrer rekursiven Struktur fehleranfalliger sind. Die Arithmetikfehler werden in Versuch 16 behandelt.

Tabelle 15-1 Vergleich der Festkomma- und Gleitkomma-Architektur Festkomma-Architektur

Gleitkomma-Architektur

Dynamik

kleiner

groBer

Prazision

konstant

exponentenabhangig

meist erforderlich

meist nicht erforderlich

einfacher groBer meist groBer

umfangreicher kleiner meist kleiner

groBer kleiner kleiner

kleiner groBer groBer

Skalierung Befehlssatz (Assembler) Optimierungsmoglichkeiten Entwicklungszeit Geschwindigkeit Stromverbrauch Preis

Sonstiges

zur Vermeidung von Arithmetikfehlem besitzen modeme FestkomSoftware von PCs und Arbeitsma-DSPs haufig spezielle Archiplatzrechnem auf den DSP tekturmerkmale, wie Sattigungsleichter portierbar kennlinien, Akkumulatoren mit vergroBerter Wortlange

15.3

FIR-Filter mit quantisierten Koeffizienten

15.3.1

Fehlermodell und Fehlerfrequenzgang

Bei einem nichtrekursiven Filter der Ordnung A^ sind die Filterkoeffizienten gleich der Impulsantwort, s. Bild 15-1. Die Quantisierung der Koeffizienten wirkt sich direkt auf die Impulsantwort und damit auf die LFbertragungsfunktion und den Frequenzgang aus. Fasst man den Effekt der Koejfizientenquantisierung als additive StorgroBe /zj[n] auf h[n] = [h[n\]Q+hd[n\

(15.1)

so resultiert das Modell einer Parallelschaltung aus dem nichtquantisierten System und des Fehlersystems durch die Koeffizientenquantisierung wie in Bild 15-2. Wegen der Additivitat linearer Systeme gilt fiir den Frequenzgang des quantisierten Systems

HQ(ei^) =

H[ei^)-H,(eJ^)

(15.2)

15.3 FIR-Filter mit quantisierten Koeffizienten

177

Bild 15-1 FIR-Filter in transversaler Struktur Bel einer Realisierung im Festkomma-Format liegen meist Zahlen im Zweierkomplement-Format vor. Sie beschrankt den Zahlenbereich auf Werte zwischen -1 und 1-LSB (LSB, Least Significant Bit). Der maximale Quantisierungsfehler betragt bei Runden LSB/2.

x[n] —

Damit kann die Abweichung vom Frequenzgang aufgrund der Koeffizientenquantisierung abgeschatzt wer- Bild 15-2 Parallelschaltung aus nichtden. Fine obere Schranke der Betragsfrequenzgangsabquantisiertem System und weichung liefert die Annahme, dass sich alle QuantisieFehlersystem rungsfehler der Koeffizienten ungiinstig uberlagem. Bei einer Quantisierung der Koeffizienten mit der Wortlange w in Bit und Runden ergibt sich demzufolge max H^[eJ^)

| < ( ^ + l ) - ^ = (A^ + l)-2'

(15.3)

Die Abschatzung erweist sich bei realen Anwendungen - wie in der Versuchsdurchfiihrung als pessimistisch. Die Fehler durch die Quantisierung der Koeffizienten kompensieren sich zum Teil. Anmerkung: Realistischere Abschatzungen in der Literatur beriicksichtigen diesen Effekt. Sie fuBen auf stochastischen Modellen fiir die Quantisierungsfehler. Wird das Zweierkomplement-Format verwendet, ist der Wertebereich der Signale und inneren GroBen des Systems zwischen -1 und 1-LSB begrenzt. Zur Vermeidung von kritischen Uberlaufen werden haufig die Betragsfrequenzgange auf den Maximalwert eins skaliert. max

H(e'-) 1 =

(15.4)

Dann wird keine Frequenzkomponente des Eingangssignals verstarkt. Die Uberlagerung unterschiedlich phasenverschobener Frequenzkomponenten im System kann jedoch zu Uberlaufen fuhren.

15.3.2

Vorbereitende Aufgaben

A 15.3-1 Bei der Versuchsdurchfiihrung sind Filterkoeffizienten zu quantisieren. Geben Sie die MATLAB-Befehlszeile an, die eine Zahl x mit -1 < jc < 1-LSB auf den Dezimalwert Xq mit Runden umrechnet. Verwenden Sie das Zweierkomplement-Format mit der Wortlange w in Bit.

178

A15.3-2

15 Reale digitale Filter: Koeffizientenquantisierung

Geben Sie fiir die Filterordnung N=20 die obere Schranke (15.3) ftir die Betragsfrequenzgangsabweichung aufgrund der Koeffizientenquantisierung in Tabelle 15-2 an. 1st die Abschatzungen beim Entwurf selektiver Filter mit vorgegebenem Toleranzschema hilfreich? Tabelle 15-2 Abschatzung der maximalen Frequenzgangsabweichung 8

Wortlange w in Bit

16

Quantisierungsintervallbreite Q max. Frequenzgangsabweichung (15.3)

A15.3-3

tJberlegen Sie, ob durch die Koeffizientenquantisierung die lineare Phase eines FIRFilters verloren geht?

15.3.3

Versuchsdurchfuhrung

n

Zum Toleranzschema in Bild 15-3 wurde ein Tiefpass mit dem MATLAB Filterdesign-und-analyse-Werkzeug f d a t o o l entworfen. Die ersten elf Koeffizienten der Impulsantwort sind in Tabelle 15-3 aufgelistet. Beachten Sie die gerade Symmetrie der Impulsantwort. Hinweise: (i) Der Einfachheit halber wurde auf die Skaherung der Koeffizienten (15.4) verzichtet und soil auch im Weiteren nicht vorgenommen werden. (ii) Die Quantisierung geschieht im Zweierkomplement-Format, siehe Programm q a n t 2 c . m . (iii) Verwenden Sie das MATLAB-Programm f d a t o o l zur Darstellung des Frequenzgangs. ISD

H{eJ^)

= 0.04

iiii 0.8

••

0.6 I-

0.4 L

illi mmmm mWy9M 5 mmm4

Wmm mmm iW jg-

Q,o/n = 0.2

,2