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 1200def.inc which defines the chip-specific constants used in this code.
The algorithm follows a very basic flow:
- Reset
- Set up ports, timer, and initial state variables.
- Enter main program loop and loop forever.
- Interrupt every 1ms. To get the timing correct, the interrupt occurs every .25ms, but the code waits 4 such cycles before executing.
- Checks the input pins, and sums up the "weighted current" coming into the chip.
- Update GK.
- Calculate intermediate variable: constMemb.
- Update membrane voltage.
- Update threshold voltage.
- Compare current membrane voltage to current threshold. Set PORTB low if voltage exceeds threshold.
- 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.
Subroutines
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.
Code
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\1200def.inc"
.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
;constGK:
; 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
.cseg
.org $0000
rjmp RESET ;RESeT
reti ;external interrupt
rjmp TIMER ;timer 0
reti ;analog comprator
;**************RESET
RESET:
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
rjmp COUNT_CURRENT
timer_nop:
;Update GK###########################################
ldi xi,constGKderVali
ldi xf,constGKderValf
tst Spike
brne gk_next
clr xi
clr xf
gk_next:
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
;constMemb#############################################
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
;Vmemb####################################################3
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
vmemb_next:
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
;threshold###############################################
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
;Spike?###########################################################
cp vmembi, thresholdi
brlo ELSE
ser Spike
rjmp SPIKEY_TEST
ELSE: clr SPike
SPIKEY_TEST:
tst Spike
breq OFF
clr temp
out PORTB, temp
rjmp timer_end
OFF: ser temp
out PORTB, temp
timer_end:
ldi tcount,time ;resets the timer counter
NOT_YET:
ldi temp, ovflCts ;preload timer=250 counts till overflow
out TCNT0, temp ;250x1024x.25 microSec = 64mSec.
out SREG, save ;restore status reg
reti
;**************count current
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)
end_e:
mov XbaseEi, save1
mov XbaseEf, save2
mov base, save3
mov iln2i, save4
mov iln2f, save5
ret
;***************** 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
mpy16s_under:
ldi mcnt16s, 0xFF
cp m16s3, mcnt16s
breq end_mpy16s
ldi m16s2, 0x80 ;underflow --> -128
clr m16s1
end_mpy16s:
mov mcnt16s, save1
ret
;***************************************************************************
;*
;* "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
;save
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
end_mpy16u:
mov mcnt16u, save1
ret
;*********************************************
;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)
div_end:
mov dcnt16s,save1 ;restore and return
mov dv16sL, save2
mov dv16sH, save3
mov dd16s1, save4
mov d16s, save5
ret