1
$\begingroup$

I'm trying to simulate a compound action potential for calibrating research instruments. The image below is a small overview of the design. enter image description here

I've already completed the data acquisition from a live animal, and processed the data in MATLAB to make a nice, noise-less signal, with 789 values in 12-bit format. I then cloned the repository to the Pi using Git. Below is the Python script I've written on the RPi. You can skip to def main in the script to see functionality.

#!/usr/bin/python

import spidev
from time import sleep
import RPi.GPIO as GPIO
import csv
import sys
import math

DEBUG = False
spi_max_speed = 20 * 1000000
V_Ref = 5000
Resolution = 2**12
CE = 0

spi = spidev.SpiDev()
spi.open(0,CE)
spi.max_speed_hz = spi_max_speed

LDAQ = 22
GPIO.setmode(GPIO.BOARD)
GPIO.setup(LDAQ, GPIO.OUT)
GPIO.output(LDAQ,GPIO.LOW)

def setOutput(val):
    lowByte = val & 0b11111111 #Make bytes using MCP4921 data sheet info
    highByte = ((val >> 8) & 0xff) | 0b0 << 7 | 0b0 << 6 | 0b1 << 5 | 0b1 << 4
    if DEBUG :
        print("Highbyte = {0:8b}".format(highByte))
        print("Lowbyte =  {0:8b}".format(lowByte))
    spi.xfer2([highByte, lowByte])

