Digitalfilter mit ATmega

Auch mit 8-Bit-Mikrocontrollern der ATmega-Reihe ist es möglich, einfache digitale Signalverarbeitungen auszuführen. Mittels kostenloser Mathematiksoftware ist die Berechnung eigener Filter möglich.

Viele von uns haben schon analoge Filterschaltungen aufgebaut, als passive LC-Filter oder aktive Filter mit Operationsverstärkern. Dank kostenloser Berechnungssoftware, wie z.B. AADE, ist die Auswahl der Schaltung und Berechnung der Bauteilwerte kein Problem mehr. Die gleichen Filter in digitaler Technik zu erstellen gilt dagegen immer noch als kompliziert, obwohl sie in professionellen Schaltungen die Analogtechnik weitgehend verdrängt haben. Ihr großer Vorteil ist der fehlende Abgleich. Dank Quarzoszillator wird ein einmal berechneter Frequenzgang ohne Präzisionsbauteile reproduzierbar eingehalten.

Die Software ist im Grunde sehr einfach aufgebaut, besteht nur aus Schiebeoperationen, Additionen und einigen Multiplikationen. Die ATmega-Controller enthalten alle dazu nötigen Funktionen, einen A/D-Wandler, einen 8-Bit Multiplizierer und mehrere D/A-Wandler durch Pulsbreitenmodulation. Am Ausgang muß nur noch ein einfacher Tiefpaß die PWM-Frequenz wegfiltern. Die Programmbeispiele wurden auf Grundlage der Applikationsschrift AVR223 mit AVRStudio entwickelt.

Am Ende der Programme befindet sich eine Koeffiziententabelle, deren Werte die Filtereigenschaften bestimmen. Um eigene Filter entwickeln zu können sollen dazu die Grundbegriffe sehr kurz erläutert und dann einfache „Kochrezepte“ vorgestellt werden.

Ein IIR-Filter 6.Ordnung

Als erstes wird ein Filter berechnet, das nur durch Änderung der Koeffizienten als Tiefpaß, Hochpaß, Bandpaß oder Bandsperre wirkt. Im Beispielprogramm soll es ein 1750 Hz Sperrfilter werden, das laute Ruftöne beim Zuhören auf Amateur-Relaisfunkstellen unterdrücken kann. Dieser Filtertyp heißt „IIR“, infinite impulse response, weil er nach Anregung mit einem Impuls und schlecht gewählten Einstellungen zum Oszillator wird. Der verwendete Aufbau aus drei hintereinandergeschalteten Teilfiltern zweiter Ordnung ist als gutmütig bekannt. Diese kaskadierten Teilfilter nennt man „SOS“ second order sections oder auch biquad sections. Der Aufbau mit einem einzigen Filter 6. Ordnung, wie er in der AVR223 von 2001 noch beschrieben ist, schwingt dagegen sehr zuverlässig.

Software zur Koeffizientenberechnung

Es gibt einige spezialisierte Programme, die aber nur eine kleine Anzahl von Digitalfilterkonfigurationen kennen, z.B. FilterFree. Wesentlich universeller sind umfangreiche Programmpakete für numerische Mathematik. Neben dem kommerziellenMatlab das auch in einer eingeschränkten Studentenversion angeboten wird, und zu dem zahlreiche Literatur existiert, z.B. ISBN 9783486584271 , gibt es auch zwei kostenlose Nachahmer: Octave und Scilab. Spezialitäten der digitalen Signalverarbeitung bilden nur einen kleinen Teil des Gesamtpakets, und die Hilfetexte dazu sind sehr knapp oder auch mal ein Jahrzehnt alt. Zum Glück sind beide Programme dem Vorbild so nahe, daß man vieles der Literatur zu Matlab entnehmen kann. Die eigentliche Koeffizientenberechnung besteht meistens aus nur wenigen Textzeilen. Man gibt wie beim analogen Filter die gewünschten Durchlaß- und Sperrfrequenzen, die Filtercharakteristik und Ordnung ein. Auf eine grafische Benutzeroberfläche muß man allerdings verzichten. Immerhin können der Frequenzgang und andere Kurven grafisch ausgegeben werden.

Transformationen und komplexe Zahlen

Die mathematischen Grundbegriffe der Digitalfilterberechnung sollen hier nur stichwortartig erwähnt werden:

Von der Fourier-Transformation hat man schon mal gehört, ihre Bedeutung ist noch einigermaßen begreifbar. Von einem periodischen Signal x(t) wie es das Oszilloskop zeigt, also im Zeitbereich, gelangt man über sie zur spektralen Darstellung x(f) im Frequenzbereich auf dem Spektrumanalysator. Das ist aber nur die halbe Information, eigentlich gehört zu jeder Harmonischen außer der angezeigten Amplitude auch noch ihre Phasenlage. x(f) ist also eine zweidimensionale Variable, Amplitude (Betrag) und Phase (Winkel) sind polare Koordinaten. Wenn wir zu kartesischen Koordinaten übergehen, und die x- Achse als reelle, die y-Achse als imaginäre Achse bezeichnen, kommen wir zur komplexen Zahlendarstellung. Für die Berechnung analoger Filter ist die Laplace-Transformation zuständig, für digitale Filter die nahe verwandte z-Transformation. Im ersten Fall wird x(t) in ein X(s) transformiert, in unserem Fall, wie der Name schon sagt in ein X(z). Zur besseren Unterscheidung wird das X groß geschrieben. Auch hier rechnen wir mit komplexen Zahlen, wir befinden uns nach der Transformation in der komplexen z-Ebene.

Übertragungsfunktion und Koeffizienten

Das Digitalfilter ist für uns zunächst ein schwarzer Kasten, in den ein Signal x(t) hineinläuft und als gefiltertes Signal y(t) wieder herauskommt. Die z-transformierten Signale heißen X(z) und Y(z). Der Quotient Y(z)/X(z) = H(z) wird Übertragungsfunktion, engl. transfer function genannt. Für bestimmte Werte von z wird der Zähler oder Nenner der Übertragungsfunktion Null, diese heißen Nullstellen und Pole der Funktion. Diese Werte der komplexen Zahl z können im Pol/Nullstellendiagramm als Punkte in der z-Ebene dargestellt werden, das hilft beispielsweise bei der Untersuchung der Stabilität des Filters gegen Schwingneigung.

Die genannten Mathematikprogramme berechnen uns aus den gewünschten Filterdaten die Übertragungsfunktion. Im einfachsten Fall können wir daraus die Zahlenwerte direkt in die Filtersoftware als Koeffiziententabelle eintragen. Die Aufteilung in Teilfilter macht die Berechnung etwas umständlicher. Die einzelnen Übertragungsfunktionen werden multipliziert, also Y(z)/X(z) = H₁(z)*H₂(z)*H₃(z). Wir müssen also H(z) in geeignete Faktoren zerlegen, auch dabei hilft uns die Mathematiksoftware.

Octave und Scilab

