SIGNAL PROCESSING SYSTEM FOR BRAIN MACHINE INTERFACE AND METHOD PERFORMED BY THE SYSTEM

Abstract

Disclosed is a signal processing system for a brain machine interface (BMI) and a method performed by such a signal processing system. The method involves receiving and filtering input data comprising cortical signal samples, calculating a threshold for detection of spikes, identifying spikes in subsequently received samples, and outputting spike data corresponding to the identified spike or spikes.

Claims

1. A signal processing method on a signal processing system in a brain machine interface (BMI), comprising the processing steps of: receiving input data comprising cortical signal samples; applying a filter to the input data to obtain filtered data samples; calculating a threshold by: initializing one or more parameters; recursively updating the one or more parameters based on a value of each filtered data sample; grouping each filtered data sample into an ordered set of groups, based on the value of the respective filtered data sample and the recursively updated one or more parameters; identifying a threshold group based on a number of filtered data samples grouped into groups distributed, around the threshold group, in the ordered set of groups; and determining the threshold based on the threshold group and the one or more parameters; identifying a spike in each of one or more subsequently received data samples based on the threshold; and outputting spike data corresponding to the spike or spikes identified in the one or more subsequently received data samples.

2. A method according to claim 1, wherein the input data is measured by a sensor in a cortex of a patient.

3. A method according to claim 2, wherein the sensor is an implanted multi-channel sensor array, the signal processing system is connected to the implanted multi-channel sensor array and the processing steps are performed on cortical signal samples captured on each respective channel of the multi-channel sensor array.

4. A method according to claim 1, wherein initializing the one or more parameters comprises setting: a maximum value (iMax); a minimum value (iMin); and a number of groups for the ordered set of groups.

5. A method according to claim 3, wherein recursively updating the one or more parameters based on a value of each filtered data sample comprises setting at least one of iMax and iMin to be, respectively, a maximum value and a minimum value across all of the filtered data samples.

6. A method according to claim 4, further comprising calculating a dynamic range of the filtered data samples (iRng) based on the formula:
iRng=iMaxiMin+1

7. A method according to claim 6, wherein grouping each filtered data sample into an ordered set of groups comprises identifying, for each filtered data sample, a respective group of the ordered set of groups according to the formula: $\mathrm{iB}=\frac{x\left(k\right)-\mathrm{iMin}}{\mathrm{iRng}}\mathrm{nBin}$ wherein iB is the bin index for the respective group, x(k) is a filtered data sample comprising k seconds of the input data, and nBin is the number of groups in the ordered set of groups.

8. A method according to claim 1, wherein applying a filter to the input data to obtain filtered data samples comprises applying a Finite Impulse Response (FIR) filter or an Infinite Impulse Response (IIR) filter to the input data.

9. A method according to claim 8, wherein applying a filter to the input data to obtain filtered data samples comprises applying an IIR filter, and wherein the input data has a first precision in which the cortical signal samples are each represented by a first number of bits, and wherein applying an IIR filter to the input data comprises applying a modified Direct Form I filter of order N, having 2N delay units, involving: scaling upscaling filter coefficients b and scaling downscaling filter coefficients a to obtain scaled filter coefficients b and a respectively; casting an input signal x(n) and the scaled filter coefficients into a second precision, wherein x(n) comprises the input data and the second precision is higher than the first precision; applying the IIR filter to the input signal x(n), using b and a as coefficients for the IIR filter, to obtain an intermediate result; and downscaling the intermediate result into the first precision.

10. A method according to claim 9, wherein scaling upscaling filter coefficients b and scaling downscaling filter coefficients a to obtain scaled filter coefficients {tilde over (b)} and comprises: applying the IIR filter to a stored set of cortical signal samples to produce a first set of filtered samples; applying a floating-point filter to the stored set of cortical signal samples to produce a second set of filtered samples; determining a number of precision bits n.sub.p resulting in a minimal difference between the first set of filtered samples and the second set of filtered samples, and scaling to 2.sup.n.sup.p to obtain the scaled filter coefficients {tilde over (b)} and .

11. A method according to claim 9, wherein the intermediate result has a first form y(n) and a second form, and wherein applying the IIR filter to the input signal x(n), using {tilde over (b)} and as coefficients for the IIR filter, obtains y(n), the method further comprising: obtaining the second form of the intermediate result by retaining a sign of y(n), casting y(n) to an unsigned integer format, wherein downscaling the intermediate result into the first precision comprises downscaling the second form of the intermediate result into the first precision.

12. A method according to claim 11, wherein the second form is SY.sub.f, which is obtained from y(n) by: retaining the sign S of y(n); casting y(n) to an unsigned integer format Y.sub.f; and bit-shifting Y.sub.f to obtain SY.sub.f.

13. A signal processing system for a brain machine interface (BMI), the signal processing system comprising: a receiver for receiving cortical signals measured by a multi-channel sensor array; a wireless transmitter for transmitting an output; memory; and at least one processor, the memory storing instructions that, when executed by the at least one processor, cause the at least one processor to, for each channel of the multi-channel sensor array: receive input data at the receiver, the input data comprising cortical signal samples; apply a filter to the input data to obtain filtered data samples; calculate a threshold by: initializing one or more parameters; recursively updating the one or more parameters based on a value of each filtered data sample; grouping each filtered data sample into an ordered set of groups, based on the value of the respective filtered data sample and the recursively updated one or more parameters; identifying a threshold group based on a number of filtered data samples grouped into groups distributed, around the threshold group, in the ordered set of groups; and determining the threshold based on the threshold group and the one or more parameters; and identify a spike in each of one or more subsequently received data samples, received at the receiver, based on the threshold, the instructions, when executed by the at least one processor, further causing the at least one processor to output the output, from the wireless transmitter to a remote system, the output comprising spike data corresponding to the spike or spikes identified in the one or more subsequently received data samples.

14. A signal processing system according to claim 13, wherein the at least one processor is caused to apply a filter to the input data to obtain filtered data samples by applying a Finite Impulse Response (FIR) filter or an Infinite Impulse Response (IIR) filter to the input data.

