|
| 1 | +# acq and init taken from: https://intantech.com/files/RHX_TCP.zip |
| 2 | +# In order to run this example script successfully, the Intan RHX software |
| 3 | +# should first be started, and through Network -> Remote TCP Control: |
| 4 | +# Command Output should open a connection at 127.0.0.1, Port 5000. |
| 5 | +# Status should read "Pending" |
| 6 | +# Waveform Output (in the Data Output tab) should open a connection at 127.0.0.1, Port 5001. |
| 7 | +# Status should read "Pending" for the Waveform Port (Spike Port is unused for this example, |
| 8 | +# and can be left disconnected) |
| 9 | +# Once these ports are opened, this script can be run to acquire ~1 second of wideband data from channel A-010, |
| 10 | + |
| 11 | +import concore |
| 12 | +import time |
| 13 | +import threading |
| 14 | +import socket |
| 15 | +import matplotlib.pyplot as plt |
| 16 | + |
| 17 | +try: |
| 18 | + showPlot = concore.params['plot'] |
| 19 | +except: |
| 20 | + showPlot = 0 |
| 21 | + |
| 22 | +try: |
| 23 | + tsamp = concore.params['tsamp'] |
| 24 | + if tsamp > 1: |
| 25 | + tsamp = 1 |
| 26 | +except: |
| 27 | + tsamp = 1 |
| 28 | + |
| 29 | +def validseq(t): |
| 30 | + validindex = 0 |
| 31 | + for i in range(0,len(t)-1): |
| 32 | + if t[i] >= t[i+1]: |
| 33 | + print('i='+str(i)+': '+str(t[i])+'>='+str(t[i+1])+' len='+str(len(t))) |
| 34 | + validindex = i+1 |
| 35 | + return validindex |
| 36 | + |
| 37 | +def extract(amplifierData): |
| 38 | + s = 0.0 |
| 39 | + for data in amplifierData: |
| 40 | + s += data |
| 41 | + #time.sleep(0.5/len(amplifierData)) #make extract really slow for testing |
| 42 | + return s/len(amplifierData) |
| 43 | + |
| 44 | +def readUint32(array, arrayIndex): |
| 45 | + variableBytes = array[arrayIndex : arrayIndex + 4] |
| 46 | + variable = int.from_bytes(variableBytes, byteorder='little', signed=False) |
| 47 | + arrayIndex = arrayIndex + 4 |
| 48 | + return variable, arrayIndex |
| 49 | + |
| 50 | +def readInt32(array, arrayIndex): |
| 51 | + variableBytes = array[arrayIndex : arrayIndex + 4] |
| 52 | + variable = int.from_bytes(variableBytes, byteorder='little', signed=True) |
| 53 | + arrayIndex = arrayIndex + 4 |
| 54 | + return variable, arrayIndex |
| 55 | + |
| 56 | +def readUint16(array, arrayIndex): |
| 57 | + variableBytes = array[arrayIndex : arrayIndex + 2] |
| 58 | + variable = int.from_bytes(variableBytes, byteorder='little', signed=False) |
| 59 | + arrayIndex = arrayIndex + 2 |
| 60 | + return variable, arrayIndex |
| 61 | + |
| 62 | +def clearall(): |
| 63 | + # Clear TCP data output to ensure no TCP channels are enabled |
| 64 | + scommand.sendall(b'execute clearalldataoutputs') |
| 65 | + time.sleep(0.1) |
| 66 | + # Send TCP commands to set up TCP Data Output Enabled for wide |
| 67 | + # band of channel A-010 |
| 68 | + scommand.sendall(b'set a-010.tcpdataoutputenabled true') |
| 69 | + time.sleep(0.1) |
| 70 | + |
| 71 | +def acq(): |
| 72 | + # Run controller for tsamp second |
| 73 | + scommand.sendall(b'set runmode run') |
| 74 | + time.sleep(tsamp) |
| 75 | + scommand.sendall(b'set runmode stop') |
| 76 | + |
| 77 | + # Read waveform data |
| 78 | + rawData = swaveform.recv(WAVEFORM_BUFFER_SIZE) |
| 79 | + if len(rawData) % waveformBytesPerBlock != 0: |
| 80 | + raise Exception('unexpected amount of data') |
| 81 | + numBlocks = int(len(rawData) / waveformBytesPerBlock) |
| 82 | + |
| 83 | + rawIndex = 0 # Index used to read the raw data |
| 84 | + |
| 85 | + for block in range(numBlocks): |
| 86 | + # Expect 4 bytes to be TCP Magic Number as uint32. |
| 87 | + # If not what's expected, raise an exception. |
| 88 | + magicNumber, rawIndex = readUint32(rawData, rawIndex) |
| 89 | + if magicNumber != 0x2ef07a08: |
| 90 | + raise Exception('Error... magic number incorrect') |
| 91 | + |
| 92 | + # Each block should contain 128 frames of data - process each |
| 93 | + # of these one-by-one |
| 94 | + for frame in range(framesPerBlock): |
| 95 | + # Expect 4 bytes to be timestamp as int32. |
| 96 | + rawTimestamp, rawIndex = readInt32(rawData, rawIndex) |
| 97 | + |
| 98 | + # Multiply by 'timestep' to convert timestamp to seconds |
| 99 | + amplifierTimestamps.append(rawTimestamp * timestep) |
| 100 | + |
| 101 | + # Expect 2 bytes of wideband data. |
| 102 | + rawSample, rawIndex = readUint16(rawData, rawIndex) |
| 103 | + |
| 104 | + # Scale this sample to convert to microVolts |
| 105 | + amplifierData.append(0.195 * (rawSample - 32768)) |
| 106 | + |
| 107 | + |
| 108 | +#initialization |
| 109 | +# Declare buffer size for reading from TCP command socket |
| 110 | +# This is the maximum number of bytes expected for 1 read. 1024 is plenty for a single text command |
| 111 | +COMMAND_BUFFER_SIZE = 1024 # Increase if many return commands are expected |
| 112 | +# Declare buffer size for reading from TCP waveform socket. |
| 113 | +# This is the maximum number of bytes expected for 1 read |
| 114 | +# There will be some TCP lag in both starting and stopping acquisition, so the exact number of data blocks may vary slightly. |
| 115 | +# At 30 kHz with 1 channel, 1 second of wideband waveform data is 181,420 byte. See 'Calculations for accurate parsing' for more details |
| 116 | +# To allow for some TCP lag in stopping acquisition resulting in slightly more than 1 second of data, 200000 should be a safe buffer size |
| 117 | +WAVEFORM_BUFFER_SIZE = 200000 # Increase if channels, filter bands, or acquisition time increase |
| 118 | +# Connect to TCP command server - default home IP address at port 5000 |
| 119 | +print('Connecting to TCP command server...') |
| 120 | +scommand = socket.socket(socket.AF_INET, socket.SOCK_STREAM) |
| 121 | +scommand.connect(('127.0.0.1', 5000)) |
| 122 | +# Connect to TCP waveform server - default home IP address at port 5001 |
| 123 | +print('Connecting to TCP waveform server...') |
| 124 | +swaveform = socket.socket(socket.AF_INET, socket.SOCK_STREAM) |
| 125 | +swaveform.connect(('127.0.0.1', 5001)) |
| 126 | +# Query runmode from RHX software |
| 127 | +scommand.sendall(b'get runmode') |
| 128 | +commandReturn = str(scommand.recv(COMMAND_BUFFER_SIZE), "utf-8") |
| 129 | +isStopped = commandReturn == "Return: RunMode Stop" |
| 130 | +# If controller is running, stop it |
| 131 | +if not isStopped: |
| 132 | + scommand.sendall(b'set runmode stop') |
| 133 | + time.sleep(0.1) # Allow time for RHX software to accept this command before the next one comes |
| 134 | + |
| 135 | +scommand.sendall(b'set TCPNumberDataBlocksPerWrite 20') |
| 136 | +time.sleep(0.1) |
| 137 | + |
| 138 | +# Query sample rate from RHX software |
| 139 | +scommand.sendall(b'get sampleratehertz') |
| 140 | +commandReturn = str(scommand.recv(COMMAND_BUFFER_SIZE), "utf-8") |
| 141 | +expectedReturnString = "Return: SampleRateHertz " |
| 142 | +if commandReturn.find(expectedReturnString) == -1: # Look for "Return: SampleRateHertz N" where N is the sample rate |
| 143 | + raise Exception('Unable to get sample rate from server') |
| 144 | +else: |
| 145 | + sampleRate = float(commandReturn[len(expectedReturnString):]) |
| 146 | + |
| 147 | +# Calculate timestep from sample rate |
| 148 | +timestep = 1 / sampleRate |
| 149 | + |
| 150 | + |
| 151 | +# Calculations for accurate parsing |
| 152 | +# At 30 kHz with 1 channel, 1 second of wideband waveform data (including magic number, timestamps, and amplifier data) is 181,420 bytes |
| 153 | +# N = (framesPerBlock * waveformBytesPerFrame + SizeOfMagicNumber) * NumBlocks where: |
| 154 | +# framesPerBlock = 128 ; standard data block size used by Intan |
| 155 | +# waveformBytesPerFrame = SizeOfTimestamp + SizeOfSample ; timestamp is a 4-byte (32-bit) int, and amplifier sample is a 2-byte (16-bit) unsigned int |
| 156 | +# SizeOfMagicNumber = 4; Magic number is a 4-byte (32-bit) unsigned int |
| 157 | +# NumBlocks = NumFrames / framesPerBlock ; At 30 kHz, 1 second of data has 30000 frames. NumBlocks must be an integer value, so round up to 235 |
| 158 | + |
| 159 | +framesPerBlock = 128 |
| 160 | +waveformBytesPerFrame = 4 + 2 |
| 161 | +waveformBytesPerBlock = framesPerBlock * waveformBytesPerFrame + 4 |
| 162 | + |
| 163 | +clearall() |
| 164 | + |
| 165 | +concore.delay = 0.01 |
| 166 | +init_simtime_u = "[0.0, 0.0, 0.0]" |
| 167 | +init_simtime_ym = "[0.0, 0.0, 0.0]" |
| 168 | + |
| 169 | +ym = concore.initval(init_simtime_ym) |
| 170 | + |
| 171 | +acq_thread = threading.Thread(target=acq, daemon=True) |
| 172 | +nxtacq_thread = threading.Thread(target=acq, daemon=True) |
| 173 | +amplifierTimestamps = [] |
| 174 | +amplifierData = [] |
| 175 | +acq_thread.start() |
| 176 | +while(concore.simtime<concore.maxtime): |
| 177 | + while concore.unchanged(): |
| 178 | + ym = concore.read(1,"ym",init_simtime_ym) |
| 179 | + acq_thread.join() |
| 180 | + acq_thread = nxtacq_thread |
| 181 | + oldAmpT = amplifierTimestamps |
| 182 | + oldAmpD = amplifierData |
| 183 | + amplifierTimestamps = [] |
| 184 | + amplifierData = [] |
| 185 | + if showPlot&1==1: |
| 186 | + plt.plot(oldAmpT,oldAmpD) |
| 187 | + plt.title("ym="+str(extract(oldAmpD))) |
| 188 | + plt.xlabel('Time (s)') |
| 189 | + plt.ylabel('Voltage (uV)') |
| 190 | + plt.show() |
| 191 | + if showPlot&2==2: |
| 192 | + plt.plot(range(0,len(oldAmpD)),oldAmpD) |
| 193 | + plt.title("ym="+str(extract(oldAmpD))) |
| 194 | + plt.xlabel('index') |
| 195 | + plt.ylabel('Voltage (uV)') |
| 196 | + plt.show() |
| 197 | + acq_thread.start() |
| 198 | + nxtacq_thread = threading.Thread(target=acq, daemon=True) |
| 199 | + validindex = validseq(oldAmpT) |
| 200 | + if validindex == 0: |
| 201 | + if oldAmpT[0] != 0: |
| 202 | + print("validseq yet starttime="+str(oldAmpT[0])) |
| 203 | + else: |
| 204 | + if oldAmpT[validindex] != 0: |
| 205 | + print(str(validindex)+" starttime="+str(oldAmpT[validindex])) |
| 206 | + if showPlot&4==4: |
| 207 | + plt.plot(oldAmpT,oldAmpD) |
| 208 | + plt.title("ym="+str(extract(oldAmpD))+" len="+str(len(oldAmpD))) |
| 209 | + plt.xlabel('Time (s)') |
| 210 | + plt.ylabel('Voltage (uV)') |
| 211 | + plt.show() |
| 212 | + if showPlot&8==8: |
| 213 | + plt.plot(range(0,len(oldAmpD)),oldAmpD) |
| 214 | + plt.title("ym="+str(extract(oldAmpD))+" len="+str(len(oldAmpD))) |
| 215 | + plt.xlabel('index') |
| 216 | + plt.ylabel('Voltage (uV)') |
| 217 | + plt.show() |
| 218 | + oldAmpT = oldAmpT[validindex:] |
| 219 | + oldAmpD = oldAmpD[validindex:] |
| 220 | + if showPlot&4==4: |
| 221 | + plt.plot(oldAmpT,oldAmpD) |
| 222 | + plt.title("ym="+str(extract(oldAmpD))+" len="+str(len(oldAmpD))) |
| 223 | + plt.xlabel('Time (s)') |
| 224 | + plt.ylabel('Voltage (uV)') |
| 225 | + plt.show() |
| 226 | + if showPlot&8==8: |
| 227 | + plt.plot(range(0,len(oldAmpD)),oldAmpD) |
| 228 | + plt.title("ym="+str(extract(oldAmpD))+" len="+str(len(oldAmpD))) |
| 229 | + plt.xlabel('index') |
| 230 | + plt.ylabel('Voltage (uV)') |
| 231 | + plt.show() |
| 232 | + ym[0] = extract(oldAmpD) |
| 233 | + print("ym="+str(ym[0])) |
| 234 | + concore.write(1,"ym",ym) |
| 235 | +acq_thread.join() |
| 236 | +print("retry="+str(concore.retrycount)) |
0 commit comments