Am besten installieren wir gleich beide Softwarepakete. Beide bieten unterschiedliche Lösungswege zur Berechnung der gewünschten SOS-Teilfilter. In Scilab werden Frequenzen als Bruchteile der Abtastfrequenz, engl. sample rate des Analog-Digitalwandlers angegeben. Der AD-Wandler im ATmega wird hier mit einer Taktfrequenz von 1/64 des Quarzoszillators betrieben und braucht pro Wandlung 13 Takte. Mit einem 15 MHz – Quarz ergibt sich also ein Abtastfrequenz von 15MHz/(64*13) = 18029 Hz. Octave rechnet mit Bruchteilen der halben Abtastfrequenz, auch Nyquist-Frequenz genannt. Die Filterordnung bezieht sich in beiden Programmen auf Tief- und Hochpassfilter. Für Bandpässe und -sperren muß die Hälfte angegeben werden, hier also 3 statt 6.

Beginnen wir mit Scilab. Das kurze Programm in Bild 1 speichern wir als Textdatei z.B. unter „Sperrfilter1750Hz.sci“ ab, dann lädt es durch Anklicken direkt in den bunten Scilab-Editor und wird mit „Ausführen – in Scilab laden“ gestartet. Ein neues Fenster mit der Frequenzgangkurve Bild 2 öffnet sich und im Scilab-Hauptfenster stehen die Übertragungsfunktionen der drei Teilfilter, sowie die Gesamtverstärkung gain. Die Filterordnung kann nicht direkt vorgegeben werden, wir müssen die vier Grenzfrequenzen und die beiden Dämpfungswerte variieren bis genau drei SOS-Teilfilter herauskommen. In Octave werden nur zwei Grenzfrequenzen verlangt, dafür können wir die Filterordnung vorgeben. Wieder wird eine Textdatei ausgeführt und liefert Übertragungsfunktionen für drei Teilfilter.

Einen Schönheitsfehler haben die berechneten Funktionen noch. In der Regelungstechnik und auch in Scilab/Octave wird mit positiven Exponenten von z gerechnet, die Literatur zur digitalen Signalverarbeitung benutzt dagegen negative. Wenn wir aber Zähler und Nenner jeder Teil-Übertragungsfunktion durch z² teilen, ändert sie sich nicht und wir haben die Formel wie sie für digitale Filter üblich ist.

Jetzt können wir die Koeffizienten direkt aus den Übertragungsfunktionen ablesen. Wir müssen sie noch skalieren, also mit einer geeigneten Zweierpotenz multiplizieren und zu einer ganzen Zahl auf/abrunden. Sie sollen hier zu einer 16 Bit-Zahl (Festkommadarstellung mit Vorzeichenbit) werden, also alle Koeffizienten im Bereich -32768...+32767 liegen.

Frequenzgang maßgeschneidert

Mit der Funktion yulewalk können wir auch kompliziertere Frequenzgangkurven erreichen. In Octave fehlt sie, drei ähnlich lautende Funktionen haben nichts mit Filtern zu tun. Die gewünschte Funktion wird mit Geradenstücken zwischen beliebigen Kurvenpunkten vorgegeben. Als Beispiel nehmen wir einen Bandpass für SSTV slow scan television, der 1974 in der "CQ-DL" vorgestellt wurde. Synchronimpuls (1200 Hz) und Bildinformation (1500-2300 Hz) sollen durchgelassen, Störungen außerhalb dieser Bereiche unterdrückt werden. Nach längerem Ausprobieren verschiedener Vorgabewerte und Reduktion der Abtastfrequenz auf die Hälfte ergibt sich eine gute Übereinstimmung mit dem Frequenzgang des Originals. Jetzt haben wir eine Übertragungsfunktion 6.Ordnung, die wieder in SOS-Teilfilter zerlegt werden muß. Diese Funktion tf2sos transfer function to second order sections fehlt in Scilab, deshalb übertragen wir die Zahlen aus Scilab nach Octave und erhalten so schließlich unsere Koeffizienten.

Literatur

Der wichtigste Artikel zur Signalverarbeitung mit Scilab ist Signal.pdf (1,2MByte) von 1998. Die dazugehörigen Programme sind etwas schwerer zu finden, hier habe ich sie daher zu SignalFiles.pdf (980kB) zusammengefasst.

Scilab-Berechnung mit eqiir

<c> // Stopband-Filter für 1750 Hz, Samplerate =18029Hz // Eckfrequenzen 1400Hz, 1700Hz, 1800 Hz, 2100Hz: omega=[2*%pi*(1400/18029),2*%pi*(1700/18029),2*%pi*(1800/18029),2*%pi*(2100/18029)]; // Maximale Durchgangswelligkeit 0,5 dB deltapass = (1-10**(-0.5/20)) ; // minimale Sperrband-Dämpfung 50 dB deltastop = 10**(-50/20); // IIR-Filter berechnen, "stopband", "elliptic": [sos,gain,zeroes,poles] = eqiir('sb','el',omega,deltapass,deltastop); // Second order sections anzeigen: sos gain //Frequenzgang berechnen und plotten: zaehlerprodukt = prod(sos(2)); nennerprodukt = prod(sos(3)); frequenzgang = gain*abs(freq(zaehlerprodukt,nennerprodukt,exp(%i*(0:0.01:%pi)))); // frequenzgang = gain*abs(freq(zaehlerprodukt,nennerprodukt,exp(%i*(0:0.01:%pi)))); n = prod(size(frequenzgang)); //plot(20*log(frequenzgang(2:n))/log(10),(n*18029/(200*%pi))); plot(20*log(frequenzgang (2:n))/log(10))

(Ergebnis:)

                           2                             2                           2
         1 - 1.6219938z + z            1 - 1.6396794z + z           1 - 1.656616z + z

sos = -------------------------- * ------------------------ * -------------------------

                               2                             2                            2
     0.9383702 - 1.5100165z + z     0.8466245 - 1.513936z + z    0.9485511 - 1.663563z + z

gain = 0.8683583 </c>

Messung des 1750 Hz-Sperrfilters

Gemessen mit "Visual Analyzer":

Frequenzgang-Planung mit Yulewalk

Oben das Programm, weitgehend aus dem Hilfetext entnommen, in der Mitte der berechnete Frequenzgang - Messung liegt noch keine vor -, unten rechts die Vorgabe mit Geradenstücken, links die Übertragungsfunktion 6.Ordnung.

Die Koeffizienten heißen im Zähler von links nach rechts b6...b0 , im Nenner a6...a0, wobei a0 immer =1 ist. Hier die Zerlegung in Teilfilter mit Octave: <c> B=[0.1211431 0.0379961 -0.0272326 -0.0361440 -0.0041574 -0.0010921 -0.0894734 ]; A=[1 -1.7641816 2.7032943 -2.3805502 1.8419165 -0.8344276 0.2958684]; [sos,g] = tf2sos(B,A) // // octave.exe:2> B=[0.1211431 0.0379961 -0.0272326 -0.0361440 -0.0041574 -0.0010921 -0.0894734 ]; // octave.exe:3> A=[1 -1.7641816 2.7032943 -2.3805502 1.8419165 -0.8344276 0.2958684]; // octave.exe:4> [sos,g] = tf2sos(B,A) // sos = // // 1.0000000 1.1503168 1.0001280 1.0000000 0.0596667 0.7154071 // 1.0000000 -0.8390889 0.7389315 1.0000000 -0.7453174 0.7120574 // 1.0000000 0.0024186 -0.9993911 1.0000000 -1.0785309 0.5808050 // // g = 0.12114 </c> In drei Zeilen untereinander die drei Teilfilter-Koeffizienten, von links nach rechts b0, b1, b2, a0 (wieder=1), a1 und a2.

