m8ta
You are not authenticated, login. 

{99}  
PMID15928412[0] Naive coadaptive Control May 2005. see notes ____References____  
{1155}  
A frequent task in the lab is to sort spikes (extracellular neural action potentials) from background noise. In the lab we are working on doing this wirelessly; to minimize power consumption, spike sorting is done before the radio. In this way only times of spikes need be transmitted, saving bandwidth and power. (This necessitates a bidirectional radio protocol, but this is a worthy sacrifice). In most sorting programs (e.g. Plexon), the raw signal is first thresholded, then waveform snippets (typically 32 samples long) are compared to a template to accept/reject them, or to sort them into different units. The comparison metric is usually the meansquared error, MSE, aka the L2 norm. This makes sense, as the spike shapes are assumed to be stereotyped (they may very well not be), and the noise white / uncorrelated (another debatable assumption). On the headstage we are working with for wireless neural recording, jumps and memory moves are expensive operations, hence we've elected to do no waveform extraction, and instead match continuously match. By using the builtin MPEG compression opcodes, we can compute the L1 norm at a rate of 4 samples / clock  very efficient. However, this was more motivated by hardware considerations an not actual spike sorting practice. Literature suggests that for isolating a fixedpattern signal embedded in noise, the best solution is instead a matched filter. Hence, a careful study of spikesorting was attempted in matlab, given the following assumptions: fixed spike shape (this was extracted from real data), and uncorrelated bandlimited noise. The later was just white noise passed through a bandpass filter, e.g. cheby1(3, 2, [500/15e3 7.5/15]) Where the passband edges are 500 Hz and 15kHz, at a sampling rate of 30kHz. (Actual rate is 31.25kHz). Since the spike times are known, we can rigorously compare the Receiver Operating Characteristic (ROC) and the area under curve (AUC) for different sorting algorithms. Four were tried: L1 (as mentioned above, motivated by the MPEG opcodes), L2 (Plexon), FIR matched filter, and IIR matched filter. The latter was very much an experiment  IIR filters are efficiently implemented on the blackfin processor, and they generally require fewer taps than their equivalent FIR implementation. To find an IIR equivalent to a given FIR matched filter (whose impulse response closely looks like the actual waveshape, just timereversed), the filter parameters were simply optimized to match the two impulse responses. To facilitate the search, the denominator was specified in terms of complex conjugate pole locations (thereby constraining the form of the filter), while the numerator coefficients were individually optimized. Note that this is not optimizing given the objective to maximize sorting quality  rather, it is to make the IIR filter impulse response as close as possible to the FIR matched filter, hence computationally light. And yet: the IIR filter outperforms the FIR matched filter, even though the IIR filter has 1/3 the coefficients (10 vs 32)! Below is the AUC quality metric for the four methods. And here are representative ROC curves at varying spike SNR ratios. The remarkable thing is that even at very low SNR, the matched IIR filter can reliably sort cells from noise. (Note that the acceptable false positive here should be weighted more highly; in the present analysis true positive and false positive are weighted equally, which is decidedly nonBayesian given most of the time there is no spike.) The matched IIR filter is far superior to the normal MSE to template / L2 norm method  seems we've been doing it wrong all along? As for reliably finding spikes / templates / filters when the SNR < 0, the tests above  which assume an equal number of spike samples and nonspike samples  are highly biased; spikes are not normally sortable when the SNR < 0. Upon looking at the code again, I realized three important things:
Including #1 above, as expected, dramatically increased the false positive rate, which is to be expected and how the filters will be used in the real world. #2 did not dramatically impact any of the discriminators, which is good. #3 alleviated the gap between the IIR and FIR filters, and indeed the FIR matched filter performance now slightly exceeds the IIR matched filer. Below, AUC metric for 4 methods. And corresponding ROC for 6 different SNR ratios (note the SNRs sampled are slightly different, due to the higher false positive rate). One thing to note: as implemented, the IIR filter requires careful matching of poles and zeros, and is may not work with 1.15 fixedpoint math on the Blackfin. The method really deserves to be tested in vivo, which I shall do shortly. More updates: See www.aicit.org/jcit/ppl/JCIT0509_05.pdf  they add an 'adjustment' function to the matched filter due to variance in the amplitude of spikes, which adds a little performance at low SNRs. $F(t) = \left[ \frac{x(t)}{k \sigma} \dot e^{1\frac{x(t)}{k \sigma}} \right]^n$ Sigma is the standard deviation of x(t), n and k determine 'zoom intensity and zoom center'. The paper is not particularly well written  there are some typos, and their idea seems unjustified. Still the references are interesting:
Update: It is not to difficult to convert FIR filters to IIR filters using simple numerical optimization. Within my client program, this is done using simulated annealing; have tested this using fminsearch in matlab. To investigate the IIRfilter fitting problem more fully, I sliced the 10dimensional optimization space along pairs of dimensions about the optimum point as found using fminsearch. The parameters are as follows:
The figure below plots the +1 beyond the optimum for each axis pair. Click for full resolution image. Note that the last parameter is discrete, hence steps in the objective function. Also note that the problem is perfectly quadratic for the numerator, as expected, which is why LMS works so well. Note that for the denominator pole locations, the volume of the optimum is small, and there are interesting features beyond this. Some spaces have multiple optima. The next figure plots +0.1 beyond the optimum for each axis vs. every other one. It shows that, at least on a small scale, the problem becomes very quadratic in all axes hence amenable to line or conjugate gradient search. Moving away from planes that pass through a found optima, what does the space look like? E.g. From a naive start, how hard is it to find at least one workable solution? To test this, I perturbed the found optimum with white noise in the parameters std 0.2, and plotted the objective function as before, albeit at higher resolution (600 x 600 points for each slice). These figures show that there can be several optima in the denominator, but again it appears that a very rough exploration followed by gradient descent should arrive at an optima.  
{1125} 
ref: 0
tags: active filter design Netherlands Gerrit Groenewold
date: 02172012 20:27 gmt
revision:0
[head]


