Source Code

The source code for this program was written in Atmel's AVR Assembler 1.30. It is a modified translation of MacGregor's ptnrn10 program (see Neuronal Modeling). The program is a total of 308 words and runs on the AT90S1200 MCU. It requires the use of which defines the chip-specific constants used in this code.

The algorithm follows a very basic flow:

  1. Reset
  2. Set up ports, timer, and initial state variables.
  3. Enter main program loop and loop forever.
  4. Interrupt every 1ms. To get the timing correct, the interrupt occurs every .25ms, but the code waits 4 such cycles before executing.
  5. Checks the input pins, and sums up the "weighted current" coming into the chip.
  6. Update GK.
  7. Calculate intermediate variable: constMemb.
  8. Update membrane voltage.
  9. Update threshold voltage.
  10. Compare current membrane voltage to current threshold. Set PORTB low if voltage exceeds threshold.
  11. Return to main loop.

The code can be physically divided into 4 parts. The first part is the interrupt vector table. Only the reset and timer interrupts are enabled. The second part of the program is the initialization code which sets up the ports, timers, and other variables. Immediately following the reset code is the 1-line main loop. This loop just jumps back to itself until power stops flowing to the chip.

Timer ISR

The third part is the timer ISR (interrupt service routine). The timer has two purposes: to sample the "input current" and for time-stepping the numerical integration in real-time. The timer ISR also controls the output signaling. All the mathematical calculations are done in this ISR. Although this is bad coding technique, it does not matter because this the only function of the code running on the chip.


The fourth part of the program are the 5 subroutines: count_current, signed 16-bit multipy, unsigned 16-bit multiply, signed divide, and exponential. The count_current routine checks the pins on port D and sums them together. If anybody knows a better way of doing, please let me know. The other four subroutines are explained in the Numerical Methods page.

To prevent subroutine data from colliding with program data, I wrote all the subroutines to follow two a few basic rules: operand register values are not preserved, any other registers that might get overwritten are saved in specific save registers, and the parameters and results reside in specific registers.

A full listing of the register file assignemnts can be found on the Registers page.


This is a full listing of point.asm that resides in the chips for the fish swim network.

;point.asm - A Derivation of Point Neuron 10 from MacGregor
;This program computes action potential output from a neural simulation
;Programmed by: Michael Katz, Cornell University, EE '01 for ELE E 491

;uses the AT90S1200 MCU
.include "c:\windows\desktop\projects\"

.device AT90S1200 ; uses the Ateml AVR 1200

;*********operation variables
.def op1i = r12		;operand1 - integer part
.def op1f = r13		;operand1 - fractional part
.def op2i = r14		;operand2 - integer part
.def op2f = r15		;operand2 - fractional part
.def resi = r16		;result1 - integer
.def resf = r17		;result1 - fraction

;**********maintenance variables
.def i = r10		;counter
.def inputz = r11	;pin input
.def save = r0		;saves status register
.def tcount = r18	;timer counter
.def temp = r20		;temporary register
.def temp2 = r21
.def xi = r22
.def xf = r23
;*********state variables
.def current = r24	;the current running through the circuit
.def thresholdi = r25	;threshold
.def thresholdf = r26
.def GKi = r28		;potassium conducatance - integer
.def GKf = r29
.def Vmembi = r30	;membrane voltage
.def Vmembf = r31

.def spike = r27	;spike hi/low

;+++++++++++Timer Stuff
.equ PreScal =2 	;constant to set timer prescale to 8
.equ ovflCts =6		;overflow counter 256-250 = 6
.equ time    =4		;count time in .25 msec increments= 1 ms

;++++++++++Bio Stuff
.equ tMeb = 5		;Mebrane = 5ms

.equ ampThreshold = 0x80 ; =.5 | threshold (0-1)
.equ tThreshold = 25	;ms
.equ initThreshold = 0	;mV above resting (dft 10mV)

.equ ampGK = 1		;potassium conducatance
.equ tGK = 90		;ms
.equ EK = -10		;eq'm potential from resting (mV)

;++++++++++++++++integration contstants
;constTH :
;	tThreshold =25, step = 10 ---> e^(-step/tThreshold) = .670320046036
;	  dividing by 1/255, constTH = 170.931611739, or 0xAB (rounding up)
.equ	constTH = 0xF5
.equ	OMconstTH = 0x0B