def main():
    with open('signal12bit.csv') as signal:
        signal_length = float(raw_input("Please input signal length in ms: "))
        delay = float(raw_input("Please input delay after signal in ms: "))
        amplitude = float(raw_input("Please input signal amplitude in mV: "))
        print "Starting Simulant with signal length %.1f ms, delay %.1f ms and amplitude %.1f mV." % (signal_length, delay, amplitude)
        if not DEBUG : print "Press ctrl+c to close."
        sleep (1) #Wait a sec before starting
        read = csv.reader(signal, delimiter=' ', quotechar='|')
        try:
            while(True):
                signal.seek(0)
                for row in read: #Loop csv file rows
                    if DEBUG : print ', '.join(row)
                    setOutput(int(row)/int((V_Ref/amplitude))) #Adjust amplitude, not super necessary to do in software
                    sleep (signal_length/(data_points*1000) #Divide by 1000 to make into ms, divide by length of data
                sleep (delay/1000)
        except (KeyboardInterrupt, Exception) as e:
            print(e)
            print "Closing SPI channel"
            setOutput(0)
            GPIO.cleanup()
            spi.close()

if __name__ == '__main__':
    main()

This script almost works as intended. Connecting the output pin of the MCP4921 to an oscilloscope shows that it reproduces the signal very well, and it outputs the subsequent delay correctly.

Unfortunately, the data points are seperated much further than I need them to be. The shortest time I can cram the signal into is about 79 ms. This is due to dividing by 789000 in the sleep function, which I know is too much to ask from Python and from the Pi, because reading the csv file takes time. However, if I try making an array manually, and putting those values out instead of reading the csv file, I can achieve a frequency over 6 kHz with no loss.

My question is this

How can I get this signal to appear at a frequency of 250 Hz, and decrease it reliably from the user's input? I've thought about manually writing the 789 values into an array in the script, and then changing the SPI speed to whatever value fits with 250 Hz. This would eliminate the slow csv reader function, but then you can't reduce the frequency from user input. In any case, eliminating the need for csv.read would help a lot. Thanks!

$\endgroup$
  • $\begingroup$ I'm confused. You said, "Connecting the output pin... to an oscilloscope shows that it reproduces the signal very well, and it outputs the subsequent delay correctly."Immediately after that, you say, "Unfortunately, the data points are separated much further than I need them to be." So, which is it? Is the signal reproduced very well with correct delays, or is it not? Can you post pictures of your oscilloscope output and highlight what's wrong with the signal? $\endgroup$ – Chuck Nov 1 '16 at 13:04
  • $\begingroup$ There are two delays, as you can see in the code. One delay is between each of the 789 data points, and one delay is after the entire signal. The delay after the signal is the one that works correctly. I usually set it to 50 ms, and the oscilloscope reading shows around 50 ms as well. The delay between data points, the very short one, is the problem. I'd post a picture of the oscilloscope reading if I could, but I can't until tomorrow morning. $\endgroup$ – Tachyon Nov 1 '16 at 13:08
0
$\begingroup$

I posted a comment to the question asking you to provide more information about what your problem is exactly.

I can kind of side-step that issue by acknowledging that it seems you have some kind of question about timing with Python on the Raspberry Pi, to which I would redirect you to this answer on StackOverflow. In summary, for Linux operating systems (like Raspian), "non-realtime Linux kernels have minimum sleep interval much closer to 1ms then 10ms but it varies in a non-deterministic manner."

Now, for a rate of 250Hz, you're asking for a 4ms delay between samples. A side/coding comment of mine at this point would be - instead of asking the user for a signal duration, why not ask for the sample rate? Then your sampling delay is the sample period, which is 1/sampleRate.

In any event, the sample period timing you're asking for doesn't appear to be outside the limits of what's achievable with Python's timing. This would seem to me to imply, then, that your hunch may be correct - things you are doing outside of the sleep itself are impacting the timing of your program.

Now the question is - why are you doing those things? Let's take a look at the essence of what you're doing:

while(True):
    signal.seek(0)
    for row in read: #Loop csv file rows
        if DEBUG : print ', '.join(row)
        setOutput(int(row)/int((V_Ref/amplitude))) #Adjust amplitude, not super necessary to do in software
        sleep (signal_length/(data_points*1000) #Divide by 1000 to make into ms, divide by length of data
    sleep (delay/1000)

Your signal never changes - it's hard-coded into the source code. Despite this, you still have to manually input the signal parameters by hand every time the code is run (not sure why you chose to do this). But all that's beside the point - what I mean to say is, once you get to the loop, you never swap out a different signal. As described, your signal isn't very big either.

My advice to you would be to move all of your data processing outside of your loop so the only things you have to do are transmit the signal and sleep. No math. For example, consider the following changes:

signal.seek(0)
signalValues = [int(row)/int((V_ref/amplitude)) for row in read]
samplePeriod = signal_length/(data_points*1000)
signalDelay = delay/1000
while(True):
    for currentVal in signalValues:
        if DEBUG : print ', '.join(currentVal)
        setOutput(currentVal)
        sleep (samplePeriod)
    sleep (signalDelay)

Hopefully you see here that all of the floating point division and CSV reading has been moved outside of the loop. None of your parameters change during the loop so there's no point in waiting to modify them on-the-fly.

Load the entire CSV (all 789 values) into a list that's held in memory. Scale them however they need to be scaled (warning on that later), then just iterate your way through the list.

About scaling your values: You ask for the inputs as float, but then skew your result by casting them to int mid-division. Specifically, you have the line:

setOutput(int(row)/int((V_Ref/amplitude)))

Maybe you are thinking that an int divided by and int forces integer division - this is not correct for Python 3.0. If you are using Python 3.0 or later, then you're very likely getting a floating point number and then doing binary bit-wise operations on it when you attempt to send it. This could result in Bad Things happening.

More to the point, if your desired output is an int, then I would suggest you avoid prematurely truncating your numbers by casting the entire result to int as opposed to casting the individual operands. Specifically, consider this alternative:

setOutput(int(row/(V_Ref/amplitude)))

Here the output is guaranteed to be an int, because you cast it that way, and you get better numeric accuracy in which integer results from your math because you waited until the end of the operation to truncate.

If you want integer division, you can get that with the // integer division operator - int(row)//int(V_Ref/amplitude) - but as with my comment about waiting to the end to cast as int, here the integer division operator guarantees an integer output; it accepts floating point inputs. So you could also have row // (V_Ref/amplitude) and should give the same result as int(row/(V_Ref/amplitude)).

So, in summary:

  1. For speed, consider converting your CSV to a list.
  2. For speed, consider pulling your static math operations (those whose operands do not change during your loops) outside of the loop.
  3. For accuracy, consider waiting to cast values to integer until after all math operations are complete.
$\endgroup$
  • $\begingroup$ Thanks for the reply! This is my first Python script, so I know it's very bad. However, I seem to have confused some people. I need the entire signal frequency to be 250 (or less, if the user wants it). With 789 values, this means I need a sample each 45 µs, or a sample rate of about 22000. I've found a page with the exact quote you gave, "non-realtime Linux kernels have minimum sleep interval much closer to 1ms then 10ms but it varies in a non-deterministic manner." To me, this implies that using sleep for the delay between the samples is impossible. $\endgroup$ – Tachyon Nov 1 '16 at 14:08
  • $\begingroup$ @Tachyon - I think I understand, you want to expand or compress the signal as required to broadcast at an arbitrary rate? In that event, like you said, the sleep delay is not going to be capable of what you want to do. In this instance, I would consider revision your tool chain - Consider using Matlab (or Python) to re-package the data back into a WAV file and then play the wav file. This means you'd be using the audio output instead of a DAC, but the handy thing there is that an audio processor uses a DAC to build the wave form. $\endgroup$ – Chuck Nov 1 '16 at 14:26
  • $\begingroup$ Also, you're free to specify the sampling rate of the WAV file to be whatever you want, which then should make it easy to vary the "width" of the signal by specifying a different sample rate. For the inter-signal break, you could either write a temporary WAV file with empty (zero) samples for the break or you could try to sleep between the signals. $\endgroup$ – Chuck Nov 1 '16 at 14:30
  • $\begingroup$ Whatever the case, I would urge you to do some bounds-checking on user inputs. If the user wants an inter-signal delay of 4ms, then you don't get any time to actually broadcast the signal. 1/(0 + 4)ms = 250Hz. You mention that the output of this is going to calibrate research instruments - you could use a unity op-amp to buffer the audio output to isolate the RPi from the instrument. $\endgroup$ – Chuck Nov 1 '16 at 14:32
  • $\begingroup$ Using an audio DAC would work, of course, but the whole point of the project is to calibrate a sensor to this analog voltage signal generated by the DAC. Yes, I thought about bounds-checking, but figured I'd add it once the system actually does what it's supposed to. It seems at this point, I'll have to remove the delay between samples, and instead adjust the SPI exchange rate. In that case I'll lose control over the signal length, but will probably be able to achieve something close to 250 Hz. $\endgroup$ – Tachyon Nov 1 '16 at 14:37
0
$\begingroup$

Figured it out earlier today, so I thought I'd post an answer here, in case someone comes upon a similar problem in the future.

The problem with the internal delay between data points cannot be solved with sleep(), for several reasons. What I ended up doing was the following

  • Move all math and function calling out of the critical loop
  • Do a linear regression analysis on the time it takes to transfer the values with no delay
  • Increase the number of datapoints in the CSV file to "plenty" (9600) in MATLAB
  • Calculate the number of points needed to meet the user's wanted signal length
  • Take evenly seperated entries from the now bigger CSV file to fit that number of points as closely as possible.
  • Calculate these values and then calculate the SPI bytes explicitly
  • Save the two byte lists, and output them directly in the critical loop

The new code, with a bit of input checking, is below

#!/usr/bin/python

import spidev
from time import sleep
import RPi.GPIO as GPIO
import sys
import csv
import ast

spi_max_speed = 16 * 1000000 # 16 MHz
V_Ref = 5000 # 5V in mV
Resolution = 2**12 # 12 bits for the MCP 4921
CE = 0 # CE0 or CE1, select SPI device on bus
total_data_points = 9600 #CSV file length

spi = spidev.SpiDev()
spi.open(0,CE)
spi.max_speed_hz = spi_max_speed

LDAQ=22
GPIO.setmode(GPIO.BOARD)
GPIO.setup(LDAQ, GPIO.OUT)
GPIO.output(LDAQ,GPIO.LOW)

def main():

    #User inputs and checking for digits
    signalLengthU = raw_input("Input signal length in ms, minimum 4: ")
    if signalLengthU.isdigit():
        signalLength = signalLengthU
    else:
        signalLength = 4

    delayU = raw_input("Input delay after signal in ms: ")
    if delayU.isdigit():
        delay = delayU
    else:
        delay = 0

    amplitudeU = raw_input("Input signal amplitude in mV, between 1 and 5000: ")
    if amplitudeU.isdigit():
        amplitude = amplitudeU
    else:
        amplitude = 5000

    #Calculate data points, delay, and amplitude
    data_points = int((1000*float(signalLength)-24.6418)/12.3291)
    signalDelay = float(delay)/1000
    setAmplitude = V_Ref/float(amplitude)

    #Load and save CSV file
    datain = open('signal12bit.csv')
    read = csv.reader(datain, delimiter=' ', quotechar='|')
    signal = []
    for row in read:
        signal.append(ast.literal_eval(row[0]))

    #Downsampling to achieve desired signal length
    downsampling = int(round(total_data_points/data_points))
    signalSpeed = signal[0::downsampling]
    listlen = len(signalSpeed)

    #Construction of SPI bytes, to avoid calling functions in critical loop
    lowByte = []
    highByte = []
    for i in signalSpeed:
        lowByte.append(int(i/setAmplitude) & 0b11111111)
        highByte.append(((int(i/setAmplitude) >> 8) & 0xff) | 0b0 << 7 | 0b0 << 6 | 0b1 << 5 | 0b1 << 4)

    print "Starting Simulant with signal length %s ms, delay %s ms and amplitude %s mV." % (signalLength, delay, amplitude)
    print "Press ctrl+c to stop."
    sleep (1)

    try:
        while(True): #Main loop
            for i in range(listlen):
                spi.xfer2([highByte[i],lowByte[i]]) #Critical loop, no delay!
            sleep (signalDelay)
    except (KeyboardInterrupt, Exception) as e:
        print e
        print "Closing SPI channel"
        lowByte = 0 & 0b11111111
        highByte = ((0 >> 8) & 0xff) | 0b0 << 7 | 0b0 << 6 | 0b1 << 5 | 0b1 << 4
        spi.xfer2([highByte, lowByte])
        GPIO.cleanup()
        spi.close()

if __name__ == '__main__':
    main()

The result is exactly what I wanted. Below is seen an example from the oscilloscope with a signal length of 5 ms; 200 Hz. Thanks for your help, guys!

Oscilloscope reading

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.