IEEE04268406 (pdf) Noise and Group Delay in Actvie Filters
 
{989}  
PMID19603074[0] Unscented Kalman filter for brainmachine interfaces.
____References____
 
{175}  
{97}  
From Uncertain Spikes to Prosthetic Control a powerpoint presentation w/ good overview of all that the Brown group has done  
{585}  
LMSbased adaptive decorrelator, xn is the noise, xs is the signal, len is the length of the signal, delay is the delay beyond which the autocorrelation function of the signal is zero but the acf of the noise is nonzero. The filter is very simple, and should be easy to implement in a DSP. function [y,e,h] = lms_test(xn, xs, len, delay) h = zeros(len, 1); x = xn + xs; for k = 1:length(x)lendelay y(k) = x(k+delay:k+len1+delay) * h ; e(k) = x(k)  y(k); h = h + 0.0004 * e(k) * x(k+delay:k+len1+delay)'; endIt works well if the noise source is predictable & stable: (black = sinusoidal noise, red = output, green = error in output) Now, what if the amplitude of the corrupting sinusoid changes (e.g. due to varying electrode properties during movement), and the changes per cycle are larger than the amplitude of the signal? The signal will be swamped! The solution to this is to adapt the decorrelating filter slowly, by adding an extra (multiplicative, nonlinear) gain term to track the error in terms of the absolute values of the signals (another nonlinearity). So, if the input signal is on average larger than the output, the gain goes up and viceversa. See the code. function [y,e,h,g] = lms_test(xn, xs, len, delay) h = zeros(len, 1); x = xn + xs; gain = 1; e = zeros(size(x)); e2 = zeros(size(x)); for k = 1:length(x)lendelay y(k) = x(k+delay:k+len1+delay) * h; e(k) = (x(k)  y(k)); h = h + 0.0002 * e(k) * x(k+delay:k+len1+delay)'; % slow adaptation. y2(k) = y(k) * gain; e2(k) = abs(x(k))  abs(y2(k)); gain = gain + 1 * e2(k) ; gain = abs(gain); if (gain > 3) gain = 3; end g(k) = gain; end If, like me, you are interested in only the abstract features of the signal, and not an accurate reconstruction of the waveform, then the gain signal (g above) reflects the signal in question (once the predictive filter has adapted). In my experiments with a length 16 filter delayed 16 samples, extracting the gain signal and filtering out outofband information yielded about +45db improvement in SNR. This was with a signal 1/100th the size of the disturbing amplitudemodulated noise. This is about twice as good as the human ear/auditory system in my tests.
It doesn't look like much, but it is just perfect for EMG signals corrupted by timevarying 60hz noise.  
{983}  
PMID21976021[0] Active tactile exploration using a brainmachinebrain interface.
____References____
 