Quelltext des ATmega48-Programms

basierend auf AVR223, zerlegt in 3 SOS-Teilfilter, mit Überlaufbegrenzung <avrasm>

• Digital filter with ATmega48P 6th-order IIR, AVR223 application note *
• 15 MHz xtal, input ADC0, PWM-output OCR1a/b, test output PD0 *
• Christoph Kessler 2008 db1uq_at_t-online_dot_de *

.nolist .include "m48pdef.inc" .list

register 0-15 (no "immediate"-operations possible)

.def MltLo = r0 ; multiplier result .def MltHi = r1 ; multiplier result .def Zero = r2 ; Zero register (used to add carry flag) .def SavSrg = r3 ; save status register during interrupt routine .def reg04 = r4 ; free .def reg05 = r5 ; free .def reg06 = r6 ; free .def reg07 = r7 ; free .def reg08 = r8 ; free .def reg09 = r9 ; free .def reg10 = r10 ; free .def reg11 = r11 ; free .def reg12 = r12 ; free .def accu0 = r13 ; Accumulator - 32bit for Interrupt routine .def accu1 = r14 ; Accumulator - 32bit for Interrupt routine .def accu2 = r15 ; Accumulator - 32bit for Interrupt routine

register 16-23 ("immediate"-operations possible)

.def accu3 = r16 ; Accumulator - 32bit for Interrupt routine .def reg17 = r17 ; free .def reg18 = r18 ; free .def reg19 = r19 ; free .def datal = r20 ; Data samples mul register .def datah = r21 ; for Interrupt routine .def coefl = r22 ; Filter coefficient mul register .def coefh = r23 ; for Interrupt routine

• register 24-31 (16Bit-registers)

.def Tmp1 = r24 ; general temporary register .def Tmp2 = r25 ; general temporary register

XL = r26 ; free

XH = r27 ; free

YL = r28 ; used by Interrupt routine

YH = r29 ; used by Interrupt routine

ZL = r30 ; init routine only, free for main program

ZH = r31 ; init routine only, free for main program

constants

1st SOS filter

.equ x11 = 0 .equ x21 = 2 .equ y11 = 4 .equ y21 = 6

2nd SOS filter

.equ x12 = 8 .equ x22 = 10 .equ y12 = 12 .equ y22 = 14

3rd SOS filter

.equ x13 = 16 .equ x23 = 18 .equ y13 = 20 .equ y23 = 22

1st SOS filter

.equ b01 = 24 .equ b11 = 26 .equ b21 = 28 .equ a11 = 30 .equ a21 = 32

2nd SOS filter

.equ b02 = 34 .equ b12 = 36 .equ b22 = 38 .equ a12 = 40 .equ a22 = 42

3rd SOS filter

.equ b03 = 44 .equ b13 = 46 .equ b23 = 48 .equ a13 = 50 .equ a23 = 52

• Macros

.MACRO load_node ldd datal,Y+@0 ; Load low byte of node ldd datah,Y+@0+1 ; Load high byte of node .ENDMACRO .MACRO update_node std Y+@0,datal ; Update low byte of node to prepare next filter run std Y+@0+1,datah ; Update high byte of node to prepare next filter run .ENDMACRO .MACRO load_coef ldd coefl,Y+@0 ; Load low byte of node ldd coefh,Y+@0+1 ; Load high byte of node .ENDMACRO .MACRO mul_move_32 muls coefh,datah ; Signed multiply, coefficient high byte and data high byte mov accu2,MltLo ; Copy result word into accumulator byte 2:3 mov accu3,MltHi ; Copy result word into accumulator byte 2:3 mul coefl,datal ; Unsigned multiply, coefficient low byte and data low byte mov accu0,MltLo ; Copy result word into accumulator byte 2:3 mov accu1,MltHi ; Copy result word into accumulator byte 2:3 mulsu coefh,datal ; Signed-unsigned multiply, coefficient high byte and data low byte sbc accu3,Zero ; Sign extention add accu1,MltLo ; Add low byte of result to accumulator byte 1 adc accu2,MltHi ; Add with carry high byte of result to accumulator byte 2 adc accu3,Zero ; Add carry to accumulator byte 3 mulsu datah,coefl ; Signed-unsigned multiply, data high byte and coefficient low byte sbc accu3,Zero ; Sign extention add accu1,MltLo ; Add low byte of result to accumulator byte 1 adc accu2,MltHi ; Add with carry high byte of result to accumulator byte 2 adc accu3,Zero ; Add carry to accumulator byte 3 .ENDMACRO .MACRO smac32 muls coefh,datah ; Signed multiply, coefficient high byte and data high byte add accu2,MltLo ; Add low byte of result to accumulator byte 2 adc accu3,MltHi ; Add with carry high byte of result to accumulator byte 3 mul coefl,datal ; Unsigned multiply, coefficient low byte and data low byte add accu0,MltLo ; Add low byte of result to accumulator byte 0 adc accu1,MltHi ; Add with carry high byte of result to accumulator byte 1 adc accu2,Zero ; Add carry to accumulator byte 2 adc accu3,Zero ; Add carry to accumulator byte 3 mulsu coefh,datal ; Signed-unsigned multiply, coefficient high byte and data low byte sbc accu3,Zero ; Sign extention add accu1,MltLo ; Add low byte of result to accumulator byte 1 adc accu2,MltHi ; Add with carry high byte of result to accumulator byte 2 adc accu3,Zero ; Add carry to accumulator byte 3 mulsu datah,coefl ; Signed-unsigned multiply, data high byte and coefficient low byte sbc accu3,Zero ; Sign extention add accu1,MltLo ; Add low byte of result to accumulator byte 1 adc accu2,MltHi ; Add with carry high byte of result to accumulator byte 2 adc accu3,Zero ; Add carry to accumulator byte 3

ret

.ENDMACRO .LISTMAC

reset/interrupt-vectors

.cseg rjmp Initial ; after RESET to main program reti ; External Interrupt Request 0 reti ; External Interrupt Request 1

reti ; External Interrupt Request 2 mega644

reti ; Pin Change Interrupt Request 0 reti ; Pin Change Interrupt Request 1 reti ; Pin Change Interrupt Request 2

reti ; Pin Change Interrupt Request 3 mega644

reti ; Watchdog Time-out Interrupt reti ; Timer/Counter2 Compare Match A reti ; Timer/Counter2 Compare Match B reti ; Timer/Counter2 Overflow reti ; Timer/Counter1 Capture Event reti ; Timer/Counter1 Compare Match A reti ; Timer/Counter1 Compare Match B reti ; Timer/Counter1 Overflow reti ; Timer/Counter0 Compare Match A reti ; Timer/Counter0 Compare Match B reti ; Timer/Counter0 Overflow reti ; SPI Serial Transfer Complete reti ; USART0, Rx Complete reti ; USART0 Data register Empty reti ; USART0, Tx Complete