;	tGK = 5, step = 10 ---> e^(-step/tGK) = .135335283237
;	  dividing by 1/255, constGK = 34.5104972254, or 0x23 (rounding up)
.equ	constGK = 0xFC
.equ	constGKderVali = 0x00	;(1 - constGK)*ampGK = 0x00.03
.equ	constGKderValf = 0x03

;step/tmemb = 1/5 = .2
.equ sdtmb = 0x33

.equ step = 1

; ******** Initialiazation
.org $0000
	rjmp RESET 	;RESeT
	reti		;external interrupt
	rjmp TIMER	;timer 0
	reti		;analog comprator

        ser     temp       	;set PORTB to output
        out     DDRB,Temp      

	clr	temp		;set PORTD to input
	out	DDRD, Temp 
	ldi	temp,0x3F
	out	PORTD, Temp

	;timer stuff+++++++++++++++++++++++++++++++
        ldi     temp,exp2(TOIE0);enable timer interrupt
        out     TIMSK, Temp
        ldi     temp, PreScal   ;prescale timer by 1024
        out     TCCR0, temp
	ldi	temp, ovflCts 
        out     TCNT0, temp   ; 125x64x.25 microSec = 2mSec.
        ldi     tcount, time     ;set up sampling interval

	;state initilizations
	clr	Vmembi
	clr	Vmembf
	ldi	thresholdi, initThreshold
	clr	thresholdf
	ldi	Spike, 0x00
	clr	GKi
	clr	GKf

	sei			;enable all interrupts

MAIN:   rjmp Main

;************timer interrupt
TIMER:  in      save, SREG      ;save the status reg0
	dec     tcount          ;keep track of # of intr
	breq	next_step	;=.96 sec? then:
        rjmp    NOT_YET         
next_step:  			;gets the current
	in	inputz, PIND	;saves the PIND inputs

	;Update GK###########################################
	ldi 	xi,constGKderVali
	ldi	xf,constGKderValf
	tst	Spike
	brne	gk_next
	clr	xi
	clr	xf
	mov	op1i,GKi
	mov	op1f,GKf
	ldi	temp,constGK
	mov	op2f,temp
	clr	op2i
	rcall	mpy16u		;GK*constGK

	add	xf, resf		;Add fractional bytes
	adc	xi, resi		;Add integer bytes with carry	

	mov	GKi,xi			;GK=constGK*GK + Spike*ampGK*(1-constGK)
	mov	GKf,xf
	inc 	xi		;GK+1
	mov	op1i,xi
	mov	op1f,xf
	ldi	temp,sdtmb
	clr	op2i
	mov	op2f,temp
	rcall	mpy16s		;(1+GK)*(step/tMemb)

	mov 	op1i, resi
	mov 	op1f, resf
	rcall 	exp		; resi.resf = e^(-(1+GK)*step/tmemb)
	mov	xi, resi
	mov 	xf, resf	;x = constMemb
	ldi temp, EK
	mov op1i, temp
	clr op1f
	mov op2i, GKi
	mov op2f, Gkf
	rcall	mpy16s		;GK*EK
	mov temp, xi
	mov temp2, xf
	add resi, current	;current + GK*EK
	neg xi
	neg xf			;1-constMemb
	mov op1i, resi
	mov op1f, resf
	mov op2i, xi
	mov op2f, xf
	rcall mpy16s		;(current + GK*EK)*(1-constMemb)
	mov op2i,GKi
	inc op2i
	mov op2f,GKf
	mov op1i,resi
	mov op1f,resf
	;(current + GK*EK)*(1-constMemb)/(1+GK)
	;modified divide routine for signed fixed point
	rcall div16s		
	mov op1i,temp	;puts constMemb in the multicand
	mov op1f,temp2
	mov temp,resi	;saves previous division
	mov temp2,resf
	mov op2i,Vmembi
	mov op2f,Vmembf
	rcall mpy16s		;Vmemb*constMemb
	mov Vmembi, resi
	mov Vmembf, resf
	add Vmembf, temp2	;Vmemb = Vmemb*constMemb + (current + GK*EK)*(1-constMemb)/(1+GK)
	adc Vmembi, temp

	ldi temp, ampThreshold
	mov op1f, temp
	clr op1i
	mov op2i, Vmembi
	mov op2f, Vmembf
	rcall mpy16s		;Vmemb*ampThreshold
	ldi temp, initThreshold 
	add resi, temp		;Vmemb*ampThreshold + initThreshold

	mov op2i, resi
	mov op2f, resf
	ldi temp, OMconstTH
	mov op1f, temp		;1-constTH
	clr op1i
	rcall mpy16s		;(Vmemb*ampThreshold+initThreshold)*(1-constTH)			
	mov temp, resi
	mov temp2,resf

	ldi xf, constTH	
	mov op1f, xf
	clr op1i
	mov op2i, thresholdi
	mov op2f, thresholdf
	rcall mpy16s		 ;threshold*constTH
	mov thresholdf, resf
	mov thresholdi, resi
	add thresholdf, temp2	;threshold = threshold*constTH + (Vmemb*ampThreshold+initThreshold)*(1-constTH)
	adc thresholdi, temp 

	cp vmembi, thresholdi
	brlo ELSE
	ser Spike
