You are not authenticated, login.
text: sort by
tags: modified
type: chronology
{421} is owned by tlh24.
hide / / print
ref: -0 tags: matlab STL boost programming C++ date: 11-06-2012 21:58 gmt revision:2 [1] [0] [head]

Recently decided to move myopen's sorting program from sqlite-based persistent state to matlab persistent state, for better interoperability with lab rigs & easier introspection. For this I wrote a class for serializing STL / boost based one, two, and three dimensional resizeable containers to and from matlab files.

The code has been tested (albeit not extensively), and therefore may be of use to someone else, even if as an example. See: http://code.google.com/p/myopen/source/browse/trunk/common_host/matStor.cpp

As you can see from the header (below), the interface is nice and concise!

#ifndef __MATSTOR_H__
#define __MATSTOR_H__

class MatStor{
        typedef boost::multi_array<float, 3> array3; 
        typedef boost::multi_array<float, 2> array2; 
        std::string m_name; //name of the file.
        std::map<std::string, std::vector<float> > m_dat1; //one dimensional
        std::map<std::string, array2> m_dat2; //two dimensions
        std::map<std::string, array3> m_dat3; //three!
        MatStor(const char* fname); 
        void save(); 
        void setValue(int ch, const char* name, float val);
        void setValue2(int ch, int un, const char* name, float val);
        void setValue3(int ch, int un, const char* name, float* val, int siz);
        float getValue(int ch, const char* name, float def);
        float getValue2(int ch, int un, const char* name, float def);
        void getValue3(int ch, int un, const char* name, float* val, int siz);

hide / / print
ref: -0 tags: feedback stability resonance butterworth matlab date: 01-22-2012 03:46 gmt revision:4 [3] [2] [1] [0] [head]

Just fer kicks, I tested what happens to low-order butterworth filters when you maladjust one of the feedback coefficients.

[B, A] = butter(2, 0.1);
[h, w] = freqz(B,A);
A(2) = A(2) * 0.9;
[h2, ~] = freqz(B,A);
hold off
hold on; plot(w,abs(h2), 'r')
title('10% change in one FB filter coef 2nd order butterworth')
xlabel('freq, rads / sec'); 
ylabel('filter response');

% do the same for a higher order filter. 
[B, A] = butter(3, 0.1);
[h, w] = freqz(B,A);
A(2) = A(2) * 0.9;
[h2, ~] = freqz(B,A);
hold on
plot(w,abs(h), 'b')
plot(w,abs(h2), 'r')
title('10% change in one FB filter coef 3rd order butterworth')
xlabel('freq, rads / sec'); 
ylabel('filter response');

The filters show a resonant peak, even though feedback was reduced. Not surprising, really; a lot of systems will show reduced phase margin and will begin to oscillate when poles are moved. Does this mean that a given coefficient (anatomical area) is responsible for resonance? By itself, of course not; one can not extrapolate one effect from one manipulation in a feedback system, especially a higher-order feedback system.

This, of course hold in the mapping of digital (or analog) filters to pathology or anatomy. Pathology is likely reflective of how the loop is structured, not how one element functions (well, maybe).

For a paper, see {1083}

hide / / print
ref: -0 tags: FIR LMS decorrelation adaptive filter matlab myopen date: 01-03-2012 03:35 gmt revision:7 [6] [5] [4] [3] [2] [1] [head]

LMS-based 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 non-zero. 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)-len-delay
	y(k) = x(k+delay:k+len-1+delay) * h ; 
	e(k) = x(k) - y(k); 
	h = h + 0.0004 * e(k) * x(k+delay:k+len-1+delay)'; 