reti ; Analog Comparator mega644

rjmp ADCint ; ADC Conversion Complete reti ; EEPROM Ready reti ; 2-wire Serial Interface reti ; Store Program Memory Read

after reset

• initialise stack,ports

Initial: ldi Tmp1,Low(RAMEND) out SPL,Tmp1 ; ldi Tmp1,High(RAMEND) out SPH,Tmp1 ; stack initialised clr Zero ; always Zero ldi Tmp1,$FF ; out PortB,Tmp1 ; PortB = % 1111 1111 out PortC,Zero ; PortC = % 0000 0000 out PortD,Tmp1 ; PortD = % 1111 1111 out DDRB,Tmp1 ; PortB direction = % 1111 1111 out DDRC,Zero ; PortC direction = % 0000 0000 out DDRD,Tmp1 ; PortD direction = % 1111 1111 ldi Tmp1,$A1 ; pos. outputs, fast PWM 8Bit sts TCCR1A,Tmp1 ; ldi Tmp1,$09 ; fast PWM 8Bit, clk/1 (no prescaling) sts TCCR1B,Tmp1 ; sts OCR1AH,Zero ; sts OCR1AL,Zero ; sts OCR1BH,Zero ; sts OCR1BL,Zero ; ldi Tmp1,$40 ; Ref=Vcc, right-adj, ADC0=input sts ADMUX,Tmp1 ;

ldi Tmp1,$EE ; Start,Auto, enable Int, Clk/64

ldi Tmp1,$EF ; Start,Auto, enable Int, Clk/128 sts ADCSRA,Tmp1 ; ldi Tmp1,$00 ; free running sts ADCSRB,Tmp1 ; ldi Tmp1,$01 ; disable ADC0 dig.inp (optional) sts DIDR0,Tmp1 ; ldi ZL,low(CoTblF<<1) ; move filter coefficents ldi ZH,high(CoTblF<<1) ; from program-flash ldi YL,low(CoTblS) ; to SRAM ldi YH,high(CoTblS) ; ldi Tmp1,30 ; count 15 * 16 Bit coefficients TbLoop: lpm Tmp2,Z+ ; fetch 1 byte from coefficient table st Y+,Tmp2 ; copy coefficient to Sram, increment Y dec Tmp1 brne TbLoop ldi YL,low(DatTbl) ; load SRAM-Pointer for Interrupt routine ldi YH,high(DatTbl) ; sei ; enable interrupts

Main program

Endless: rjmp Endless ;

ADC-Interrupt routine

• compute filter and update PWM output

ADCint: sbi PortD,0 ; just to measure INT-Time in SavSrg,SREG ; save status register

• 1st SOS filter

load_coef b21 ; b21*x1[t-2] load_node x21 ; mul_move_32 ; load_coef b11 ; b11*x1[t-1] load_node x11 ; update_node x21 ; smac32 ; load_coef b01 ; b01*x1[t] lds datal,ADCL ; Load new sample lds datah,ADCH ; into data multiply register dec datah ; convert 10 bit unsigned $03FF -> $01FF dec datah ; to 16 Bit signed $0000 -> $FE00

lsl datal rol datah lsl datal rol datah lsl datal rol datah lsl datal rol datah lsl datal rol datah lsl datal rol datah

update_node x11 ; smac32 ; load_coef a21 ; a21*y1[t-2] load_node y21 ; smac32 ; load_coef a11 ; a11*y1[t-1] load_node y11 ; update_node y21 ; smac32 ; cpi accu3,$20 ; accu < $20000000 ? brlo NoLimit1 ; cpi accu3,-$20 ; accu >= -$20000000 ? brsh NoLimit1 ; clr accu2 ; accu2 = $00 cpi accu3,0 ; brge PosLim1 ; ldi accu3,$80 ; accu3 = $80 rjmp Limit1 ; PosLim1: ldi accu3,$7F ; accu3 = $7F com accu2 ; accu2 = $FF rjmp Limit1 ; NoLimit1: lsl accu1 ; Scaling to gain factor 2^15 rol accu2 ; rol accu3 ; lsl accu1 ; Scaling to gain factor 2^16 rol accu2 ; rol accu3 ; Limit1: std Y+y11,accu2 ; Updating y1[t-1] node to prepare next filter run std Y+y11+1,accu3 ; and save it for following SOS filter

• 2nd SOS filter

load_coef b22 ; b22*x2[t-2] load_node x22 ; mul_move_32 ; load_coef b12 ; b12*x2[t-1] load_node x12 ; update_node x22 ; smac32 ; load_coef b02 ; b02*x2[t] load_node y11 ; result of 1st SOS update_node x12 ; smac32 ; load_coef a22 ; a22*y2[t-2] load_node y22 ; smac32 ; load_coef a12 ; a12*y2[t-1] load_node y12 ; update_node y22 ; smac32 ; cpi accu3,$20 ; accu < $20000000 ? brlo NoLimit2 ; cpi accu3,-$20 ; accu >= -$20000000 ? brsh NoLimit2 ; clr accu2 ; accu2 = $00 cpi accu3,0 ; brge PosLim2 ; ldi accu3,$80 ; accu3 = $80 rjmp Limit2 ; PosLim2: ldi accu3,$7F ; accu3 = $7F com accu2 ; accu2 = $FF rjmp Limit2 ; NoLimit2: lsl accu1 ; Scaling to gain factor 2^15 rol accu2 ; rol accu3 ; lsl accu1 ; Scaling to gain factor 2^16 rol accu2 ; rol accu3 ; Limit2: std Y+y12,accu2 ; Updating y2[t-1] node to prepare next filter run std Y+y12+1,accu3 ; and save it for following SOS filter

• 3rd SOS filter

load_coef b23 ; b23*x3[t-2] load_node x23 ; mul_move_32 ; load_coef b13 ; b13*x3[t-1] load_node x13 ; update_node x23 ; smac32 ; load_coef b03 ; b03*x3[t] load_node y12 ; result of 2nd SOS update_node x13 ; smac32 ; load_coef a23 ; a23*y3[t-2] load_node y23 ; smac32 ; load_coef a13 ; a13*y3[t-1] load_node y13 ; update_node y23 ; smac32 ; cpi accu3,$20 ; accu < $20000000 ? brlo NoLimit3 ; cpi accu3,-$20 ; accu >= -$20000000 ? brsh NoLimit3 ; clr accu2 ; accu2 = $00 cpi accu3,0 ; brge PosLim3 ; ldi accu3,$80 ; accu3 = $80 rjmp Limit3 ; PosLim3: ldi accu3,$7F ; accu3 = $7F com accu2 ; accu2 = $FF rjmp Limit3 ; NoLimit3: lsl accu1 ; Scaling to gain factor 2^15 rol accu2 ; rol accu3 ; lsl accu1 ; Scaling to gain factor 2^16 rol accu2 ; rol accu3 ; Limit3: std Y+y13,accu2 ; Updating y3[t-1] node to prepare next filter run std Y+y13+1,accu3 ;

• PWM - D/A-output