{957}  
PMID19923243 Complex Spatiotemporal Tuning in Human UpperLimb Muscles
 
{727}  
quote: For a given amount of 10kHz rejection, a 1ksps $\Sigma \Delta$ ADC (a) is equivalent to a 1Msps SAR ADC (b). Filtering determines the ADC's effective throughput, not the sampling rate (c). exactly! (click to see full res version. from an excellent article, "Managing highvoltage lithiumion batteries in HEVs" by Michael Kultgen. more notes:
other stuff in that issue:
 
{438}  
I needed to evaluate different fixedpoint filters & implement them in a way which would moreorless easily transfer to handcoded blackfin assembly. Hence, I wrote a small program in C to test two likely alternatives: a 5th order elliptic filter, with 82db stopband rejection, or a 4th order elliptic, with 72 db stopband rejection. see {421} also, and this useful reference. (also at {584}) UPDATE: this filter is not stable/ suitable for the 1.15 signed fraction format of the blackfin processor. As in the reference, you must use a noncanonic form which puts the numerator before the denominator. I had a great deal of difficulty trying to determine why the output of the Direct Form II filter was giving crap results with poles/zeros designed for lowpass operation. The reason appears to be that the Direct FormII filter requires dynamic range > +1 for the w's (first feedback section). The numerator / b coefficients for these highpass filters is usually [1 2 1]  a highpass  which returns 0 when passed saturated data. See the signal flowchart below. The program for implementing these filters as cascaded biquad/ triquads follows. The method of generating the filter coefficients is in the comments. Note that the filter coefficients for the triquad exceed and absolute magnitude of 2, hence the triquad implementation has a fixedpoint format of s.2.13 (one bit sign, 2 bits whole part, 13 bits fractional part). The biquad has s.1.14 format  one bit for the whole part. Also note that I incorporated a scale factor of x2 into the first stage, and x4 into the second stage, as the output of the ADC is only 12 bits and we can afford to expand it to the full range. The filter responses, as designed (in floating point). 5th order is on the left, 4th on the right. The filter output to white noise, in the time domain. 5th order has, in this impementation, a 'softer shoulder', which i do not think is appropriate for this application. The same output in fourier space, confirming the softer shoulder effect. this was actually kinda unespected .. I thought an extra pole/zero would help! I guess the triquad & lower res coefficient quantization is less efficient.. ? #include <stdio.h> // gcc Wall filter_test.c o filter_test //need to test fixedpoint filtering in C (where it is easy) //before converting to handcoded assembly. short w1_1[2]= {0,0}; short w1_2[3]= {0,0,0}; short w2_1[2]= {0, 0}; short w2_2[2]= {0, 0}; /* filter 1: 5th order. biquad then a triquad. [B1, A1] = ellip(5,0.2,84, 6/31.25); %84db = 14 bits. rb = roots(B1); ra = roots(A1); pa1 = poly(ra(1:2)); pa2 = poly(ra(3:5)); pb1 = poly(rb(1:2)); pb2 = poly(rb(3:5)); b1_1 = round(pb1*sqrt(B1(1))*2^15) a1_1 = round(pa1* 2^14)* 1 b1_2 = round(pb2*sqrt(B1(1))*2^15) a1_2 = round(pa2*2^13)*1 % triquad. */ //biquad1: short b1_1[3] = {1199, 87, 1199}; short a1_1[2] = {24544 14085}; //triquad1: short b1_2[4] = {1199, 2313, 2313, 1199}; short a1_2[3] = {18090, 14160, 3893}; /*filter 2: can be implemented by 2 biquads. [B3, A3] = ellip(4,0.5,72, 6/31.25); rb = roots(B3); ra = roots(A3); pa1 = poly(ra(1:2)); pa2 = poly(ra(3:4)); pb1 = poly(rb(1:2)); pb2 = poly(rb(3:4)); b2_1 = round(pb1*sqrt(B3(1))*2^15) a2_1 = round(pa1* 2^14)* 1 b2_2 = round(pb2*sqrt(B3(1))*2^16) %total gain = 8. a2_2 = round(pa2*2^14)*1 */ //biquad2: short b2_1[3] = {2082, 914, 2082}; short a2_1[2] = {24338, 13537}; //biquad2 (2): short b2_2[3] = {4163, 6639, 4163}; short a2_2[2] = {24260, 9677}; int madd_sat(int acc, short a, short b){ //emulate the blackfin //signedinteger mode of operating. page 567 in the programming reference. int c; long long l, lo, hi; lo = (long long)(2147483647); hi = (long long)(2147483647); c = (int)a * (int)b; l = (long long)acc + (long long)c; if(l < lo) l = lo; if(l > hi) l = hi; acc = (int)l; return acc; } //need to deal with samples one at a time, as they come in, from different channels //need to quantize the coeficients //need to utilize a biquad topology / filter structure. // (this is simple  just segregate the poles. order 5 filter req. biquad & triquad. short filter_biquad(short in, short* a1, short* b1, short* w1){ int acc; short out, w; //in varies from 0 to 0x0fff  ADC is straight binary. (0 to 4095). unsigned. acc = madd_sat(0, in, 16384); acc = madd_sat(acc, w1[0], a1[0]); acc = madd_sat(acc, w1[1], a1[1]); w = (short)(acc >> 14); //hopefully this does signextended shift. acc = madd_sat(0, w, b1[0]); acc = madd_sat(acc, w1[0], b1[1]); acc = madd_sat(acc, w1[1], b1[0]); //symmetry. out = (short)(acc >> 14); w1[1] = w1[0]; w1[0] = w; return out; } short filter_triquad(short in, short* a2, short* b2, short* w2){ int acc; short out, w; //this one is a bit different, as the coefficients are > 2 in the denom. // hence have to normalize by a different fraction. //multiply the numerator by 4 for a net gain of 8. acc = madd_sat(0 , in, 8192); acc = madd_sat(acc, w2[0], a2[0]); acc = madd_sat(acc, w2[1], a2[1]); acc = madd_sat(acc, w2[2], a2[2]); w = acc >> 13; acc = madd_sat(0 , w, b2[0]); acc = madd_sat(acc, w2[0], b2[1]); acc = madd_sat(acc, w2[1], b2[1]); acc = madd_sat(acc, w2[2], b2[0]); //symmetry. out = (short)(acc >> 13); w2[2] = w2[1]; w2[1] = w2[0]; w2[0] = w; return out; } int main( int argv, char* argc[]){ //read the samples from stdin. int in, out; while(scanf("%d", &in)){ //compare 5th and 4th order filters directly. in = 2048; out = filter_biquad(in, a1_1, b1_1, w1_1); out = filter_triquad(out, a1_2, b1_2, w1_2); printf("%d ", out); out = filter_biquad(in, a2_1, b2_1, w2_1); out = filter_biquad(out, a2_2, b2_2, w2_2); printf("%d ", out); } return 0; } // gcc Wall filter_test.c o filter_test // cat filter_test_in.dat  ./filter_test > filter_test_out.dat Below, some matlab code to test the filtering. % make some noisy data. x = randn(2000, 1); x = x .* 1000; x = x + 2048; i = find(x >= 4096) x(i) = 4095; i = find(x < 0); x(i) = 0; x = round(x); fid = fopen('filter_test_in.dat', 'w'); for k = 1:length(x) fprintf(fid, '%d ', x(k)); end fprintf(fid, 'quit'); fclose(fid); Here is a blackfin assembly implementation of the filter, Note the filter weights (i0) and delays (i1, i2) need to be on two different cache/sram banks, or the processor will stall. Also note that , as per the second diagram above, the delays for biquad 1 output are synonymous to the delays for biquad 2 input, hence they are only represented once in memory. /* directform 1 biquad now, form II saturates 1.15 format. operate on the two samples in parallel (both in 1 32bit reg). r0 x(n)  the input from the serial bus. r1 x(n1) (yn1)  pingpong the delayed registers. r2 x(n2) (yn2)  do this so save read cycles. r3 y(n1) (xn1) r4 y(n2) (xn2) r5 b0 b1  (low high) r6 a0 a1 i0 reads the coeficients into the registers; it loops every 32 bytes (16 coef, 4 biquads) i1 reads the delays. it loops every 640 bytes = 10 delays * 4 bytes/delay * 16 stereo channels. only increments. i2 writes the delays, loops every 640 bytes. also only increments. if i1 and i2 are dereferenced in the same cycle, the processor will stall  each of the 1k SRAM memory banks has only one port. format of delays in memory: [x1(n1) , x1(n2) , x2(n1) aka y1(n1) , x2(n2) aka y1(n2) , x3(n1) aka y2(n1) , x3(n2) aka y2(n2) , x4(n1) aka y3(n1) , x4(n2) aka y3(n2) , y4(n1) , y4(n2) ] that's 10 delays, 4 bytes each. */ r5 = [i0++]  r1 = [i1++]; a0 = r0.l * r5.l , a1 = r0.h * r5.l  r6 = [i0++]  [i2++] = r0; a0 += r1.l * r5.h, a1 += r1.h * r5.h  r2 = [i1++] ; a0 += r2.l * r5.l, a1 += r2.h * r5.l  r3 = [i1++] ; a0 += r3.l * r6.l, a1 += r3.h * r6.l  r4 = [i1++] ; r0.l = (a0 += r4.l * r6.h), r0.h = (a1 += r4.h * r6.h) (s2rnd)  [i2++] = r1; r5 = [i0++]  [i2++] = r0; a0 = r0.l * r5.l, a1 += r0.h * r5.l  r6 = [i0++]  [i2++] = r3; a0 += r3.l * r5.h, a1 += r3.h * r5.h  r1 = [i1++]; a0 += r4.l * r5.l, a1 += r4.h * r5.l  r2 = [i1++]; a0 += r1.l * r6.l, a1 += r1.h * r6.l; r0.l = (a0 += r2.l * r6.h), r0.h = (a1 += r2.h * r6.h) (s2rnd); r5 = [i0++]  [i2++] = r0; a0 = r0.l * r5.l, a1 = r0.h * r5.l  r6 = [i0++]  [i2++] = r1; a0 += r1.l * r5.h, a1 += r1.h * r5.h  r3 = [i1++]; a0 += r2.l * r5.l, a1 += r2.h * r5.l  r4 = [i1++]; a0 += r3.l * r6.l, a1 += r3.h * r6.l; r0.l = (a0 += r4.l * r6.h), r0.h = (a1 += r4.h * r6.h) (s2rnd); r5 = [i0++]  [i2++] = r0; a0 = r0.l * r5.l, a1 += r0.h * r5.l  r6 = [i0++]  [i2++] = r3; a0 += r3.l * r5.h, a1 += r3.h * r5.h  r1 = [i1++]; a0 += r4.l * r5.l, a1 += r4.h * r5.l  r2 = [i1++]; a0 += r1.l * r6.l, a1 += r1.h * r6.l; r0.l = (a0 += r2.l * r6.h), r0.h = (a1 += r2.h * r6.h) (s2rnd); [i2++] = r0; //save the delays. [i2++] = r1; //normally this would be pipelined.  
{584}  
images/584_1.pdf  this is a highly useful paper, but is no longer available on the author's website. For the good of all, i post it here. These may also have utility:
 