15. A signal processing system according to claim 14, wherein applying a filter to the input data to obtain filtered data samples comprises applying an IIR filter, and wherein the input data has a first precision in which the cortical signal samples are each represented by a first number of bits, and the at least one processor is caused to apply an IIR filter to the input data by applying a modified Direct Form I filter of order N, having 2N delay units, by: scaling upscaling filter coefficients b and scaling downscaling filter coefficients a to obtain scaled filter coefficients b and a respectively; casting an input signal x(n) and the scaled filter coefficients into a second precision, wherein x(n) comprises the input data and the second precision is higher than the first precision; applying the IIR filter to the input signal x(n), using b and a as coefficients for the IIR filter, to obtain an intermediate result; and downscaling the intermediate result into the first precision.

16. A signal processing system according to claim 15, wherein the at least one processor is caused to scale upscaling filter coefficients b and scale downscaling filter coefficients a to obtain scaled filter coefficients {tilde over (b)} and by: applying the IIR filter to a stored set of cortical signal samples to produce a first set of filtered samples; applying a floating-point filter to the stored set of cortical signal samples to produce a second set of filtered samples; determining a number of precision bits n.sub.p resulting in a minimal difference between the first set of filtered samples and the second set of filtered samples, and scaling to 2.sup.n.sup.p to obtain the scaled filter coefficients {tilde over (b)} and .

17. A signal processing system according to claim 15, wherein the intermediate result has a first form y(n) and a second form, and wherein the at least one processor is caused to: apply the IIR filter to the input signal x(n), using {tilde over (b)} and as coefficients for the IIR filter, to obtain y(n); and obtain the second form of the intermediate result by retaining a sign of y(n), casting y(n) to an unsigned integer format, wherein downscaling the intermediate result into the first precision comprises downscaling the second form of the intermediate result into the first precision.

18. A signal processing system according to claim 17, wherein the second form is SY.sub.f, and the at least one processor obtains SY.sub.f from y(n) by: retaining the sign S of y(n); casting y(n) to an unsigned integer format Y.sub.f; and bit-shifting Y.sub.f to obtain SY.sub.f.

19. A signal processing system according to claim 18, wherein the at least one processor is caused to bit-shift Y.sub.f by left shifting Y.sub.f to scale by 2.sup.n.sup.r.sup.n.sup.p, wherein n.sub.r is a total number of bits used to represent x(n).

20. A brain machine interface comprising a signal processing system according to claim 13.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0122] Preferred embodiments of the invention are hereafter described, by way of non-limiting example only, with reference to the accompanying drawings, in which:

[0123] FIG. 1 is a schematic diagram of a method for DSP on a signal processing system, with the steps of the method being performed in an exemplary order;

[0124] FIG. 2 is an embodiment of a system for performing the method of FIG. 1;

[0125] FIG. 3 is a schematic block diagram of an embodiment of an overall system integrating a signal processing engine for signal processing system in accordance with present embodiments with the decoding/classification engine;

[0126] FIG. 4 illustrates a precision-preserving fixed-point filter with overflow/underflow handling (PPFF-FwHd) in accordance with method 100;

[0127] FIG. 5 is an illustration of the bit-shifting process performed during filtering according to FIG. 4;

[0128] FIG. 6 is a block diagram illustrating a threshold calculation method in accordance with method 100, specifically a RVCT-aTrFx method;

[0129] FIG. 7 illustrates spike detection in accordance with method 100; and

[0130] FIGS. 8 and 9 show comparisons of the filtered data of a channel and thresholds obtained.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0131] The following detailed description sets out embodiments of a signal processing system (which may be mounted or attached to a scalp of the subject (e.g. a patient) and connect with a multi-channel (e.g. microelectrode) array implanted in the motor cortex of the subject, or the signal processing system may itself be implantablee.g. with the microelectrode array), system (which may be, or comprise, the signal processing system) and methods executed by the signal processing system that employ a fixed-point implementation to compute neural spike features from cortical signals to facilitate DSP implementation in a brain machine interface (BMI).

[0132] The purpose of some embodiments is to enable the low-bandwidth, wireless transmission of the spike firing rate for BMI decoding. This can allow the user to regain free-range movements with the removal of the wired connection between the brain and the processing computer. Unlike existing approaches that require offline training to obtain a threshold or other parameters for online processing, the methods proposed herein compute the threshold by dynamic value-count tracking, which allows for adaptive thresholding by iteratively updating a buffer as described below. The threshold computation is based on recursive value-count tracking that approximates the median-based threshold with proven robustness and non-sensitivity to noise. The final spike detection is triggered by the crossing-transitions to reduce the computational power. This method possesses the lowest, or very low, computational load making it suitable for a signal processing system in, for example, an on-chip implementation.

[0133] In order to transmit the neural spike signals from the brain cortex to a computer wirelessly, fixed-point DSP has been found to be easy to integrate with applications. However, the overflow and underflow, and round-off noise has to be considered for each operation of fixed-point implementation. Most existing methods for signal processing systems, particularly those that are mounted to the patient or are otherwise on-chip implementations, require an offline training process to obtain the threshold for online processing. As taught herein, novel methods are proposed for performing fully online signal processing and wireless transmission of the neural spike data. In particular, proposed is a recursive value-count tracking-based method to calculate the threshold adaptively, an efficient spike detection method triggered by crossing-transition, and a precision-preserving fixed-point filter that is able to handle the overflow/underflow automatically with the preserved bits.

[0134] FIG. 1 is a schematic diagram of a signal processing method 100 for a signal processing system in a BMI. The method 100 broadly comprises the processing steps of: [0135] 102: receiving input data comprising cortical signal samples; [0136] 104: applying a filter to the input data to obtain filtered data samples; [0137] 106: calculating a threshold; [0138] 108: identifying a spike in each of one or more subsequently received data samples based on the threshold; and [0139] 110: outputting spike data corresponding to the spike or spikes identified in the one or more subsequently received data samples.