subi accu3,$80 ; signed to unsigned $8000 ->$0000, $7FFF ->$FFFF sts OCR1AH,Zero ; sts OCR1AL,Tmp1 ; fast PWM 8 bit input signal to OCR1A sts OCR1BH,Zero ; sts OCR1BL,accu3 ; fast PWM 8 bit result to OCR1B cbi PortD,0 ; just to measure INT-Time out SREG,SavSrg ; restore status register reti ;

• coefficient table in flash memory

CoTblF: /*

1750 Hz Sperrfilter

.dw 16384 ; b01 *16384 ( 1 ) .dw -26575 ; b11 *16384 (-1,6219938) .dw 16384 ; b21 *16384 ( 1 ) .dw 24740 ; -a11 *16384 ( 1,5100165) .dw -15374 ; -a21 *16384 (-0,9383702)

.dw 16384 ; b02 *16384 ( 1 ) .dw -26865 ; b12 *16384 (-1,6396794) .dw 16384 ; b22 *16384 ( 1 ) .dw 24804 ; -a12 *16384 ( 1,5139300) .dw -13871 ; -a22 *16384 (-0,8466245)

.dw 16384 ; b03 *16384 ( 1 ) .dw -27142 ; b13 *16384 (-1,6566160) .dw 16384 ; b23 *16384 ( 1 ) .dw 27256 ; -a13 *16384 ( 1,6635630) .dw -15541 ; -a23 *16384 (-0,9485511) */

SSTV-Bandpass nach DJ6HP 015

.dw 16384 ; b01 *16384 ( 1 ) .dw 18847 ; b11 *16384 ( 1.1503168) .dw 16386 ; b21 *16384 ( 1.0001280) .dw -978 ; -a11 *16384 (-0.0596667) .dw -11721 ; -a21 *16384 (-0.7154071)