ELSE:	clr SPike

	tst Spike
	breq OFF
	clr temp
	out PORTB, temp
	rjmp timer_end

OFF:	ser temp 
	out PORTB, temp

	ldi	tcount,time	 ;resets the timer counter
	ldi     temp, ovflCts  	;preload timer=250 counts till overflow
        out     TCNT0, temp	;250x1024x.25 microSec = 64mSec.
        out     SREG, save      ;restore status reg

;**************count current
	clr	current
	sbrs	inputz,0
	inc	current
	sbrs	inputz,1
	inc	current
	sbrs	inputz,2
	inc	current
	sbrs	inputz,3
	inc	current
	sbrs	inputz,4
	inc	current
	sbrs	inputz,5
	inc	current
	sbrs	inputz, 6	;PIN6 has the negative 20amp weight for hyperpolarization
	subi	current, 20
	rjmp	timer_nop

;* "exp" - 8-bit negative exponential
;* This subroutine takes the exponential of negative op1i.opt1f
;* The result is placed in resi.resf
;* The program requires the use of mpy16um  

;***** Subroutine Register Variables
.def op1i = r12
.def op1f = r13
.def op2i = r14
.def op2f = r15
.def XbaseEi = r22
.def XbaseEf = r23
.def resi = r16
.def resf = r17
.def iln2i = r25
.def iln2f = r26
.def base = r24
.def save1 = r1
.def save2 = r2
.def save3 = r3
.def save4 = r4
.def save5 = r5

;ln2 = .69314718056 --> (8bit) 176.75231043
.equ ln2 = 0xB1
;1/ln2 = 1.44269504089 = 1 + 112.87235427 = 0x01.71

;the linear interpolation of 2^x, 0 slope corresponding to x = .25 = 0x40, y = 1.1892 
.equ p2 = 0xD3 ; .8284 --> slope corresponding to x = .5  = 0x80, y = 1.4142
.equ p3 = 0xE8 ; .9091 --> slope corresponding to x = .75 = 0xC0, y = 1.6818
.equ p4 = 0x00; 1 --> slope corresponding to x = 1, y =2 

exp:	mov save1, XbaseEi ; saves the inputs
	mov save2, XbaseEf
	mov save3, base
	mov save4, iln2i
	mov save5, iln2f

	mov XbaseEi, op1i

	clr resi	;sets the output to 0
	clr resf

;	neg XbaseEi		;negates the integer part

	;check to make sure there isn't an underflow
	cpi XbaseEi, 0x08
	brge end_e

	;multiply x by 1/ln2
	ldi iln2i, 0x01
	ldi iln2f, 0x71
	mov op1i, XbaseEi
	mov op2i, iln2i
	mov op2f, iln2f
	rcall mpy16u

;convert the integer part
	ldi base, 0x80		;divides by 2 for each integer value of x
	tst resi
	breq frac_part
loop:	lsr base		;div by 2
	dec resi
	tst resi		;=0?
	brne loop

frac_part: ;convert the fractional part
	ldi iln2f, 0x00
	sub iln2f, resf		;1-decimal part of x

	clr XbaseEi