[0140] The method 100 may be performed on a system 200, as shown in FIG. 2, including a signal processing system 202 (which may be a chip), e.g. for DSP. The signal processing system 202 comprises a transmitter 206 for transmitting an output, memory 208, at least one processor (for ease, hereinafter the processor 210) and a receiver 214 connectable to a multi-channel array 204 for sensing cortical signals. The memory 208 stores instructions that, when executed by the processor 210, cause the processor 210 to perform the method using the array 204, transmitter 206 and memory 208, receiver 214 for each channel of the multi-channel array 204. The memory 208, or a different memory or database, may also store cortical signals. The processing steps of method 100 are performed on cortical signal samples captured on each respective channel of the multi-channel array 204. The signal processing system 202 may incorporate the array 204 but, in the present embodiment, is connected to the array 204 via a hardwired connection or a wireless connection as indicated by the broken lines.

[0141] As will be appreciated by the skilled person, wired or wireless transmission may be achieved using any suitable method and protocol. As mentioned above, the connection between system 212 and array 204 may be hardwirede.g. by a wire passing through the scalp to the array 204 implanted in the cortex of the subject, or where the system 202 and array 204 form part of the same device they may both be implanted on the cortex together to facilitate communication of measured cortical signalsor by wireless connectione.g. by providing an antenna, or micro antenna, on the array 204 and an antenna of the system 202, the antenna or micro antenna on the array potentially only being part of a transmitter rather than a transceiver. The system 202 may similarly communicate with the remote system 212 wirelessly, using the same or different antenna to that used to communicate wirelessly with the array 204, or in a hardwired configurationthe latter is less preferable since it requires the subject to be near the remote system 212 during signal decoding. The skilled person will understand such arrangements and variations thereon, all of which are intended to form part of the present disclosure.

[0142] The system 200 of the embodiment shown in FIG. 2 also includes a spike decoding system (i.e. remote system 212), in wired or wireless communication with the transmitter 206 (thus making the transmitter a wired or wireless transmitter, respectively), for receiving spike data from the transmitter 206 and decoding the spike data to determine control commands.

[0143] Many of the methods provided in accordance with method 100, for implementation by system 200, in BMI decoding are described below and consist of methods for recursive value-count tracking-based adaptive threshold calculation for fixed-point implementation (RVCT-aTrFx), crossing-transition triggered spike detection (CTT-SPD) and precision-preserving fixed-point filter with automatic overflow/underflow handling (PPFF-FwHd). The objective as described above is to facilitate the low-bandwidth wireless transmission of the neural spike features from the brain cortex to a computer for online decoding to control an assistive device. Such implementation is important to allow the users to be untethered. As a result, the subjects can be relieved from restricted movements, infection, and other noise and electrical interference.

[0144] The system 200 employs a method for calculating the threshold by dynamically tracking value-countsee block 312 of FIG. 3 and FIG. 6of the cortical signal samples for the purpose of detecting spikes in BMI for the fixed-point implementation.

[0145] As will be described in more detail, the important features in threshold calculation method are: [0146] 1) Recursively tracking of the dynamic ranges of samples so that an efficient value-count based representation can be established. The bin or group number is selected to store the samples based on the dynamic ranges and maximum values of the samples. [0147] 2) Computing the basic threshold (i.e., estimation of the standard deviation of background noise) based on the distributions of the samples in the value-count profile. [0148] 3) Calculating the threshold adaptively at a predefined time interval with buffer management. In some instances, adaptive calculation may not be required.

[0149] As described with reference to block 312 of FIG. 3 and the exemplary embodiment in FIG. 6, the above steps were achieved by: [0150] a) Dynamically tracking the samples to find the ranges of the values of the samples by iteratively updating the minimum and maximum, till all the samples were processed (e.g. taking each sample as a single data point and identifying a maximum and minimum value of the data points for a given set of samples). [0151] b) Recursively tracking the value-count by throwing the samples into the proper group or bin based on the value of the given sample. The number of samples contained in each bin were dynamically added to the accumulated counter, until the counter value was larger than the defined distribution criterione.g. the median (e.g., 50% of total number of samples)and the relevant bin (or its bin index) was selected. [0152] c) Taking the corresponding value of the selected bin index as the basic threshold. [0153] d) Scaling the basic threshold based on the parameters given which can be obtained from offline or asynchronous training based on history data to optimize the classification accuracy. [0154] e) Calculating the threshold based on the latest k seconds of data in the buffer which was refreshed for every n seconds, where k and n were determined based on the buffer condition and processing speed requirements.

[0155] The system 200 also employs a method for detecting spikes, detection being triggered by crossing-transition as described with reference to block 314 of FIG. 3, and FIG. 7. As will be described in more detail below, the important features for this method are: [0156] 1) Defining a trigger criterion to detect the spikes only when a crossing-transition exists in the processing sample and its previous immediate neighbouring sample (PV-INS). This significantly reduces the computational complexity since the crossings only occupy a small portion of the original samples. [0157] 2) Defining a local fluctuation measure to identify the entering point of the edge of spikes.

[0158] In one exemplary embodiment, this process involved: [0159] a) Identifying a crossing when the sum of the sample (i.e. sample value) and the threshold was less than zero. [0160] b) Examining its PV-INS when a crossing transition was detected (i.e., only one of the sample and its PV-INS was a crossing). If the value of PV-INS was larger than the processing sample, an entering point of the edge of a spike was detected, hence a spike was detected.

[0161] The system 200 also employs a method for filtering the raw cortical signal based on fixed-point computation to facilitate DSP implementationsee block 310 of FIG. 3, and FIGS. 4 and 5. The important features of this filtering process are: [0162] 1) Choosing the precision bit used to convert the filter coefficients and represent the input data (n.sub.p) by minimizing the difference between the filtered data obtained from fixed-point filter and that from floating-point filter. [0163] 2) Preserving the data precision by double-scaling technique. The data samples (through multiplying by the filter coefficients) were firstly scaled up to a higher precision and the signs were protected. After filtering, the intermediate results were scaled down to the original data precision for processing the next sample. [0164] 3) Automatically handling the underflow and overflow using the preserved bits, i.e., the difference between the determined precision bit (n.sub.p), and the actual bits representing the input data (n.sub.a), generally n.sub.p>n.sub.a.