.dw 16384 ; b02 *16384 ( 1 ) .dw -13748 ; b12 *16384 (-0.8390889) .dw 12107 ; b22 *16384 ( 0.7389315) .dw 12211 ; -a12 *16384 ( 0.7453174 .dw -11666 ; -a22 *16384 (-0.7120574)

.dw 16384 ; b03 *16384 ( 1 ) .dw 40 ; b13 *16384 ( 0.0024186) .dw -16374 ; b23 *16384 (-0.9993911) .dw 17671 ; -a13 *16384 ( 1.0785309) .dw -9516 ; -a23 *16384 (-0.5808050)

• data memory in SRAM

.dseg .org $0100 DatTbl: .dw $0000 ; x11 .dw $0000 ; x21 .dw $0000 ; y11 .dw $0000 ; y21

.dw $0000 ; x12 .dw $0000 ; x22 .dw $0000 ; y12 .dw $0000 ; y22

.dw $0000 ; x13 .dw $0000 ; x23 .dw $0000 ; y13 .dw $0000 ; y23

• coefficient table in SRAM

CoTblS: .dw $0000 ; b01 *16384 .dw $0000 ; b11 *16384 .dw $0000 ; b21 *16384 .dw $0000 ; -a11 *16384 .dw $0000 ; -a21 *16384

.dw $0000 ; b02 *16384 .dw $0000 ; b12 *16384 .dw $0000 ; b22 *16384 .dw $0000 ; -a12 *16384 .dw $0000 ; -a22 *16384

.dw $0000 ; b03 *16384 .dw $0000 ; b13 *16384 .dw $0000 ; b23 *16384 .dw $0000 ; -a13 *16384 .dw $0000 ; -a23 *16384

</avrasm>

Spezielle Digitalfilter

Neben den üblichen Filtern, Tiefpass, Hochpass, Bandpass und Bandsperre, gibt es einige weitere Typen. In der Literatur findet man u.a. noch den Allpass, Phasenschieber, Integrator und Differentiator.

Der Hilbert-Transformator, ein Breitband-Phasenschieber

In den Lehrbüchern wird er völlig zu Unrecht eher am Rande abgehandelt. Tatsächlich findet man ihn in allen Schaltungen zur Datenübertragung, zu Modulation oder Demodulation mit I/Q-Mischern, im „Software-defined-radio“.

Auch ein Hilbert-Transformator kann mithilfe von Scilab und Octave berechnet werden. Allerdings gibt es keine fertige Software, man muß sich durch einige Literatur hindurchfressen, z.B. ISBN 0471619957 Sanjit Mitra, Handbook for DSP (1998), Kapitel „Special Filter Designs“ von Phillip Regalia.

Eine Fundgrube sind die Matlab-Texte von Prof. Schüssler, die vorab 1998 zum zweiten Band seines Lehrbuchs auf seiner Webseite veröffentlicht wurden. Inzwischen ist der erste Band in der fünften Auflage erschienen ISBN 9783540782506 , der zweite Band wird immer noch angekündigt, aber die Webseite ist nach seinem Tod 2007 abgeschaltet. Zum Glück gibt es das Webarchiv,wo die Texte noch zu finden sind.

Zu dieser Zeit waren die Matlab-Programme noch nicht in „Toolboxen“ untergebracht, deren Funktion von außen nicht mehr sichtbar ist. Daher lassen sie sich relativ einfach in Scilab oder Octave übersetzen.

Drei Versionen von Hilbert-Transformatoren werden hier berechnet:

• FIR-Hilbert, aufwändig, aber streng linearphasig, der bekannteste Typ, da schon lange mit dem Matlab-Befehl remez / firpm und Parameter 'hilbert' berechenbar. Scilab und Octave kennen diesen Parameter nicht.
• näherungsweise linearphasiger IIR-Hilbert, damit für Digitalsignale besser geeignet als der folgende:
• minimalphasiger IIR-Hilbert, kompakt, aber mit großen Gruppenlaufzeitschwankungen, eher für Audioanwendungen.

Zumindest der kompakte IIR-Hilbert läßt sich auch in einem ATmega unterbringen. Er besteht aus einer Verzögerung und zwei Allpässen, die wieder aus SOS-Teilfiltern aufgebaut sind. Pro Teilfilter gibt es nur einen Koeffizienten, sodaß man mindestens acht Teilfilter verwenden kann.

Ein Hilbert-Kochrezept (noch unvollständig)

Folgen wir dem Zahlenbeispiel und Berechnungsweg im „Handbook for DSP“.

Zuerst berechnen wir einen elliptischen Halbband-Tiefpass. Halbband besagt, dass die beiden Grenzfrequenzen symmetrisch zur halben Nyquist- Frequenz also 0,25 * Abtastfrequenz liegen, siehe Bild oben rechts. Die genauen Bedingungen lauten in der Schreibweise des Buches (p=passband, s=stopband):

• Omega_s + Omega_p = Pi
  (auf Omega_sample normierte Skala: Samplerate=2*Pi)
• (1-Delta₁)² +Delta₂²=1

Als Zahlenwerte werden Omega_s = 0,55*Pi und Delta₂=0,01 gewählt, und damit ein elliptischer IIR-Tiefpass 7.Ordnung berechnet. Damit erhält er folgende Zahlenwerte der Pole:

0  0,436688  0,743707 0,927758

Mit Scilab erhalten wir nur fast dieselben Zahlen, siehe Bild unten links (Die Differenz ist hoffentlich auf die unterschiedlichen IIR-Berechnungsverfahren zurückzuführen, nicht auf meinen Rechenfehler). Die Pole müssten für ein ideales Halbbandfilter genau auf der imaginären Achse liegen, der Realteil genau Null sein. Diese Zahlen werden noch quadriert (konjugiert komplexe Pole zusammengefasst) und bilden schließlich die Koeffizienten des Hilbert-Transformators.

Scilab-Berechnung des Halbbandfilters

Die Übertragungsfunktion des gesuchten Hilbert-Transformators soll wie bisher H(z) heißen. Zur Unterscheidung nennen wir die bis jetzt berechnete Tiefpassfunktion G(z).

Wieder wird sie in SOS-Teilfilter aufgeteilt. Hier allerdings in die Summe von zwei parallel geschaltete Allpass-Funktionen A₁ und A₂ (plus eine Verzögerung um einen Sampletakt, in der Übertragungsfunktion als „z^-1“ ausgedrückt), was mit vielen Übertragungsfunktionen möglich ist.

Ein Allpass läßt alle Signalfrequenzen mit unveränderter Amplitude durch, verändert nur die Phase. Seine SOS-Teilfunktion ist besonders einfach aufgebaut, zwei Koeffizienten sind Null, zwei sind gleich Eins und die beiden restlichen identisch. Damit benötigt man nur eine Multiplikation pro Teilfilter.

Die Tiefpass-Übertragungsfunktion erhält so die Form:

G(z) = ½ * ( A₁ (z²) + z^-1 * A₂ (z²))

Zum Übergang auf den Hilbert-Transformator ersetzen (substituieren) wir das z durch -jz und multiplizieren das ganze mit 2:

H(z)= 2*G(-jz)

H(z) = A₁ (-z^-2) + j* z^-1 * A₂ (-z^-2)

Die Ausgänge der beiden Allpässe sind jetzt gegeneinander um 90 Grad phasenverschoben, und können auf einen I/Q-Modulator gegeben werden. Umgekehrt kann man auch die beiden Eingänge auftrennen, und zwei um 90 Grad versetzte Signale aus einem I/Q-Demodulator einspeisen und die Ausgangssignale addieren oder subtrahieren, je nach gewünschtem Seitenband.

Weitere Literatur zu Hilbert

• Schüssler / Steffen "Halfband Filters and Hilbert Transformers"
  Irgendwo im Web hatte ich den vollständigen Artikel als PDF gefunden, leider die Adresse nicht notiert, Google findet nur kostenpflichtige Downloads.
• Dutta / Kumar On Digital Differentiators, Hilbert Transformers, and Half-Band Low-Pass Filters
  dank Schrifterkennungsprogramm leicht lädiert
• Miroslav Lutovac / Ljiljana Mili Half-band IIR Filter Design Using Matlab, weitere Artikel der Autoren
• Göckler / Damjanovic A Family of Efficient Complex Halfband Filters

ATmega-Programm Hilbert-Transformer ( fast fertig, schwingt leider)

mit den Koeffizienten von Olli Niemitalo Da der Adressenbereich des ATmega nicht groß genug ist, wird die Tabelle im SRAM in zwei Hälften geteilt. <avrasm>

• Hilbert transformer with ATmega48P, 90 degree broadband phase shifter *
• 15 MHz xtal, input ADC0, PWM-output OCR1a/b, test output PD0 *
• Christoph Kessler 2008 db1uq_at_t-online_dot_de *

.nolist .include "m48pdef.inc" .list

register 0-15 (no "immediate"-operations possible)

.def mult_l = r0 ; multiplier result .def mult_h = r1 ; multiplier result .def zero = r2 ; zero register (used to add carry flag) .def savsrg = r3 ; save status register during interrupt routine .def reg04 = r4 ; free .def reg05 = r5 ; free .def reg06 = r6 ; free .def reg07 = r7 ; free .def reg08 = r8 ; free .def reg09 = r9 ; free .def reg10 = r10 ; free .def reg11 = r11 ; free .def reg12 = r12 ; free .def accu0 = r13 ; Accumulator - 32bit for Interrupt routine .def accu1 = r14 ; Accumulator - 32bit for Interrupt routine .def accu2 = r15 ; Accumulator - 32bit for Interrupt routine

register 16-23 ("immediate"-operations possible)

.def accu3 = r16 ; Accumulator - 32bit for Interrupt routine .def reg17 = r17 ; free .def reg18 = r18 ; free .def reg19 = r19 ; free .def datal = r20 ; Data samples mul register .def datah = r21 ; for Interrupt routine .def coefl = r22 ; Filter coefficient mul register .def coefh = r23 ; for Interrupt routine

• register 24-31 (16Bit-registers)

.def tmp1 = r24 ; general temporary register .def tmp2 = r25 ; general temporary register

XL = r26 ; free

XH = r27 ; free

YL = r28 ; used by Interrupt routine

YH = r29 ; used by Interrupt routine

ZL = r30 ; init routine only, free for main program

ZH = r31 ; init routine only, free for main program

constants

1st table

;"allpass"0: .equ x10 = 0 ; allpass 2: .equ x12 = 2 .equ x22 = 4 .equ y02 = 6 .equ y12 = 8 .equ y22 = 10 ; allpass 4: .equ x14 = 12 .equ x24 = 14 .equ y04 = 16 .equ y14 = 18 .equ y24 = 20 ; allpass 6: .equ x16 = 22 .equ x26 = 24 .equ y06 = 26 .equ y16 = 28 .equ y26 = 30 ; allpass 8: .equ x18 = 32 .equ x28 = 34 .equ y08 = 36 .equ y18 = 38 .equ y28 = 40

.equ k2 = 42 .equ k4 = 44 .equ k6 = 46 .equ k8 = 48

2nd table

.equ k1 = 0 .equ k3 = 2 .equ k5 = 4 .equ k7 = 6 ; allpass 1: .equ x11 = 8 .equ x21 = 10 .equ y01 = 12 .equ y11 = 14 .equ y21 = 16 ; allpass 3: .equ x13 = 18 .equ x23 = 20 .equ y03 = 22 .equ y13 = 24 .equ y23 = 26 ; allpass 5: .equ x15 = 28 .equ x25 = 30 .equ y05 = 32 .equ y15 = 34 .equ y25 = 36 ; allpass 7: .equ x17 = 38 .equ x27 = 40 .equ y07 = 42 .equ y17 = 44 .equ y27 = 46

• Macros

.MACRO load_node ldd datal,Y+@0 ; Load low byte of node ldd datah,Y+@0+1 ; Load high byte of node .ENDMACRO .MACRO update_node std Y+@0,datal ; Store low byte of node std Y+@0+1,datah ; Store high byte of node .ENDMACRO .MACRO add_node ldd coefl,Y+@0 ; Load 2nd node low ldd coefh,Y+@0+1 ; Load 2nd node high add datal,coefl ; Add low bytes adc datah,coefh ; Add with carry high bytes .ENDMACRO .MACRO mult_32 ldd coefl,Y+@0 ; Load low byte of coefficient ldd coefh,Y+@0+1 ; Load high byte of coefficient muls coefh,datah ; Signed multiply, coefficient high byte and data high byte mov accu2,mult_l ; Copy result word into accumulator byte 2:3 mov accu3,mult_h ; Copy result word into accumulator byte 2:3 mul coefl,datal ; Unsigned multiply, coefficient low byte and data low byte mov accu0,mult_l ; Copy result word into accumulator byte 2:3 mov accu1,mult_h ; Copy result word into accumulator byte 2:3 mulsu coefh,datal ; Signed-unsigned multiply, coefficient high byte and data low byte sbc accu3,zero ; Sign extention add accu1,mult_l ; Add low byte of result to accumulator byte 1 adc accu2,mult_h ; Add with carry high byte of result to accumulator byte 2 adc accu3,zero ; Add carry to accumulator byte 3 mulsu datah,coefl ; Signed-unsigned multiply, data high byte and coefficient low byte sbc accu3,zero ; Sign extention add accu1,mult_l ; Add low byte of result to accumulator byte 1 adc accu2,mult_h ; Add with carry high byte of result to accumulator byte 2 adc accu3,zero ; Add carry to accumulator byte 3 .ENDMACRO .MACRO subtr_node ldd datal,Y+@0 ; Load 2nd node low ldd datah,Y+@0+1 ; Load 2nd node high sub accu0,datal ; Subtract low byte sbc accu1,datah ; Subtract with carry high byte sbc accu2,zero ; sbc accu3,zero ; .ENDMACRO .MACRO limit_store cpi accu3,$40 ; accu < $40000000 ? brlo NoLimit1 ; cpi accu3,-$40 ; accu >= -$40000000 ? brsh NoLimit1 ; clr accu2 ; accu2 = $00 cpi accu3,0 ; brge PosLim1 ; ldi accu3,$80 ; accu3 = $80 rjmp Limit1 ; PosLim1: ldi accu3,$7F ; accu3 = $7F com accu2 ; accu2 = $FF rjmp Limit1 ; NoLimit1: lsl accu1 ; Scaling to gain factor 2^16 rol accu2 ; rol accu3 ; Limit1: std Y+@0,accu2 ; std Y+@0+1,accu3 ; .ENDMACRO .LISTMAC

reset/interrupt-vectors

.cseg rjmp Initial ; after RESET to main program reti ; External Interrupt Request 0 reti ; External Interrupt Request 1 reti ; Pin Change Interrupt Request 0 reti ; Pin Change Interrupt Request 1 reti ; Pin Change Interrupt Request 2 reti ; Watchdog Time-out Interrupt reti ; Timer/Counter2 Compare Match A reti ; Timer/Counter2 Compare Match B reti ; Timer/Counter2 Overflow reti ; Timer/Counter1 Capture Event reti ; Timer/Counter1 Compare Match A reti ; Timer/Counter1 Compare Match B reti ; Timer/Counter1 Overflow reti ; Timer/Counter0 Compare Match A reti ; Timer/Counter0 Compare Match B reti ; Timer/Counter0 Overflow reti ; SPI Serial Transfer Complete reti ; USART0, Rx Complete reti ; USART0 Data register Empty reti ; USART0, Tx Complete rjmp ADCint ; ADC Conversion Complete reti ; EEPROM Ready reti ; 2-wire Serial Interface reti ; Store Program Memory Read

after reset

• initialise stack,ports

Initial: ldi tmp1,Low(RAMEND) out SPL,tmp1 ; ldi tmp1,High(RAMEND) out SPH,tmp1 ; stack initialised clr zero ; always zero ldi tmp1,$FF ; out PortB,tmp1 ; PortB = % 1111 1111 out PortC,zero ; PortC = % 0000 0000 out PortD,tmp1 ; PortD = % 1111 1111 out DDRB,tmp1 ; PortB direction = % 1111 1111 out DDRC,zero ; PortC direction = % 0000 0000 out DDRD,tmp1 ; PortD direction = % 1111 1111 ldi tmp1,$A1 ; pos. outputs, fast PWM 8Bit sts TCCR1A,tmp1 ; ldi tmp1,$09 ; fast PWM 8Bit, clk/1 (no prescaling) sts TCCR1B,tmp1 ; sts OCR1AH,zero ; sts OCR1AL,zero ; sts OCR1BH,zero ; sts OCR1BL,zero ; ldi tmp1,$40 ; Ref=Vcc, right-adj, ADC0=input sts ADMUX,tmp1 ; ldi tmp1,$EE ; Start,Auto, enable Int, Clk/64

ldi tmp1,$EF ; Start,Auto, enable Int, Clk/128

sts ADCSRA,tmp1 ; ldi tmp1,$00 ; free running sts ADCSRB,tmp1 ; ldi tmp1,$01 ; disable ADC0 dig.inp (optional) sts DIDR0,tmp1 ; ldi ZL,low(CoTblF<<1) ; move filter coefficents ldi ZH,high(CoTblF<<1) ; from program-flash ldi YL,low(CoTblS) ; to SRAM ldi YH,high(CoTblS) ; ldi tmp1,16 ; count 8 * 2-Byte coefficients TbLoop: lpm tmp2,Z+ ; fetch 1 byte from coefficient table st Y+,tmp2 ; copy coefficient to Sram, increment Y dec tmp1 brne TbLoop ldi YL,low(DatTbl1) ; load SRAM-Pointer for Interrupt routine ldi YH,high(DatTbl1) ; sei ; enable interrupts

Main program

Endless: rjmp Endless ;

ADC-Interrupt routine

• compute filter and update PWM output

ADCint: sbi PortD,0 ; just to measure INT-Time in savsrg,SREG ; save status register

allpass 2

load_node x10 ; add_node y22 ; y22 +---+ y12 +---+ mult_32 k2 ; o--<-|1/z|<-o--|1/z|<---+ subtr_node x22 ; | +---+ +---+ | limit_store y02 ; v | load_node x12 ; x10 +---+ +---+ +---+ | update_node x22 ; o-->| + |-->|*K2|--->| - |--->o---> load_node x10 ; | +---+ +---+ + +---+ y02 update_node x12 ; | - A lds datal,ADCL ; | +---+ +---+ | lds datah,ADCH ; +----|1/z|--o->|1/z|->-o dec datah ; | +---+ x12 +---+ x22 dec datah ; | lsl datal ; +-----------+ rol datah ; +---+ | lsl datal ; ADC >--|1/z|--+ rol datah ; +---+ lsl datal ; rol datah ; lsl datal ; convert 10 bit unsigned rol datah ; to 15 Bit signed lsl datal ; $03FF -> $01FF -> $3FE0 rol datah ; $0000 -> $FE00 -> $C000 lsl datal ; rol datah ; update_node x10 ; save new sample for 2nd allpass-chain load_node y12 ; update_node y22 ; load_node y02 ; update_node y12 ;

allpass 4

add_node y24 ; mult_32 k4 ; y24 +---+ y14 +---+ subtr_node x24 ; o--<-|1/z|<-o--|1/z|<---+ limit_store y04 ; | +---+ +---+ | load_node x14 ; v | update_node x24 ; y02 +---+ +---+ +---+ | load_node y02 ; ->o-->| + |-->|*K4|--->| - |--->o---> update_node x14 ; | +---+ +---+ + +---+ y04 load_node y14 ; | - A update_node y24 ; | +---+ +---+ | load_node y04 ; +----|1/z|--o->|1/z|->-o update_node y14 ; +---+ x14 +---+ x24

allpass 6

add_node y26 ; mult_32 k6 ; y26 +---+ y16 +---+ subtr_node x26 ; o--<-|1/z|<-o--|1/z|<---+ limit_store y06 ; | +---+ +---+ | load_node x16 ; v | update_node x26 ; y04 +---+ +---+ +---+ | load_node y04 ; ->o-->| + |-->|*K6|--->| - |--->o---> update_node x16 ; | +---+ +---+ + +---+ y06 load_node y16 ; | - A update_node y26 ; | +---+ +---+ | load_node y06 ; +----|1/z|--o->|1/z|->-o update_node y16 ; +---+ x16 +---+ x26

allpass 8

add_node y28 ; mult_32 k8 ; y28 +---+ y18 +---+ subtr_node x28 ; o--<-|1/z|<-o--|1/z|<---+ limit_store y08 ; | +---+ +---+ | load_node x18 ; v | update_node x28 ; y06 +---+ +---+ +---+ | load_node y06 ; ->o-->| + |-->|*K8|--->| - |--->o---> update_node x18 ; | +---+ +---+ + +---+ y08 load_node y18 ; | - A update_node y28 ; | +---+ +---+ | load_node y08 ; +----|1/z|--o->|1/z|->-o update_node y18 ; +---+ x18 +---+ x28 load_node x10 ; ldi YL,Low(DatTbl2) ;

allpass 1

add_node y21 ; mult_32 k1 ; y21 +---+ y11 +---+ subtr_node x21 ; o--<-|1/z|<-o--|1/z|<---+ limit_store y01 ; | +---+ +---+ | load_node x11 ; v | update_node x21 ; x10 +---+ +---+ +---+ | load_node x10 ; ->o-->| + |-->|*K1|--->| - |--->o---> update_node x11 ; | +---+ +---+ + +---+ y01 load_node y11 ; | - A update_node y21 ; | +---+ +---+ | load_node y01 ; +----|1/z|--o->|1/z|->-o update_node y11 ; +---+ x11 +---+ x21

allpass 3

add_node y23 ; mult_32 k3 ; y23 +---+ y13 +---+ subtr_node x23 ; o--<-|1/z|<-o--|1/z|<---+ limit_store y03 ; | +---+ +---+ | load_node x13 ; v | update_node x23 ; y01 +---+ +---+ +---+ | load_node y01 ; ->o-->| + |-->|*K3|--->| - |--->o---> update_node x13 ; | +---+ +---+ + +---+ y03 load_node y13 ; | - A update_node y23 ; | +---+ +---+ | load_node y03 ; +----|1/z|--o->|1/z|->-o update_node y13 ; +---+ x13 +---+ x23

allpass 5

add_node y25 ; mult_32 k5 ; y25 +---+ y15 +---+ subtr_node x25 ; o--<-|1/z|<-o--|1/z|<---+ limit_store y05 ; | +---+ +---+ | load_node x15 ; v | update_node x25 ; y03 +---+ +---+ +---+ | load_node y03 ; ->o-->| + |-->|*K5|--->| - |--->o---> update_node x15 ; | +---+ +---+ + +---+ y05 load_node y15 ; | - A update_node y25 ; | +---+ +---+ | load_node y05 ; +----|1/z|--o->|1/z|->-o update_node y15 ; +---+ x15 +---+ x25

allpass 7

add_node y27 ; mult_32 k7 ; y27 +---+ y17 +---+ subtr_node x27 ; o--<-|1/z|<-o--|1/z|<---+ limit_store y07 ; | +---+ +---+ | load_node x17 ; v | update_node x27 ; y05 +---+ +---+ +---+ | load_node y05 ; ->o-->| + |-->|*K7|--->| - |--->o---> update_node x17 ; | +---+ +---+ + +---+ y07 load_node y17 ; | - A update_node y27 ; | +---+ +---+ | load_node y07 ; +----|1/z|--o->|1/z|->-o update_node y17 ; +---+ x17 +---+ x27

• PWM - D/A-output

output y07 to PWM-DAC

test

output y01

ldd accu3, Y+y01+1 ;

subi accu3,$80 ; signed to unsigned $8000 ->$0000, $7FFF ->$FFFF sts OCR1AH,zero ; sts OCR1AL,accu3 ; fast PWM 8 bit y07 output to OCR1A ldi YL,Low(DatTbl1) ; back to 1st data table

output y08 to PWM-DAC

load_node y08 ; load y08

test outpt y02

load_node y02 ; load y02

subi datah,$80 ; signed to unsigned $8000 ->$0000, $7FFF ->$FFFF sts OCR1BH,zero ; sts OCR1BL,accu3 ; fast PWM 8 bit y08 output to OCR1B cbi PortD,0 ; just to measure INT-Time out SREG,savsrg ; restore status register reti ;

• coefficient table in flash memory

CoTblF: .dw 15709 ; k2 *32768 (0,4794008655888) .dw 28712 ; k4 *32768 (0,8762184935393) .dw 32001 ; k6 *32768 (0,9765975895082) .dw 32686 ; k8 *32768 (0,9974992559356) .dw 5301 ; k1 *32768 (0,1617584983677) .dw 24020 ; k3 *32768 (0,7330289323415) .dw 30977 ; k5 *32768 (0,9453497003291) .dw 32460 ; k7 *32768 (0,9905991566845)

• data in SRAM

.dseg .org $0100 DatTbl1: ; "allpass" 0 .dw $0000 ; x10 ; allpass 2 .dw $0000 ; x12 .dw $0000 ; x22 .dw $0000 ; y02 .dw $0000 ; y12 .dw $0000 ; y22 ; allpass 4 .dw $0000 ; x14 .dw $0000 ; x24 .dw $0000 ; y04 .dw $0000 ; y14 .dw $0000 ; y24 ; allpass 6 .dw $0000 ; x16 .dw $0000 ; x26 .dw $0000 ; y06 .dw $0000 ; y16 .dw $0000 ; y26 ; allpass 8 .dw $0000 ; x18 .dw $0000 ; x28 .dw $0000 ; y08 .dw $0000 ; y18 .dw $0000 ; y28 CoTblS: .dw $0000 ; k2 .dw $0000 ; k4 .dw $0000 ; k6 .dw $0000 ; k8 DatTbl2: .dw $0000 ; k1 .dw $0000 ; k3 .dw $0000 ; k5 .dw $0000 ; k7 ; allpass 1 .dw $0000 ; x11 .dw $0000 ; x21 .dw $0000 ; y01 .dw $0000 ; y11 .dw $0000 ; y21 ; allpass 3 .dw $0000 ; x13 .dw $0000 ; x23 .dw $0000 ; y03 .dw $0000 ; y13 .dw $0000 ; y23 ; allpass 5 .dw $0000 ; x15 .dw $0000 ; x25 .dw $0000 ; y05 .dw $0000 ; y15 .dw $0000 ; y25 ; allpass 7 .dw $0000 ; x17 .dw $0000 ; x27 .dw $0000 ; y07 .dw $0000 ; y17 .dw $0000 ; y27

</avrasm>