;check the range
	ldi XbaseEf, p1
	cpi iln2f, 0x40
	brlo zoiks		;if x<.25, slope = p1
	ldi XbaseEf, p2
	cpi iln2f, 0x80
	brlo zoiks		;if x<.5, slope = p2
	ldi XbaseEf, p3
	cpi iln2f, 0xC0
	brlo zoiks		;if x<.75, slope = p3
	ldi XbaseEi, 0x01
	ldi XbaseEf, p4		;if x==.75, slope = 1.p4

zoiks:	clr op2i
	mov op2f, iln2f
	mov op1i,XbaseEi
	mov op1f,XbaseEf
	rcall mpy16u		;x*slope

	inc resi		;x*slope+1, to fit linear interp of 2^x, -1 < x < 0
	mov op2i, resi
	mov op2f, resf
	clr iln2i
	mov op1i, iln2i
	mov op1f, base
	rcall mpy16u		;2^integer*(liner interp of decimal)

	mov XbaseEi, save1
	mov XbaseEf, save2
	mov base, save3
	mov iln2i, save4
	mov iln2f, save5

;*****************  resi.resf = e^x, -8 < x < 0

;* "mpy16s" - 16x16 Bit Signed Multiplication
;* This subroutine multiplies signed the two 16-bit register variables 
;* mp16sH:mp16sL and mc16sH:mc16sL.
;* The result is placed in m16s3:m16s2:m16s1:m16s0.
;* The routine is an implementation of Booth's algorithm. If all 32 bits
;* in the result are needed, avoid calling the routine with
;* -32768 ($8000) as multiplicand
;* Number of words	:16 + return
;* Number of cycles	:210/226 (Min/Max) + return
;* Low registers used	:None
;* High registers used  :7 (mp16sL,mp16sH,mc16sL/m16s0,mc16sH/m16s1,
;*			    m16s2,m16s3,mcnt16s)	

;***** Subroutine Register Variables

.def	mc16sL	=r13		;multiplicand low byte
.def	mc16sH	=r12		;multiplicand high byte
.def	mp16sL	=r15		;multiplier low byte
.def	mp16sH	=r17		;multiplier high byte
.def	m16s0	=r15		;result byte 0 (LSB)
.def	m16s1	=r17		;result byte 1
.def	m16s2	=r16		;result byte 2
.def	m16s3	=r9		;result byte 3 (MSB)
.def	mcnt16s	=r31		;loop counter
.def 	save1	=r1

;***** Code
mpy16s:	mov	save1, mcnt16s
	mov	mp16sH, op2i

	clr	m16s3		;clear result byte 3
	sub	m16s2,m16s2	;clear result byte 2 and carry
	ldi	mcnt16s,16	;init loop counter
m16s_1:	brcc	m16s_2		;if carry (previous bit) set
	add	m16s2,mc16sL	;    add multiplicand Low to result byte 2
	adc	m16s3,mc16sH	;    add multiplicand High to result byte 3
m16s_2:	sbrc	mp16sL,0	;if current bit set
	sub	m16s2,mc16sL	;    sub multiplicand Low from result byte 2
	sbrc	mp16sL,0	;if current bit set
	sbc	m16s3,mc16sH	;    sub multiplicand High from result byte 3
	asr	m16s3		;shift right result and multiplier
	ror	m16s2
	ror	m16s1
	ror	m16s0
	dec	mcnt16s		;decrement counter
	brne	m16s_1		;if not done, loop more	

	tst 	m16s3
	breq	end_mpy16s
	brlt	mpy16s_under	;negative number --> underflow
	ser 	m16s1
	ldi	m16s2, 0x7F	;overflow --> 127.99999
	rjmp 	end_mpy16s

	ldi	mcnt16s, 0xFF
	cp	m16s3, mcnt16s
	breq	end_mpy16s
	ldi 	m16s2, 0x80	;underflow --> -128
	clr	m16s1

	mov	mcnt16s, save1

;* "mpy16u" - 16x16 Bit Unsigned Multiplication
;* This subroutine multiplies the two 16-bit register variables 
;* mp16uH:mp16uL and mc16uH:mc16uL.
;* The result is placed in m16u3:m16u2:m16u1:m16u0.
;* Number of words	:14 + return
;* Number of cycles	:153 + return
;* Low registers used  :7 (mp16uL,mp16uH,mc16uL/m16u0,mc16uH/m16u1,m16u2,
;*                          m16u3,mcnt16u)	
;* High registers used	:None

;***** Subroutine Register Variables