[0165] Importantly, the proposed system and method can be integrated with any decoding/classification method by transmitting the computed spike features (spike data) wirelessly to a computer where the engine of decoding/classification and control of prosthetic device resides. Such low-bandwidth features will allow efficient decoding the neural spiking signal in real time to control an assistive device.

[0166] The purpose of the system 200 and method 100 is to enable the low-bandwidth, wireless transmission of the spike firing rate for BMI decoding. This will allow the user to regain free-range movements with the removal of the wired connection between the brain and the processing computer. Unlike most existing approaches that require offline training to obtain a threshold or other parameters for online processing, our proposed method computed the threshold by dynamically value-count tracking, which allows for adaptive thresholding by iteratively updating the buffer. The threshold computation was based on recursive value-count tracking (tracking the number of samples with each respective value from a set of all possible values) which approximated the median-based threshold with proved robustness and non-sensitiveness to noise. The final spike detection was triggered by the crossing-transitions to reduce the computational power. This method possessed the lowest computational load making it suitable for on-chip implementation.

Overall System

[0167] An example overall system 300 integrating the present method 100 of signal processing for, for example, DSP chip implementation with the decoding/classification engine, to control mobile robots such as controlling a wheelchair is shown in FIG. 3. The system 300 is represented as functional blocks that can be implemented by the memory 208 (e.g. incorporating a buffer), transmitter 206, array 204 and processor 210. The decoding can ultimately be performed by remote system 212, corresponding to computer 318.

[0168] FIG. 3 illustrates the proposed overall system integrating the presently proposed signal processing engine (e.g. for DSP chip implementation) with the decoding/classification engine stored in a remote computer 318such as remote system 212.

[0169] The raw data (i.e. measurements) obtained from the multi-channel array 204, presently microwire electrode arrays 302 (i.e., data recording device) were amplified and the amplified raw data were then analog-to-digital converted at block 304, and subsequently multiplexed at block 306. Hence, proper de-multiplexing techniques were employed to obtain data channel by channel. The amplification and multiplexing blocks 304, 306 may be on the same chip or device as the multi-channel array 204, and/or may constitute components of signal processing system 202 (e.g. in the same housing as signal processing system 202).

[0170] The initial data (filtered data samples) for calculating the threshold (i.e., Dt-Tr) were extracted from the beginning of data (e.g., 0.5 second of data (referred to with reference to FIG. 6 as k seconds of data) by skipping the unstable part at the very beginning) and buffered in a buffer. Dt-Tr can be iteratively updated with more data (i.e. subsequently received cortical signal samples).

[0171] The buffer can be refreshed with new data after a predefined time interval, i.e., by the buffer management block. Hence, the threshold calculation can operate in adaptive mode.

[0172] Subsequent data were sequentially read to calculate firing rate (Dt-FR). Thus, while method 100 appears to propose a specific sequence of steps, those steps may be skipped or performed as needs bee.g. the threshold need not be calculated for every set of input data.

[0173] The data was pre-processed at block 308 before the filtering process was executed at block 310 (corresponding to FIG. 4). Pre-processing may involve any appropriate data manipulation, such as centering the data to remove DC drift (i.e., the mean of the data was not centred at zero), which can be done by deducting the mean of data for the chosen time segment.

[0174] The pre-processed data outputted by block 308 were then cast to the format with the same precision as the input (denoted as X.sub.cTr). X.sub.cTr were filtered with a precision-preserving fixed-point filter with overflow/underflow handling, discussed with reference to block 310 and FIG. 4, to obtain X.sub.fTr. Thresholds were calculated based on a recursive value-count tracking-based adaptive threshold calculation method, discussed with reference to block 312 and FIG. 6, to obtain T.sub.r.

[0175] The same procedure was carried out for subsequent data to calculate the firing rate, i.e., Dt-FR. The data were pre-processed the data to obtain X.sub.cFr and filtered the data to obtain X.sub.fFr. The spikes were detected based on a crossing-transition triggered spike detection (CTT-SPD) method discussed with reference to block 314 and FIG. 3. A certain decoding bin width (e.g. 20 ms, 50 ms, 100 ms or otherwise as desired) was employed to refresh the counter timely, e.g., every 100 ms to output the firing rate (FR) for all the channels. Note that the processed samples can be adaptively buffered and used to calculate the mean (or other distribution characteristic as discussed with reference to FIG. 6) of the data.

[0176] The output FR can be wirelessly transmitted (316) to a computer (212, 318) where a decoding or classification engine resides. The final decoding results can then be used to control a mobile robot such as wheelchair (320) to allow the user to achieve mobility.

Embodiment of Step 104Precision-Preserving Fixed-Point Filter with Overflow/Underflow Handling (PPFF-FwHd)

[0177] The main concept behind the present filter design for BMI is to keep the computational complexity low to facilitate signal processing system implementation, such as on-chip implementation, while preserving precision of data. An infinite impulse response (IIR) filter was chosen due to the fact it could meet given specifications with much lower order compared with that of finite impulse response (FIR) filters which may also be used in some implementations. This may significantly conserve the number of computations required per time step.

[0178] A disadvantage of IIR filters is the nonlinear phase. However, the non-linear phase may be eliminated by offline (or asynchronous) implementation with non-causal, zero-phase filtering.

[0179] The filter 400 proposed in FIG. 4 employs a modified Direct Form I (MDF-I) design, with 2N delay units for both input and output signals for a filter of order N. The filter 400 also comprises a Shifter engine that operates directly on the input filter coefficients and intermediate results.

[0180] The precision bits (n.sub.p) were firstly determined offline (or asynchronously) based on historical data. The determination was made to minimize the errors or differences between the filtered data obtained from using a fixed-point filter and that of floating-point filter.

