1
/********************************************************************************
3
* this file is part of: *
4
* libeep, the project for reading and writing avr/cnt eeg and related files *
6
********************************************************************************
8
* LICENSE:Copyright (c) 2003-2009, *
9
* Advanced Neuro Technology (ANT) B.V., Enschede, The Netherlands *
10
* Max-Planck Institute for Human Cognitive & Brain Sciences, Leipzig, Germany *
12
********************************************************************************
14
* This library is free software; you can redistribute it and/or modify *
15
* it under the terms of the GNU Lesser General Public License as published by *
16
* the Free Software Foundation; either version 3 of the License, or *
17
* (at your option) any later version. *
19
* This library is distributed WITHOUT ANY WARRANTY; even the implied warranty *
20
* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
21
* GNU Lesser General Public License for more details. *
23
* You should have received a copy of the GNU Lesser General Public License *
24
* along with this program. If not, see <http://www.gnu.org/licenses/> *
26
*******************************************************************************/
37
char RCS_raw3_h[] = RCS_RAW3_H;
38
char RCS_raw3_c[] = "$RCSfile: raw3.c,v $ $Revision: 2415 $";
42
/* #ifdef RAW3_CHECK obsolete */
43
/* indicator flags for the citical compression methods */
44
short ERR_FLAG_16 = 0;
46
short ERR_FLAG_EPOCH = 0;
49
void raw3_reset_ERR_FLAG_EPOCH(void) { ERR_FLAG_EPOCH = 0; }
50
short raw3_ERR_FLAG_EPOCH(void) { return ERR_FLAG_EPOCH; }
52
unsigned char CheckVerbose = 0;
53
void raw3_setVerbose(int onoff)
55
assert(onoff == 0 || onoff == 1);
56
CheckVerbose = (unsigned char) onoff;
60
known residual prediction/encoding methods
61
these codes are stored in data files - never change this!
64
#define RAW3_COPY 0 /* no residuals, original 16-bit values stored */
65
#define RAW3_TIME 1 /* 16 bit residuals from first deviation */
66
#define RAW3_TIME2 2 /* 16 bit residuals from second deviation */
67
#define RAW3_CHAN 3 /* 16 bit residuals from difference of first deviation */
68
/* and first dev. of neighbor channel */
70
#define RAW3_COPY_32 8 /* the same for 32 bit sample data */
71
#define RAW3_TIME_32 9
72
#define RAW3_TIME2_32 10
73
#define RAW3_CHAN_32 11
79
find the number of bits needed to store the passed (signed!) value
87
static int nbits[128] = {
88
/* 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */
90
1, 2, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5,
91
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
92
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
93
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
95
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
96
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
97
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
98
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8
103
/* positive number which needs the same number of bits */
104
if (y < 0) y = -(y+1);
108
else if (y < 0x00008000)
109
return nbits[y >> 8] + 8;
110
else if (y < 0x00800000)
111
return nbits[y >> 16] + 16;
113
return nbits[y >> 24] + 24;
116
int huffman_size(int *bitfreq, int n, int *nbits, int *nexcbits)
118
int lmin = 1000000000;
123
/* don't store with one bit per sample point */
124
bitfreq[2] += bitfreq[1]; bitfreq[1] = 0;
126
/* determine the max bits to store */
128
while(bitfreq[j] == 0) j--;
132
calculate the length needed for coding in i bits with exceptions for all
133
the residuals out of range
135
for (i = 2; i <= *nexcbits; i++) {
137
/* all residuals out of bitrange are exceptions */
139
for (j = i + 1; j <= *nexcbits; j++) {
144
we need one value in bitrange for exception coding, which will cause
146
e.g. for 4-bit it's every 16th value; we divide it because the lower edge value
147
we are going to use for this code is less frequent than the centers
150
nexc += (n / (1 << i)) / 2;
152
lcur = (n - 1) * i + nexc * (*nexcbits);
154
lcur = (n - nexc) * i + nexc * (*nexcbits + i);
156
| exceptions (code + full value)
160
/* note the best i for return */
173
printf("%d\t%d\n", *nbits, *nexcbits);
179
/* --------------------------------------------------------------------
180
raw3 compressed data vectors
182
TIME/TIME2/CHAN:| method | nbits | nexcbits | x[0] | x[1] .. x[n-1]
183
4 4 4 16 nbits or nbits + nexcbits
185
..._32 :| method | nbits | nexcbits | x[0] | x[1] .. x[n-1]
186
4 6 6 32 nbits or nbits + nexcbits
188
COPY: | method | dummy | x[0] .. x[n-1]
191
COPY_32: | method | dummy | x[0] .. x[n-1]
196
int huffman(sraw_t *in, int n, int method, int nbits, int nexcbits,
202
sraw_t inval, loval, hival;
204
int exc, check_exc = nbits != nexcbits;
205
int nbits_1 = nbits - 1;
206
int nexcbits_1 = nexcbits - 1;
208
/*static int xxxx = 0; */
210
static unsigned char setoutbit[8] = {
211
0x01, 0x02, 0x04, 0x08,
212
0x10, 0x20, 0x40, 0x80
215
static unsigned int setinbit[32] = {
216
0x00000001, 0x00000002, 0x00000004, 0x00000008,
217
0x00000010, 0x00000020, 0x00000040, 0x00000080,
218
0x00000100, 0x00000200, 0x00000400, 0x00000800,
219
0x00001000, 0x00002000, 0x00004000, 0x00008000,
220
0x00010000, 0x00020000, 0x00040000, 0x00080000,
221
0x00100000, 0x00200000, 0x00400000, 0x00800000,
222
0x01000000, 0x02000000, 0x04000000, 0x08000000,
223
0x10000000, 0x20000000, 0x40000000, 0x80000000
226
/* first 4 bits for the method */
227
out[0] = method << 4;
229
printf("method: %d, %d, %d\n", method, out[0]>>4, out[0]);
232
/* the rest depends on the method */
233
if (method != RAW3_COPY && method != RAW3_COPY_32) {
235
/* 32 bit versions (6 bits to store number of bits, 32 bit start val.) */
237
out[0] |= nbits >> 2;
238
out[1] = (nbits << 6) | nexcbits;
239
out[2] = in[0] >> 24;
240
out[3] = in[0] >> 16;
248
/* 16 bit versions (4 bits to store number of bits, 16 bit start val.) */
251
out[1] = (nexcbits << 4) | ((in[0] >> 12) & 0x0f);
252
out[2] = (in[0] >> 4) & 0xff;
253
outval = (in[0] & 0x0f) << 4;
258
hival = 1 << nbits_1;
262
for (nin = 1; nin < n; nin++) {
266
/* out of range ? - note and force exception code */
267
if (check_exc && (inval <= loval || inval >= hival)) {
272
/* copy the value to output (bitwise) */
273
for (bitin = nbits_1; bitin >= 0; bitin--) {
274
if (inval & setinbit[bitin])
275
outval |= setoutbit[bitout];
278
out[nout++] = outval;
284
/* exception alert - write the full value */
287
for (bitin = nexcbits_1; bitin >= 0; bitin--) {
288
if (inval & setinbit[bitin])
289
outval |= setoutbit[bitout];
292
out[nout++] = outval;
300
/* don't forget to write the last bits */
302
/* printf("bitout: %d\n", bitout);*/
303
out[nout++] = outval;
307
else if (method == RAW3_COPY) {
309
for (nin = 0; nin < n; nin++) {
311
out[nout++] = (inval >> 8);
317
else if (method == RAW3_COPY_32) {
319
for (nin = 0; nin < n; nin++) {
320
out[nout++] = (in[nin] >> 24);
321
out[nout++] = (in[nin] >> 16);
322
out[nout++] = (in[nin] >> 8);
323
out[nout++] = in[nin];
328
printf("\nnout: %d\n", nout);
335
int dehuffman16(unsigned char *in, int n, int *method, sraw_t *out)
337
int nin = 0, nout = 0;
338
int nbit, nbit_1, nexcbit, nexcbit_1, check_exc;
347
static sraw_t negmask[33] = {
348
0xffffffff, 0xfffffffe, 0xfffffffc, 0xfffffff8,
349
0xfffffff0, 0xffffffe0, 0xffffffc0, 0xffffff80,
350
0xffffff00, 0xfffffe00, 0xfffffc00, 0xfffff800,
351
0xfffff000, 0xffffe000, 0xffffc000, 0xffff8000,
352
0xffff0000, 0xfffe0000, 0xfffc0000, 0xfff80000,
353
0xfff00000, 0xffe00000, 0xffc00000, 0xff800000,
354
0xff000000, 0xfe000000, 0xfc000000, 0xf8000000,
355
0xf0000000, 0xe0000000, 0xc0000000, 0x80000000,
359
static sraw_t posmask[33] = {
360
0x00000000, 0x00000001, 0x00000003, 0x00000007,
361
0x0000000f, 0x0000001f, 0x0000003f, 0x0000007f,
362
0x000000ff, 0x000001ff, 0x000003ff, 0x000007ff,
363
0x00000fff, 0x00001fff, 0x00003fff, 0x00007fff,
364
0x0000ffff, 0x0001ffff, 0x0003ffff, 0x0007ffff,
365
0x000fffff, 0x001fffff, 0x003fffff, 0x007fffff,
366
0x00ffffff, 0x01ffffff, 0x03ffffff, 0x07ffffff,
367
0x0fffffff, 0x1fffffff, 0x3fffffff, 0x7fffffff,
371
static sraw_t setbit[16] = {
372
0x0001, 0x0002, 0x0004, 0x0008,
373
0x0010, 0x0020, 0x0040, 0x0080,
374
0x0100, 0x0200, 0x0400, 0x0800,
375
0x1000, 0x2000, 0x4000, 0x8000
379
printf("dehuffman16: ");
382
*method = (in[0] >> 4) & 0x0f;
384
if (*method != RAW3_COPY) {
386
/* using 4-bit coding nbit=16 is coded as zero */
390
fprintf(stderr,"\nlibeep: critical compression method encountered "
391
"(method %d, 16 bit)\n", *method);
396
nexcbit = (in[1] >> 4) & 0x0f;
397
/* using 4-bit coding nexcbit=16 is coded as zero */
398
if (nexcbit == 0) nexcbit = 16;
400
printf("\nnbit: %d, nexcbit: %d, method: %d\n",nbit,nexcbit,*method);
402
nexcbit_1 = nexcbit - 1;
403
excval = -(1 << (nbit_1));
404
check_exc = (nbit != nexcbit);
406
out[0] = ((sraw_t) in[1] & 0x0f) << 12
407
| (sraw_t) in[2] << 4
408
| (((sraw_t) in[3] >> 4) & 0x0f);
409
if (out[0] & 0x8000) out[0] |= 0xffff0000;
414
/* we need to read 2 or 3 bytes from input to get all input value bits */
417
/* index of first byte in inbytes containing bits of current value */
418
hibytein = bitin >> 3;
420
/* max 2 inbytes are needed here for all the bits */
421
iwork = (in[hibytein] << 8) + in[hibytein + 1];
422
swork = iwork >> (((hibytein + 2) << 3) - bitin - nbit);
424
/* handle the sign in the work buffer */
425
if (swork & setbit[nbit_1]) {
426
swork |= negmask[nbit];
429
swork &= posmask[nbit];
434
/* exception ? - forget what we got and read again */
435
if (swork == excval && check_exc) {
436
/* index of first byte containing bits of current exc. value */
437
hibytein = bitin >> 3;
439
/* max 3 inbytes are needed for all the exc. bits */
440
iwork = (in[hibytein] << 16) + (in[hibytein + 1] << 8) + in[hibytein + 2];
441
swork = iwork >> (((hibytein + 3) << 3) - bitin - nexcbit);
442
/* handle the sign in the work buffer */
443
if (swork & setbit[nexcbit_1]) {
444
swork |= negmask[nexcbit];
447
swork &= posmask[nexcbit];
458
/* index of first byte in inbytes containing bits of current value */
459
hibytein = bitin >> 3;
461
/* max 3 inbytes are needed for all the bits */
462
iwork = (in[hibytein] << 16) + (in[hibytein + 1] << 8) + in[hibytein + 2];
463
swork = iwork >> (((hibytein + 3) << 3) - bitin - nbit);
465
/* handle the sign in the work buffer */
466
if (swork & setbit[nbit_1]) {
467
swork |= negmask[nbit];
470
swork &= posmask[nbit];
475
/* exception ? - forget what we got and read again */
476
if (swork == excval && nbit != nexcbit) {
477
/* index of first byte containing bits of current exc. value */
478
hibytein = bitin >> 3;
480
/* max 3 inbytes are needed for all the bits */
481
iwork = (in[hibytein] << 16) + (in[hibytein + 1] << 8) + in[hibytein + 2];
482
swork = iwork >> (((hibytein + 3) << 3) - bitin - nexcbit);
483
/* handle the sign in the work buffer */
484
if (swork & setbit[nexcbit_1]) {
485
swork |= negmask[nexcbit];
488
swork &= posmask[nexcbit];
497
if (bitin & 0x07) nin++;
502
fprintf(stderr,"\nlibeep: critical compression method encountered "
503
"(method 0 RAW3_COPY)\n");
507
for (nout = 0; nout < n; nout++) {
508
out[nout] = (((sraw_t) in[nin]) << 8) | in[nin + 1];
509
/* read s16 instead of u16 */
510
if(out[nout] > 32767) out[nout] -= 65536;
517
printf("%d | %d | %d |%d | %d\n", *method, nbit, nexcbit, nin, out[0]);
518
for(i = 0; i < n; i++) {
519
fprintf(stderr, "%d ", out[i]);
521
fprintf(stderr, "\n");
527
int dehuffman32(unsigned char *in, int n, int *method, sraw_t *out)
529
int nin = 0, nout = 0;
530
int nbit, nbit_1, nexcbit, nexcbit_1, check_exc;
532
char inval; sraw_t outval, excval;
535
static unsigned char setinbit[8] = {
536
0x01, 0x02, 0x04, 0x08,
537
0x10, 0x20, 0x40, 0x80
540
static unsigned int setoutbit[32] = {
541
0x00000001, 0x00000002, 0x00000004, 0x00000008,
542
0x00000010, 0x00000020, 0x00000040, 0x00000080,
543
0x00000100, 0x00000200, 0x00000400, 0x00000800,
544
0x00001000, 0x00002000, 0x00004000, 0x00008000,
545
0x00010000, 0x00020000, 0x00040000, 0x00080000,
546
0x00100000, 0x00200000, 0x00400000, 0x00800000,
547
0x01000000, 0x02000000, 0x04000000, 0x08000000,
548
0x10000000, 0x20000000, 0x40000000, 0x80000000
551
static unsigned int negmask[33] = {
552
0xffffffff, 0xfffffffe, 0xfffffffc, 0xfffffff8,
553
0xfffffff0, 0xffffffe0, 0xffffffc0, 0xffffff80,
554
0xffffff00, 0xfffffe00, 0xfffffc00, 0xfffff800,
555
0xfffff000, 0xffffe000, 0xffffc000, 0xffff8000,
556
0xffff0000, 0xfffe0000, 0xfffc0000, 0xfff80000,
557
0xfff00000, 0xffe00000, 0xffc00000, 0xff800000,
558
0xff000000, 0xfe000000, 0xfc000000, 0xf8000000,
559
0xf0000000, 0xe0000000, 0xc0000000, 0x80000000,
563
*method = (in[0] >> 4) & 0x0f;
565
if (*method != RAW3_COPY_32) {
566
nbit = ((in[0] << 2) & 0x3c) | ((in[1] >> 6) & 0x03);
568
nexcbit = in[1] & 0x3f;
569
/* if (nexcbit == 0) nexcbit = 64; */
570
nexcbit_1 = nexcbit - 1;
572
out[0] = ((sraw_t) in[2] << 24)
573
| ((sraw_t) in[3] << 16)
574
| ((sraw_t) in[4] << 8)
583
excval = -(1 << (nbit_1));
584
check_exc = (nbit != nexcbit);
588
/* transfer bitwise from in to out */
589
if (inval & setinbit[bitin]) {
590
outval |= setoutbit[bitout];
591
/* handle the sign in out according to sign bit in in */
592
if (bitout == nbit_1) {
593
outval |= negmask[nbit];
606
/* is a residual complete ? */
609
/* its an exception ? - forget this and read the large value which follows */
610
if (outval == excval && check_exc) {
612
for (bitout = nexcbit_1; bitout >= 0; bitout--) {
613
if (inval & setinbit[bitin]) {
614
outval |= setoutbit[bitout];
615
if (bitout == nexcbit_1) {
616
outval |= negmask[nexcbit];
628
/* now we really have the value, store it and set up for the next one */
629
out[nout++] = outval;
634
/* count last incomplete input byte */
635
if (bitin != 7) nin++;
641
for (nout = 0; nout < n; nout++) {
642
out[nout] = ((sraw_t) in[nin] << 24)
643
| ((sraw_t) in[nin+1] << 16)
644
| ((sraw_t) in[nin+2] << 8)
645
| ((sraw_t) in[nin+3] );
651
fprintf(stderr, "%d\t%d\t%ld\t%d\t%d\n", nbit, nexcbit, nin, *method, out[0]);
652
for (i = 0; i < n; i++) {
653
fprintf(stderr, "%d ", out[i]);
655
fprintf(stderr, "\n");
661
int dehuffman(unsigned char *in, int n, int *method, sraw_t *out)
663
/* method bit 3 indicates 16 or 32 bit compression */
664
if (in[0] & (unsigned char) 0x80) {
665
return dehuffman32(in, n, method, out);
668
return dehuffman16(in, n, method, out);
673
/* ---------------------------------------------------------------------
674
take native raw data vectors (neighbors)
675
calc the residuals using the supported methods
676
find the best among them and compress
677
write the compressed output to buffer
679
int compchan(raw3_t *raw3, sraw_t *last, sraw_t *cur, int n, char *out)
681
int i, imin, mi, short_method = RAW3_COPY;
685
sraw_t *res, *restime=NULL;
689
/* don't care about short vectors in each loop - force copy instead */
692
short_method = RAW3_COPY;
694
for (i = 0; i < n; i++) {
695
if (cur[i] >= 32768 || cur[i] < -32768) {
696
short_method = RAW3_COPY_32;
700
length = huffman(cur, n, short_method, 0, 0, (unsigned char *) out);
703
/* calc residuals and their bit distribution
704
using different methods
705
select the best one and compress
708
for (mi = 0; mi < RAW3_METHODC; mi++) {
710
/* prepare access to residual data table for this method */
714
memset(hst, 0, 33 * sizeof(int));
719
rc->method = RAW3_TIME;
721
for (sample = 1; sample < n; sample++)
722
hst[bitc(res[sample] = cur[sample] - cur[sample - 1])]++;
724
rc->length = huffman_size(hst, n - 1, &rc->nbits, &rc->nexcbits);
729
rc->method = RAW3_TIME2;
731
hst[bitc(res[1] = restime[1])]++;
732
for (sample = 2; sample < n; sample++)
733
hst[bitc(res[sample] = restime[sample] - restime[sample - 1])]++;
735
rc->length = huffman_size(hst, n - 1, &rc->nbits, &rc->nexcbits);
739
rc->method = RAW3_CHAN;
741
for (sample = 1; sample < n; sample++)
742
hst[bitc(res[sample] = restime[sample] - last[sample] + last[sample - 1])]++;
744
rc->length = huffman_size(hst, n - 1, &rc->nbits, &rc->nexcbits);
749
/* find the best residual method */
751
for (mi = 0; mi < RAW3_METHODC; mi++) {
753
printf("\nmethod: %d, nbits: %d, nexcbits: %d, length: %d\n",
754
raw3->rc[mi].method, raw3->rc[mi].nbits, raw3->rc[mi].nexcbits,
755
raw3->rc[mi].length);
757
if (raw3->rc[mi].length < lmin) {
759
lmin = raw3->rc[mi].length;
762
rc = &raw3->rc[imin];
764
/* need 32 bit storage ? */
765
if (rc->nexcbits > 16 || rc->res[0] < -32768 || rc->res[0] >= 32768 )
768
short_method = RAW3_COPY_32;
770
/* is the best compression better than store only ? */
771
if (rc->nbits < 32 && lmin + 3 < n * 4) {
772
length = huffman(rc->res, n, rc->method, rc->nbits, rc->nexcbits,
773
(unsigned char *) out);
777
length = huffman(cur, n, short_method, 0, 0,
778
(unsigned char *) out);
784
printf("\nmethod: %d, nbits: %d, nexcbits: %d, length: %d\n",
785
rc->method, rc->nbits, rc->nexcbits, length);
786
else printf("\nmethod: %d\n", short_method);
792
int decompchan(raw3_t *raw3, sraw_t *last, sraw_t *cur, int n, char *in)
794
int method, sample, sample_1, sample_2;
796
sraw_t *res = raw3->rc[0].res;
798
/* restore the residuals */
800
length = dehuffman((unsigned char *) in, n, &method, res);
802
/* build the values using residuals and method */
804
switch (method & 0x07) {
808
for (sample = 1; sample < n; sample++) {
809
cur[sample] = cur[sample_1] + res[sample];
816
cur[1] = cur[0] + res[1];
819
for (sample = 2; sample < n; sample++) {
821
2 * cur[sample_1] - cur[sample_2] + res[sample];
830
for (sample = 1; sample < n; sample++) {
831
cur[sample] = cur[sample_1] + last[sample] - last[sample_1] + res[sample];
837
memcpy(cur, res, n * sizeof(sraw_t));
841
#if !defined(WIN32) || defined(__CYGWIN__)
842
fprintf(stderr, "raw3: unknown compression method!\n");
850
int compepoch_mux(raw3_t *raw3, sraw_t *in, int length, char *out)
854
sraw_t *tmp, *chanbase, *cur, *last;
856
int outsize = 0, outsizealt;
860
memset(last, 0, length * sizeof(sraw_t));
862
for (chan = 0; chan < raw3->chanc; chan++) {
864
/* extract one vector from MUX buffer */
865
chanbase = &in[raw3->chanv[chan]];
867
for (sample = 0; sample < length; sample++) {
868
cur[sample] = chanbase[samplepos];
869
samplepos += raw3->chanc;
871
/* calculate compression and its statistics */
872
outsizealt = outsize;
873
outsize += compchan(raw3, last, cur, length, &out[outsize]);
874
decompchan(raw3,last,cur,length,&out[outsizealt]);
875
/* prepare for reading next channel */
876
tmp = cur; cur = last; last = tmp;
882
int decompepoch_mux(raw3_t *raw3, char *in, int length, sraw_t *out)
886
sraw_t *tmp, *chanbase, *last, *cur;
893
memset(last, 0, length * sizeof(sraw_t));
895
for (chan = 0; chan < raw3->chanc; chan++) {
898
insize += decompchan(raw3, last, cur, length, &in[insize]);
900
/* mangle into MUX buffer */
901
chanbase = &out[raw3->chanv[chan]];
903
for (sample = 0; sample < length; sample++) {
904
chanbase[samplepos] = cur[sample];
905
samplepos += raw3->chanc;
908
/* prepare for reading next channel */
909
tmp = cur; cur = last; last = tmp;
915
/* #ifdef RAW3_CHECK obsolete */
917
int decompepoch_mux_RAW3_CHECK(eeg_t *EEG, raw3_t *raw3, char *in, int length, sraw_t *out)
921
sraw_t *tmp, *chanbase, *last, *cur;
928
memset(last, 0, length * sizeof(sraw_t));
930
for (chan = 0; chan < raw3->chanc; chan++) {
933
insize += decompchan(raw3, last, cur, length, &in[insize]);
935
/* #ifdef RAW3_CHECK obsolete */
938
int n = 0; /* count chars in output line */
942
/*fflush(stderr); fflush(stdout);*/
944
if( EEG && eep_get_chanc(EEG) > 0 ) {
945
/* print label of affected channel */
946
label = eep_get_chan_label(EEG, raw3->chanv[chan]);
947
sprintf(buf,"channel %s", label );
948
fprintf(stderr, buf);
951
/* the following channels will probably be also affected */
952
/* wrap lines after 77 chars */
953
if(chan<raw3->chanc-1) {
954
fprintf(stderr," ("); n+=2;
955
for(j=chan+1; j < raw3->chanc; j++) {
956
label = eep_get_chan_label(EEG, raw3->chanv[j]);
957
if( (n+strlen(label)) < 78) {
958
fprintf(stderr," %s", label );
959
n += strlen(label) + 1;
961
fprintf(stderr,"\n%s", label);
965
fprintf(stderr," )\n");
966
} else fprintf(stderr,"\n");
972
if( EEG && eep_get_chanc(EEG) > 0 ) {
973
label = eep_get_chan_label(EEG, raw3->chanv[chan]);
974
fprintf(stderr,"channel %s\n", label );
981
/* mangle into MUX buffer */
982
chanbase = &out[raw3->chanv[chan]];
984
for (sample = 0; sample < length; sample++) {
985
chanbase[samplepos] = cur[sample];
986
samplepos += raw3->chanc;
989
/* prepare for reading next channel */
990
tmp = cur; cur = last; last = tmp;
996
void compchanv_mux(sraw_t *buf, int length,
997
short chanc, short *chanv)
999
int chan, i, j, jmax, sample;
1001
float sumi, sumii, sumj, sumjj, sumij, x, y;
1003
/* create a [chanc x chanc] square matrix */
1004
rvv = (float **) malloc(chanc * sizeof(float *));
1005
for (i = 0; i < chanc; i++)
1006
rvv[i] = (float *) malloc(chanc * sizeof(float));
1008
/* calc a full matrix of correlation coefficients */
1009
for (i = 0; i < chanc; i++) {
1010
for (j = 0; j <= i; j++) {
1016
/* we need sums, squaresums, cosums */
1017
sumi = sumii = sumj = sumjj = sumij = 0.0;
1018
for (sample = 0; sample < length; sample++) {
1019
sumi += x = buf[i + chanc * sample];
1021
sumj += y = buf[j + chanc * sample];
1025
sumi /= length; sumii /= length;
1026
sumj /= length; sumjj /= length;
1029
/* calc the correlation coeff. from sums, avoid fp exceptions here */
1030
x = (sumii - sumi * sumi) * (sumjj - sumj * sumj);
1031
y = x > 0.0 ? sqrt(x) : 0.0;
1032
rvv[i][j] = y > 1e-6 ? (sumij - sumi * sumj) / y : 0.0;
1033
rvv[j][i] = rvv[i][j];
1038
/* find a channel-similar channel-... path trough the matrix */
1040
for (chan = 1; chan < chanc; chan++) {
1041
/* mark used channels in matrix - set column to a unreachable min correlation*/
1042
for (i = 0; i < chanc; i++)
1043
rvv[i][chanv[chan - 1]] = -2.0;
1045
/* find next chan to use - most similar to last -> max. correlation in row*/
1047
i = chanv[chan - 1];
1048
for (j = 0; j < chanc; j++) {
1049
if (rvv[i][j] > rmax) {
1059
for(i = 0; i < chanc; i++)
1060
free((char *) rvv[i]);
1064
raw3_t *raw3_init(int chanc, short *chanv, int length)
1067
raw3_t *raw3 = (raw3_t *) malloc(sizeof(raw3_t));
1072
raw3->chanc = chanc;
1073
raw3->chanv = (short *) malloc(chanc * sizeof(sraw_t));
1075
for (i = 0; i < 3; i++)
1076
raw3->rc[i].res = (sraw_t *) malloc(length * sizeof(sraw_t));
1077
raw3->last = (sraw_t *) malloc(length * sizeof(sraw_t));
1078
raw3->cur = (sraw_t *) malloc(length * sizeof(sraw_t));
1080
if (!raw3->cur || !raw3->last || !raw3->chanv) {
1084
memcpy(raw3->chanv, chanv, chanc * sizeof(short));
1089
void raw3_free(raw3_t *raw3)
1094
if (raw3->chanv) free(raw3->chanv);
1095
for (i = 0; i < 3; i++)
1096
if (raw3->rc[i].res) free(raw3->rc[i].res);
1097
if (raw3->last) free(raw3->last);
1098
if (raw3->cur) free(raw3->cur);
added
removed
src/modules/data/ext/libeep/libcnt/raw3.c