.def	mc16uL	=r13		;multiplicand low byte
.def	mc16uH	=r12		;multiplicand high byte
.def	mp16uL	=r15		;multiplier low byte
.def	mp16uH	=r17		;multiplier high byte
.def	m16u0	=r15		;result byte 0 (LSB)
.def	m16u1	=r17		;result byte 1
.def	m16u2	=r16		;result byte 2
.def	m16u3	=r9		;result byte 3 (MSB)
.def	mcnt16u	=r31		;loop counter
.def	save1 = r1

;***** Code

mpy16u:	mov	save1, mcnt16u
	mov	mp16uH, op2i

	clr	m16u3		;clear 2 highest bytes of result
	clr	m16u2
	ldi 	mcnt16u, 16	;init loop counter
	lsr	mp16uH
	ror	mp16uL

m16u_1:	brcc	noad8		;if bit 0 of multiplier set
	add	m16u2,mc16uL	;add multiplicand Low to byte 2 of res
	adc	m16u3,mc16uH	;add multiplicand high to byte 3 of res
noad8:	ror	m16u3		;shift right result byte 3
	ror	m16u2		;rotate right result byte 2
	ror	m16u1		;rotate result byte 1 and multiplier High
	ror	m16u0		;rotate result byte 0 and multiplier Low
	dec	mcnt16u		;decrement loop counter
	brne	m16u_1		;if not done, loop more

	;overflow handling
	tst	m16u3
	breq	end_mpy16u
	ser	m16u2		;if overflow, result = 0xFF.FF
	ser	m16u1

	mov	mcnt16u, save1

;NewDiv.asm --> div16s -> based on a 32-bit by 16 divide scheme
;highly specialized for point.asm
;48 words, 405 cycles, 101.25us

;***** Subroutine Register Variables

.def	d16s	=r30		;sign register
.def	drem16sL=r9		;remainder low byte		
.def	drem16sH=r11		;remainder high byte
.def	drem16s3=r10
.def	dres16s2=r16		;result middle byte
.def	dres16s1=r18		;result highest byte
.def	dres16s3=r17		;result lowest byte
.def	dd16s2	=r16		;dividend low byte
.def	dd16s1	=r18		;dividend high byte
.def	dd16s3	=r17		;dividend zero byte
.def	dv16sL	=r19		;divisor low byte
.def	dv16sH	=r20		;divisor high byte
.def	dcnt16s	=r31		;loop counter
.def 	save1	=r1
.def	save2	=r2
.def	save3	=r3
.def	save4 	=r4
.def	save5	=r5

;***** Code

div16s:	mov	save1, dcnt16s	
	mov	save2, dv16sL
	mov	save3, dv16sH
	mov	save4, dd16s1
	mov	save5, d16s
	mov	dd16s2, op1f
	mov	dd16s1, op1i
	mov	dv16sL, op2f
	mov	dv16sH, op2i
	clr	d16s		;clear sign register
	tst	dd16s1		;if dividend negative, set sign register
	brge	div_0
	ser	d16s
	com	dd16s1		;    change sign of dividend
	com	dd16s2		
	subi	dd16s2,low(-1)
	sbci	dd16s2,high(-1)

div_0:	clr	dd16s3		;make 24-bit #

	clr	drem16sL	;sets remainder to 0
	clr	drem16sH

	ldi	dcnt16s, 0x18	;sets counter at 24

div_1:	clc

	rol	dres16s3	;shift dividend left
	rol	dres16s2
	rol	dres16s1
	rol	drem16sL	;shift remainder left + carry in from dividend
	rol	drem16sH

	sub	drem16sL, dv16sL ;remainder = remainder - divisor
	sbc	drem16sH, dv16sH

	brcc	div_2		;negative?
	add	drem16sL, dv16sL ; restore remaidner
	adc	drem16sH, dv16sH
	rjmp	div_3

div_2:	ori	dres16s3, 0x01	;positive --> Q0 = 1

div_3:	dec 	dcnt16s		;count--
	brne	div_1
	tst	d16s
	brge	div_end
	com	dres16s2	;        change sign of result
	com	dres16s3
	subi	dres16s3,low(-1)
	sbci	dres16s2,high(-1)

	mov	dcnt16s,save1	;restore and return
	mov	dv16sL,	save2
	mov	dv16sH, save3
	mov	dd16s1, save4
	mov	d16s,	save5