[0181] Typically n.sub.p<n.sub.r (the unused bits n.sub.ub=n.sub.rn.sub.p), where n.sub.r is the total number of bits used to represent input data, e.g., 16 bits. However, the number of precision bits (n.sub.p) is typically larger than that of the actual precision of input data (n.sub.a), e.g., a first number of bit such as 10 bits. The preserved bits may therefore be denoted as n.sub.pv=n.sub.pn.sub.a. The preserved bits can be used to automatically handle the underflow and overflow problem, where the summation would otherwise result in a number outside the range capable of being specified by the number of precision bits.

[0182] A double-scaling process was then employed. The first scaling process is to scale the filter coefficients up to the precision bits (i.e., multiply by 2.sup.n.sup.p). This will only be performed once before the filtering process. The second scaling process is firstly to scale the intermediate results up to the power of 2.sup.n.sup.r.sup.n.sup.p and scale down to the total number of bits used to represent the input data (n.sub.r). Hence, overall the data are scaled down to 2.sup.n.sup.p. This ensures that the input data and output filtered results are always of the same precision. Thus precision is preserved.

[0183] The proposed method 400 has therefore been named precision-preserving fixed-point filter with overflow/underflow handling (PPFF-FwHd) for DSP implementation in BMI. This is illustrated in FIG. 4. In the present example, the filtering process was carried out in a higher precision buffer (e.g., int32) to avoid overflow as well.

Detailed Filtering Process

[0184] The PPFF-FwHd may perform the following steps:

[0185] Upscaling filter coefficients b and downscaling filter coefficients a were scaled to obtain scaled filter coefficients {tilde over (b)} and respectively (402). The filter coefficients (a and b) were thus scaled to 2.sup.n.sup.p to obtain and {tilde over (b)}. and {tilde over (b)} were then cast to the same precision as the input data, e.g., int16. In some embodiments, this step can be considered part of the filter 400, though it has been performed in advance of the steps represented schematically in FIG. 4the proposed precision-preserving fixed-point filter with overflow/underflow handling (PPFF-FwHd).

[0186] The input signal x(n) and the scaled filter coefficients are cast into a second precision (note the precision of the input signal is specified above as being of a first precision), wherein x(n) comprises the input data received at step 102. Moreover, the second precision is higher than the first precision. In other words, the input signal (x(n)) and scaled filter coefficients ( and {tilde over (b)}) were cast into a higher precision of integer format (e.g., int32) to avoid overflow during computation.

[0187] A buffer (Y) is then initialised to store the intermediate results with size y=zeros(1,length(B.sub.t)) where B.sub.t=[zeros(t,1); x(:)], t=max(L.sub.b, L.sub.a1), and L.sub.b and L.sub.a are where the length of filter coefficients a and b, respectively; x(:) represents all the input samples.

[0188] For each sample x(n), the filter (e.g. IIR filter) is applied using {tilde over (b)} and as coefficients to obtain an intermediate result. The intermediate result may be obtained based on Equation (1):

y(n)={tilde over (b)}.sub.0*x(n)+{tilde over (b)}.sub.1*x(n1)+ . . . +{tilde over (b)}.sub.nb*x(nnb).sub.1*y(n1) . . . .sub.na*y(nna)(1)

where nb=L.sub.b1, na=L.sub.a1

[0189] In the present embodiment, the intermediate result has a first form y(n), described above. The intermediate result also has a second form obtained by retaining the sign of y(n) and casting y(n) to an unsigned integer format. Thus, the sign S of y(n) is kept and y(n) is cast to an unsigned integer format, Y.sub.f=cast(y(n)*s,u int32). In this case, y(n) is converted to positive when it is negative. It is the second form of the intermediate result that is subsequently downscaled into the first precisione.g. the precision of the input data.

[0190] The second form may further involve bit-shifting Y.sub.fe.g. left-shifting to scale Y.sub.f by 2.sup.n.sup.r.sup.n.sup.p-bit to obtain sY.sub.f (406).

[0191] The intermediate results sY.sub.f are then downscaled to the first precisioni.e. the same precision as the precision of the input data (n.sup.r, e.g., int16). The most significant bits of y(n) are then kept for the calculation of next samplethis may be achieved by right shifting out the extra bits from the intermediate results (e.g. 32 bits) to downscale back to the input data precision (e.g. 16 bits), with the least significant bits thus being discarded.

[0192] These steps are reiterated until all the samples are processed.

[0193] An illustration of the bit shift during filtering process is shown in FIG. 5.

[0194] An exemplary embodiment of the IIR filter employed in accordance with FIGS. 4 and 5 is the Chebyshev Type II band-pass filter of order 2. The parameters of one example are set in Table I. Considering the fact that the spikes are high frequency fluctuations, these frequencies are selected to eliminate interference from low frequency and power line noise. The floating coefficients and its fixed-point counterparts (converted based on data precision of 14 bits, sampling rate of data is 12,987 Hz) are presented in Table II. The Chebyshev Type II filters have flat passbands with wide transition bands, in contrast to that of the Chebyshev Type I and elliptic filters with fast roll off and ripples in passband.

TABLE-US-00001 TABLE I Parameter settings for Cheby2 band-pass filter Item Parameter Values 1 Low stop 80 Hz 2 Low cut-off 300 Hz 3 High cut-off 5000 Hz 4 High stop 6000 Hz 5 Stopband ripple 30 dB 6 Passband loss 6 dB

TABLE-US-00002 TABLE II Coefficients for Cheby2 band-pass filter with the parameters defined in Table I. Filter Filter type coefficients Coefficients values floating B 0.485235041197720 0.013531256974626 0.941960738398651 0.013531256974625 0.485235041197720 A 1.000000000000000 0.964164865350100 0.355831162559359 0.108373032453001 0.257375778830363 Fixed B 7950 222 15433 222 7950 point A 16384 15797 5830 1776 4217