It 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 vice-versa. 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)-len-delay
	y(k) = x(k+delay:k+len-1+delay) * h; 
	e(k) = (x(k) - y(k)); 
	h = h + 0.0002 * e(k) * x(k+delay:k+len-1+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;
	g(k) = gain; 

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 out-of-band information yielded about +45db improvement in SNR. This was with a signal 1/100th the size of the disturbing amplitude-modulated 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 time-varying 60hz noise.

hide / / print
ref: Laubach-2003.03 tags: cluster matlab linux neurophysiology recording on-line data_analysis microstimulation nicolelis laubach date: 12-17-2011 00:38 gmt revision:4 [3] [2] [1] [0] [head]

IEEE-1215970 (pdf)

  • 2003
  • M. Laubach
  • Random Forests - what are these?
  • was this ever used??

follow up paper: http://spikelab.jbpierce.org/Publications/LaubachEMBS2003.pdf

  • discriminant pusuit algorithm & local regression basis (again what are these? lead me to find the lazy learning package: http://iridia.ulb.ac.be/~lazy/


Laubach, M. and Arieh, Y. and Luczak, A. and Oh, J. and Xu, Y. Bioengineering Conference, 2003 IEEE 29th Annual, Proceedings of 17 - 18 (2003.03)

hide / / print
ref: notes-0 tags: filter DSP finite precision 16bit blackfin matlab date: 07-22-2008 03:25 gmt revision:7 [6] [5] [4] [3] [2] [1] [head]

I needed to evaluate different fixed-point filters & implement them in a way which would more-or-less easily transfer to hand-coded 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 non-canonic 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 fixed-point 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 fixed-point filtering in C (where it is easy)
//before converting to hand-coded 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.
short b1_1[3] = {1199, 87, 1199}; 
short a1_1[2] = {24544 -14085};
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
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
	//signed-integer 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 sign-extended 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); 
", 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)); 
fprintf(fid, 'quit'); 

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(n-1) (yn-1)	-- ping-pong the delayed registers. 
	r2	x(n-2) (yn-2)	-- do this so save read cycles. 
	r3	y(n-1) (xn-1)
	r4	y(n-2) (xn-2)
	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(n-1) , 
	 x1(n-2) ,  
	 x2(n-1) aka y1(n-1) , 
	 x2(n-2) aka y1(n-2) , 
	 x3(n-1) aka y2(n-1) , 
	 x3(n-2) aka y2(n-2) , 
	 x4(n-1) aka y3(n-1) , 
	 x4(n-2) aka y3(n-2) , 
	 y4(n-1) , 
	 y4(n-2) ] 
	 --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.

hide / / print
ref: notes-0 tags: DSP filter quantize lowpass elliptic matlab date: 09-02-2007 15:20 gmt revision:0 [head]

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 & un-quantizes the filter coeff, then compares the frequency responses:
function [Bq, Aq, Bcoef, Acoef] = filter_quantize(B, A) 
% quantize filter coeficients & un-quantize so as to get some idea to
% the *actual* fixed-point 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:order-1
	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; 
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:order-1
	Arootsq(o*2+1: o*2+2) = roots(Aq2(o+1, :)); 
	Brootsq(o*2+1: o*2+2) = roots(Bq2(o+1, :)); 
Aq = poly(Arootsq); 
Bq = poly(Brootsq).*B(1); 
[H, W] = freqz(B, A); 
[Hq, Wq] = freqz(Bq, Aq); 
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.

hide / / print
ref: bookmark-0 tags: optimization function search matlab linear nonlinear programming date: 08-09-2007 02:21 gmt revision:0 [head]


very nice collection of links!!

hide / / print
ref: bookmark-0 tags: neural_networks machine_learning matlab toolbox supervised_learning PCA perceptron SOM EM date: 0-0-2006 0:0 revision:0 [head]

http://www.ncrg.aston.ac.uk/netlab/index.php n.b. kinda old. (or does that just mean well established?)

hide / / print
ref: notes-0 tags: matlab mysql hack date: 0-0-2006 0:0 revision:0 [head]

LD_PRELOAD=/lib/libgcc_s.so.1 matlab this allows you to command-line call mysql or other programs that use linux's standard libgcc.

hide / / print
ref: bookmark-0 tags: information entropy bit rate matlab code date: 0-0-2006 0:0 revision:0 [head]


  • concise, well documented, useful.
  • number of bins = length of vector ^ (1/3).
  • information = sum(log (bincounts / prior) * bincounts) -- this is just the divergence, same as I do it.