{523}  
 
{94} 
ref: bookmark0
tags: particle_filter unscented monte_carlo MCMC
date: 12112007 16:46 gmt
revision:2
[1] [0] [head]


 
{493}  
PMID17614134[0] Equalization filters for multiplechannel electromyogram arrays.
____References____  
{441}  
So, in order to measure how quantizing filter coeficients affects filter response, I quantized the coefficients of a 8th order bandpass filter designed with: [B1, A1] = ellip(4,0.8,70, [600/31.25e3 6/31.25]);here is a function that quantizes & unquantizes the filter coeff, then compares the frequency responses: function [Bq, Aq, Bcoef, Acoef] = filter_quantize(B, A) % quantize filter coeficients & unquantize so as to get some idea to % the *actual* fixedpoint filter performance. % assume that everything in broken into biquads. base = 10; Aroots = roots(A); Broots = roots(B); order = length(Aroots)/2; % the number of biquads. scale = B(1).^(1/order); % distribute the gain across the biquads. for o = 0:order1 Acoef_biquad(o+1, :) = poly(Aroots(o*2+1 : o*2+2)); Bcoef_biquad(o+1, :) = poly(Broots(o*2+1 : o*2+2))*scale; end Bcoef = round(Bcoef_biquad .* 2^base); Acoef = round(Acoef_biquad .* 2^base); % now, reverse the process. Bq2 = Bcoef ./ 2^base; Aq2 = Acoef ./ 2^base; for o = 0:order1 Arootsq(o*2+1: o*2+2) = roots(Aq2(o+1, :)); Brootsq(o*2+1: o*2+2) = roots(Bq2(o+1, :)); end Aq = poly(Arootsq); Bq = poly(Brootsq).*B(1); [H, W] = freqz(B, A); [Hq, Wq] = freqz(Bq, Aq); figure plot(W, db(abs(H)), 'b') hold on plot(W, db(abs(Hq)), 'r') axis([0 pi 100 0])The result: high frequency is not much affected but low frequency is strongly affected. But this is at a quatization to 10 bits  quantization to 15 bits lead to reasonably good performance. I'm not sure if this conclusively indicates / counterindicates downsampling prior to highpassing for my application, but i would say that it does, as if you downsample by 2 the highpass cutoff frequency will be 2x larger hence the filter will be less senitive to quantization errors which affect low frequencies.  
{410}  
 