Embodiment of Step 106Recursive Value-Count Tracking-Based Adaptive Threshold Calculation for Fixed-Point Implementation (RVCT-aTrFx)

[0195] A recursive value-count tracking-based fixed-point thresholding (RVCT-aTrFx) method was designed to accommodate the narrow range of the filtered integer data (i.e. values of the filtered data samples). By recursively finding the dynamic range of the data, a fixed number of groups (i.e. bins) can be defined and used to store the value-count of filtered data (the number of filtered data samples having each value in the range of possible values).

[0196] A threshold can subsequently be calculated based on the defined criteria by considering distributions of the value-count profile. An exemplary embodiment of such criteria is the median, which is the middle value that separates the higher part from the lower part of the samples. The bin index for such a median value (e.g., 50.sup.th percentilei.e. there are up to 50% of samples with higher values, and up to 50% of samples with lower values) can be found and its sample value taken as the initial threshold.

[0197] The final threshold may then calculated by properly scaling the initial threshold to approximate the standard deviation (e.g. a single deviation, though multiple standard deviations may be used in some embodiments) of the background noise. This process is recursively applied to all data segments. An adaptive threshold mode is proposed below, where the threshold is calculated adaptively at a predefined time interval with proper buffer management. For example, the threshold can be calculated based on the latest k seconds of data which are refreshed every n seconds, where k and n can be determined based on the buffer condition and processing speed requirements.

[0198] A block diagram 600 illustrating an exemplary embodiment of the proposed RVCT-aTrFx method is shown in FIG. 6 with the detail steps described as follows, in accordance with method 100 of FIG. 1. The block diagram 600 comprises an upper part 602 and a lower part 604. The threshold calculation in the lower part 604 can be done for all the adaptively buffered time segments for each predefined time interval (e.g. k seconds) discussed above.

[0199] Firstly, filtered data is received per step 104, from the filtered described with reference to FIGS. 3 to 5. Diagram 600 then performs step 106 of method 100, namely calculation of a threshold or thresholds.

[0200] Block 606 involves initialisation of one or more parameters. This may involve setting initial values for a maximum value (iMax), a minimum value (iMin), and a number of groups (interchangeably referred to as bins), that may form an ordered set of groups of non-overlapping, ordered values or values ranges (e.g. spike values). A range of the data (iRng) may also be specified, or determined from iMax and iMin, or set to 0, when required. In the present embodiment, at block 606, the method 100 initializes the maximum (iMax), minimum (iMin) and range of the data (iRng) (e.g., 32768, 32767, 0 for input data of int16). The number of bins is also initialised based on available valuesi.e. the available values for the filtered data (e.g., nBin=256). This may be determined based on historical data (from the subject in question or from a pool of subjects) or a first predetermined period of data. An accumulating counter aCnt may also be initialised (e.g. set to 0).

[0201] The absolute value for each sample may then be taken for subsequent processing.

[0202] At blocks 608, 610 and 612 the one or more parameters are updated based on values of the filtered data samples (i.e. a value of each respective filtered data sample). In the present embodiment, iMax and iMin are updated at block 608 by comparing each of them with the value (e.g. absolute values specified above) of the each sample, until all samples are processeddecision block 610. The dynamic range is obtained by the formula iRng=iMaxiMin+1 when all samples are processedblock 612where iMin and iMax can be cast to the next lower or higher power of 2 for efficient implementation.

[0203] Thus iMax and iMin are set to be, respectively, a maximum value and a minimum value across all of the filtered data samples for the processed k seconds of datai.e. all the adaptively buffered time segments for each predefined time interval.

[0204] Blocks 614 and 616 involve grouping the filtered data samples into binse.g. an ordered set of groupsbased on a value in each filtered data sample and parameter(s) recursively updated at blocks 608, 610 and 612. Blocks 614 and 616 effectively perform recursive value-count tracking. This process throws all the values (e.g. filtered data samples) into proper bins by calculating the bin index for each sample (e.g., x(k)) by

[00002]$\mathrm{iB}=\frac{x\left(k\right)-\mathrm{iMin}}{\mathrm{iRng}}\mathrm{nBin}$

in other words, for each filtered data sample, a respective bin is identified by normalizing the value of each filtered data sample over iRng.

[0205] Once the relevant bin has been identifiedi.e. the bin associated with the relevant bin indexthe counter for the relevant bin is incremented by Cnt(iB)=Cnt(iB)+1.

[0206] Blocks 618 and 620 use the bin counters from block 616 to identify a threshold bin based on a number of filtered data samples grouped into bins distributed around the threshold bin. Presently, the median is usedi.e. the threshold bin is the bin containing the median sample valuee.g. 50% of samples have either the same sample value of lower, and 50% of the samples have either the same sample value or higher where the 50% in each case can include the median sample. To achieve this, the count (Cnt) for all bins is added to the accumulating counter aCnt, i.e., aCnt=aCnt+Cnt(iB), until aCnt>Pr*Ls, where Ls denotes the total number of samples, and Pr denotes the chosen percentage to approximate the distribution (e.g., 0.5, or 50%, for median).

[0207] For the following calculations, we assume the bin index of the bin containing the median (or relevant value as determined by Pr) is id.

[0208] Block 622 involves determining the threshold based on the threshold group or bin (id) and the one or more parameters.

[0209] The threshold is calculated by the formula

[00003]$\mathrm{Tr}=C.Math..Math.1\left(\frac{\mathrm{id}\mathrm{iRng}}{\mathrm{nBin}}+\mathrm{iMin}\right)/C.Math..Math.2$

where C1 and C2 are the constants used to scale the basic threshold obtained. C1 and C2 were obtained by optimizing the classification accuracy based on history data, normally C1>1&C1<10, C2<1 & C2>0.

[0210] The buffer containing the last k seconds of samples is then refreshed and blocks 606 to 622 re-executed until all of the samples are processed.

[0211] It is noted that the data used for detecting dynamic ranges (Dt-Dr) and calculation of the threshold (Dt-Tr) can be the same or different. The computational load can be saved if the dynamic range is determined once and applied for all subsequent datae.g. buffered samples. However, this may compromise some precision in the determined threshold.

