Digitalfilter mit ATmega

Wechseln zu: Navigation, Suche

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[Bearbeiten]

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 theoretisch eine unendliche Impulsantwort liefert. 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, neigt auf Grund der endlichen Genauigkeit der Koeffizienten eher zum Schwingen.

Software zur Koeffizientenberechnung[Bearbeiten]

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 kommerziellen Matlab 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[Bearbeiten]

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[Bearbeiten]

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) = H1(z)*H2(z)*H3(z). Wir müssen also H(z) in geeignete Faktoren zerlegen, auch dabei hilft uns die Mathematiksoftware.

Octave und Scilab[Bearbeiten]

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 ö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 z2 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[Bearbeiten]

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.

Es gibt noch weitere Verfahren, um "Wunschfrequenzgänge" zu berechnen:
prony in Matlab ( in Scilab nicht vorhanden) hier im Forum erwähnt - in ISBN 9780792395591 ab Seite 35 steht ein Matlab-m-file dazu, Update: es gibt inzwischen (seit2013) eine Signal processing toolbox für Scilab, die auch "prony" (und ähnlich auch "stmcb") enthält. und
FDLS (frequency domain least square) von Greg Berchin: Matlab-File und PDFs dazu, Diskussion auf fpgacentral.com Grundlagenartikel "Precise Filter Design" Auch als Kapitel 7 in ISBN 9781118316917 oder ? ISBN 9781118278383 (2.Auflage 2012) oder ISBN 9780470131572 (2007) ab Seite 59 zu finden

Von der Übertragungsfunktion zum Programm[Bearbeiten]

Unsere SOS-Teilfilter haben jetzt (je nach Vorzeichen der Exponenten) diese Übertragungsfunktion:


H(z)=\frac{b_{2}\cdot z^{-2}+b_{1}\cdot z^{-1}+b_{0}}{a_{2}\cdot z^{-2}+a_{1}*z^{-1}+1} = 
{\frac{b_{2}+b_{1}\cdot z+b_{0}\cdot z^{2}}{a_{2}+a_{1}\cdot z+z^{2}}}

Jetzt wird von z-Bereich wieder in den Zeitbereich rücktransformiert. Diese Funktion y(t)= f(x(t))und der dazugehörige Programmablauf wird üblicherweise in Form eines Blockdiagramms dargestellt:

IIR Block.png

Wir sehen drei Blocktypen: ein Verzögerung um einen Sampletakt Δt, Addition mehrere Summanden, und Multiplikationen mit Konstanten, die aus der Übertragungsfunktion abgelesen werden. Spezielle Signalprozessoren enthalten übrigens für die beiden letzteren einen sogenannten MAC multiply-accumulate Befehl. Im Atmega-Programm wird es in einem Assembler-Macro namens SMAC32 durchgeführt.

Literatur[Bearbeiten]

  • 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.
  • Bücher zu Scilab gibt es nur wenige, und diese gehen nicht auf das Thema Signalverarbeitung ein.
  • Für Januar 2010 angekündigt: ISBN 3527407243 oder 9783527407248 Signals and Systems Using SCILAB, Vorabinfo: kurzer Verlagstext PDF Seite 33 und bei Google-books.
  • ISBN 9781584502814 (2004) Digital signal processing fundamentals mit einem sehr kleinen Anhang "Scilab Tutorial", sehr oberflächlich
  • ISBN 9780817640095 (1999) Engineering and scientific computing with Scilab mit einem Kapitel zu "Signal Processing", für den Einstieg in Scilab sehr zu empfehlen. Auch als unveränderter Nachdruck von 2007 noch neu erhältlich.

Scilab-Berechnung mit eqiir[Bearbeiten]

// 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

Messung des 1750 Hz-Sperrfilters[Bearbeiten]

Gemessen mit "Visual Analyzer":

1750HzNotch.png

Frequenzgang-Planung mit Yulewalk[Bearbeiten]

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.

Massgeschneiderter Frequenzgang mit yulewalk

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:

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

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[Bearbeiten]

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

 
;***************************************************************************
;* 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
;*************************************************************************