{176}  
PMID15010499[0] Recursive Bayesian Decoding of Motor Cortical Signals by Particle Filtering
____References____  
{15} 
ref: bookmark0
tags: monte_carlo MCMC particle_filter probability bayes filtering biblography
date: 002007 0:0
revision:0
[head]


http://wwwsigproc.eng.cam.ac.uk/smc/papers.html  sequential monte carlo methods. (bibliography)
 
{98}  
http://hardm.ath.cx/pdf/unscentedKalmanFilter.pdf  the square root transform. contains a nice tabulation of the original algorithm, which i what I use. http://hardm.ath.cx/pdf/unscentedKalmanFilter2000.pdf  the original, with examples of state, parameter, and dual estimation http://en.wikipedia.org/wiki/Kalman_filter  wikipedia page, also has the unscented kalman filter http://www.cs.unc.edu/~welch/kalman/media/pdf/Julier1997_SPIE_KF.pdf  Julier and Ulhmann's original paper. a bit breif. http://www.cs.ubc.ca/~murphyk/Papers/Julier_Uhlmann_mar04.pdf  Julier and Ulhmann's invited paper, quite excellent.  
{37} 
ref: bookmark0
tags: Unscented sigma_pint kalman filter speech processing machine_learning SDRE control UKF
date: 002007 0:0
revision:0
[head]


 
{100}  
http://www.cs.unc.edu/~welch/media/pdf/maybeck_ch1.pdf  great explanation!! really sensible!  
{95}  
PMID15188861 Modeling and Decoding Motor Cortical Activity Using a Switching Kalman Filter  
{92}  
with the extended kalman filter, from '92: http://ftp.ccs.neu.edu/pub/people/rjw/kalmanijcnn92.ps with the unscented kalman filter : http://hardm.ath.cx/pdf/NNTrainingwithUnscentedKalmanFilter.pdf 