Embodiment of Step 108Crossing-Transition Triggered Spike Detection (CTT-SPD)

[0212] The basic concept behind spike detection is finding a local valley in the sample data. In the present embodiment, to produce a computationally-efficient method, detection will only be triggered if a transition in threshold crossing is found, i.e., if the difference in values between one the processing sample and its previous immediate neighboring sample (PV-INSthe sample immediately before the current sample being processed) is a crossing (crosses the threshold calculated at step 106). The method is hence named crossing-transition triggered spike detection (CTT-SPD).

[0213] The crossings can be obtained once the threshold for each channel is obtained. In one embodiment, a minimum-sized window of two samples is considered to further expedite the processing speed. To identify whether or not a sample is a spike, the sample is compared with its PV-INS. The condition for identifying a crossing are defined as (a) if a crossing-transition is detected and (b) the value of PV-INS is larger than that of the processing sample. The sample that satisfies conditions a) and b) is identified as a spike, i.e., a local valley has been found. This would eliminate false spikes. The detailed steps for computing the number of spikes for a channel is described as follows.

[0214] To perform spike detection, the number of spikes (N.sub.sk) is initialised for the processing channel k, N.sub.sk(k)=0, and pre-processed samples.

[0215] The filtered signal (e.g., F) is obtained based on the PPFF-FwHd filter described with reference to FIGS. 4 and 5.

[0216] The zero crossings are obtained, those crossings satisfying I.sub.c=F.sub.s+T.sub.rp<0, where T.sub.rp is obtained by replicating the threshold T.sub.r for channel k to the same length as F.sub.s.

[0217] The method 100 thereafter identifies whether or not the processing sample is a spike. For the i.sup.th sample, we define C.sub.n1=I.sub.c (i)I.sub.c (i1)whether or not a crossing-transition existsand C.sub.n2=F.sub.s (i1)>F.sub.s(i)whether or not a valley exists. If both conditions hold, i.e., C.sub.n1& C.sub.n2=1, N.sub.sk(k)=N.sub.sk(k)+1 the entering point (i.e. crossing of the thresholdsee FIG. 7) of the edge of the spike is found. A spike is thus detected.

[0218] This process is reiterated until all the samples are processed, and then until all the channels are processed.

[0219] The spike data (e.g. firing rates) are outputted in accordance with step 110, for every time bin width (in this context, the term bin width may refer to the shifting time window used in the decoding process, rather than the decoding bin width used to refresh the counter used in determining the firing rate). The extracted features, e.g. the firing rate, may be used as an input to a linear decoding function/method (linear combination of the features, e.g., linear discriminant analysis) or non-linear decoding function/method (e.g., nonlinear kernel-based support vector machine, multi-layer neural network etc.) to obtain a control signal for controlling the assistive device. Recursive Kalman filtering may also be used as an effective decoder to estimate the internal state of a brain signal given the sequence of noisy observations (features). The decoding function/method may be performed, in some embodiments, by the signal processing system (e.g. system 202 of FIG. 2) or in a device in which the components of system 202 are housed (e.g. a device placed on or attached to the scalp of the subject) or may be located remotely as described with reference to FIG. 2.

[0220] The optimal parameters should be obtained based on training data. In one embodiment, the parameters are set as: C1=5, C2=0.6745, and the time bin width is set as 100 ms.

[0221] FIG. 7 illustrates the detected entering point of the spike, point A. Point B was not detected since condition C.sub.n2 was not satisfied. Thus, only the inward crossing of the threshold C, commencing a valley or trough, is identified.

Decoding Spike Data

[0222] Lastly, shown in FIG. 2 is a remote system (spike decoding system) 212. The remote system 212 includes a receiver with either a wired connection to signal processing system 202, or in wireless communication with signal processing system 202. The remote system 212 receives the spike data generated at step 110, and decodes that spike data using any appropriate method. Accordingly, the method 100 may further comprise decoding the spike data to determine a control commandthe decoding will occur in remote system 212 in the present embodiment.

Experimental Evaluation

[0223] Experiments were conducted to validate the effectiveness of the proposed signal processing, i.e., filter-thresholding-spike detection, for fixed-point DSP implementation in BMI decoding. The results were compared with that obtained using floating point in terms of filtered data, thresholds, number of spikes and accuracies in decoding.

Data Set

[0224] All procedures and experiments were performed in compliance with the Institutional Animal Care and Use Committee (IACUC). The data set employed was obtained by training a non-human primate (NHP) to drive a mobile robot to move forward, turn left or right, or stay still using a joystick for each single trial, cued by a liquid reward from the trainer. A reward was given for the successfully performed task within 15 seconds. Three floating microwire arrays with 96 channels were used to acquire the spike signals at the left primary motor cortex, with sampling rate of 12,987 Hz. The analog joystick signals were recorded, where four classes representing the right, left and forward movements, and stop actions were defined with appropriate voltage thresholds. Data were collected on 9 non-consecutive days spanning 3 months, with each session consisting of about 20 trials for each class.

Comparison of the Filtered Data and Thresholds Obtained

[0225] Firstly, a comparison was made between the filtered data and thresholds obtained based on the floating-point and the present signal processing engine based on fixed-point implementation, with the filtered data for a short time segment and thresholds shown in FIG. 8the solid line represents the filtered signal using a floating point method, and the dashed line represents the filtered signal using the proposed fixed point method. Note that the median threshold were employed in the floating-point implementation. The evaluation was based on one selected channel with prominent spikes and a time segment of 0.5 second was employed for threshold calculation.

[0226] FIG. 9 shows the comparison of threshold obtained using the floating point method (solid line) and presently proposed fixed point method (broken line).

[0227] It can be seen from FIGS. 8 and 9 that the obtained filtered data and thresholds using proposed signal processing engine for fixed-point implementation match closely with that obtained using floating point implementation. This validates the effectiveness of the proposed method.

Comparison of Classification Accuracies

[0228] Table III compares the results of the accuracy of classifying the four-class data by using a floating-point processing engine and the presently described signal processing engine for fixed-point implementationthe results are for 1 NHP across 3 months for 9 non-consecutive days. In each row of Table III, the accuracy was the mean and standard error of all the accuracies of different sessions for that particular day. The last row is the mean accuracies across all the sessions over all days.

TABLE-US-00003 TABLE III Accuracies (%) obtained using floating-point and proposed digital signal processing engine for fixed-point implementation. Avg. Acc across sessions of a day (mean standard error) Experiment day floating-point Proposed fixed-point 1 87.83 1.00 77.08 2.66 2 91.90 2.27 90.31 3.12 3 87.62 1.94 87.05 0.73 4 94.36 1.55 84.40 6.53 5 91.51 1.70 80.68 5.89 6 82.12 3.63 70.86 6.91 7 78.44 2.19 87.07 1.58 8 87.31 2.91 87.57 0.90 9 80.28 5.39 81.08 3.62 Overall Avg. 86.82 1.83 82.90 2.04 Acc.

[0229] It can be observed from Table III that the proposed signal processing engine based on fixed-point implementation (presently for a signal processing system in a DSP chip implementation) yielded comparable classification accuracies compared with that of a floating-point implementation. Statistically paired t-test at 5% significance demonstrated no significant differences between the accuracies across all of the days based on the floating-point and proposed signal processing engine (p-value: 0.1343). These results revealed that the accuracies obtained for two settings were not significantly different from each other, which validated the effectiveness of the presently proposed method.

[0230] Table III shows: [0231] Data. An NHP to drive a mobile robot to move forward, turn left or right, or stay still using a joystick. [0232] Spike signals acquired at the left primary motor cortex, with sampling rate of 12,987 Hz. Data collected on 9 non-consecutive days spanning 3 months, each session consisted of about 20 trials for each class.

[0233] Exemplary methods that use neural activity to decode movement planning and execution for control the prosthetic devices will generally focus on the robustness of the decoding algorithms across days and sessions, as well as the decoding accuracy.

[0234] A combination of plan and peri-movement signals may be employed to improve the reconstruction of end-point movements. This is especially useful in the situation where pre-selecting the type of neural activity was not possible. Adaptive point-process filtering and maximum likelihood filtering may be employed in an estimator for plan (movement goal variable) and peri-movement signals (movement execution variable). The system may also be a closed loop systeme.g. a closed loop BMI systemthat decodes the neural signal from the brain, in a known manner, to control an external actuator of system 200 to provide sensory feedback. However, this will usually be performed in remote system 212 since providing decoding functions in a signal processing system as described herein would result in higher power consumptionthis is particularly important for on-chip implementations.

[0235] Mental states may be coded based on, for example, reinforcement learning, Kalman filtering, and the transition between different brain states. Functional mapping between the mental states and behaviour in view of the feedback may be learned via reinforcement learning i.e., to learn both the association between the states and control actions by the machine/system 212 and the relationship between brain activity and BMI control by the user. The system 200 may employ a neural network comprising a gamma structure front-end with time-varying gamma weight kernels, and a greedy policy backend. Rewards may also be used to update functional mapping. In general, the user will not be aware of which control actions to take, but must discover which control action yields the greatest reward.

[0236] In the decoder in system 200 (particularly remote system 212), Kalman filtering was explored for continuous decoding of the neural signal (e.g. spike data, such as firing rate and/or time of sample, outputted at step 112), in which control of a prosthetic device was achieved using intention estimating kinematics of original and modified BMIs. Velocities were modelled using the Kalman filter, where the velocity and position kinematics were encoded as intentions and feedback, respectively. The modified BMI was performed by changing the direction of velocity vectors towards an end target of the executed movements of the prosthetic device.

[0237] In order to improve the speed and robustness, the neural signals associated with direction and speed of movement were separated by the decoder and different weights were assigned. The first method is to fit the BMI (being a system incorporating the signal processing system, or the method employed by that system, depending on context) against a subject's (experimentally, the subject was a monkey) intended cursor kinematics, which is achieved by rotating the original cursor movement velocity vectors towards the target. In the second method, the velocity is modelled as the user intention and position as the feedback. This second method allows the subject to control velocity of the cursor with a measured neural signal while taking into account the influence of the current cursor position.

[0238] Parallel control of a continuous and a discrete action state decoders may be used to enable decoding of both task-related and discrete action states. A continuous decoder may be based on kinematics of the prosthetic device represented by neural signals from the motor cortex as processed by signal processing system 202. A discrete decoder is used to learn distributions and transitions among different states related to task, velocity of prosthetic device and engagement of user or subject. Internal dynamics of the neural activity was modelled by the neural dynamical structure and was incorporated into the control to restore the motor function. The estimated neural state was a linear combination of the dynamical evolution of previous neural state estimate, and newly observed spike counts in the spike data.

[0239] In order to increase the performance and robustness of BMI, two decoders were proposed, where the first decoder was for intention execution, and the second decoder was for error detection in a close-loop error fashion. The rationale lies in the assumption that neural correlates of BMI-based behavior errors exist in some brain regions, e.g., the motor cortex such as the dorsal premotor cortex or primary motor cortex. The error detection could detect, prevent and correct the errors and free the user or subject from manual correction.

[0240] The reference to any prior art in this specification is not, and should not be taken as, an acknowledgment or any form of suggestion that the prior art forms part of the common general knowledge.

[0241] In this specification and the statements made herein, unless stated otherwise, the word comprise and its variations, such as comprises and comprising, imply the inclusion of a stated integer, step, or group of integers or steps, but not the exclusion of any other integer or step or group of integers or steps.

[0242] References in this specification to any prior publication, information derived from any said prior publication, or any known matter are not and should not be taken as an acknowledgement, admission or suggestion that said prior publication, or any information derived from this prior publication or known matter forms part of the common general knowledge in the field of endeavour to which the specification relates.