TEAP (Toolbox for Emotion Analysis using Physiological Signals) doc
scpopen.m
Go to the documentation of this file.
1 function [HDR]=scpopen(arg1,CHAN,arg4,arg5,arg6)
2 % SCPOPEN reads and writes SCP-ECG files
3 %
4 % SCPOPEN is an auxillary function to SOPEN for
5 % opening of SCP-ECG files for reading ECG waveform data
6 %
7 % Use SOPEN instead of SCPOPEN
8 %
9 % See also: fopen, SOPEN,
10 
11 
12 % $Id: scpopen.m 2899 2012-02-21 00:25:35Z schloegl $
13 % (C) 2004,2006,2007,2008 by Alois Schloegl <a.schloegl@ieee.org>
14 % This is part of the BioSig-toolbox http://biosig.sf.net/
15 %
16 % BioSig is free software: you can redistribute it and/or modify
17 % it under the terms of the GNU General Public License as published by
18 % the Free Software Foundation, either version 3 of the License, or
19 % (at your option) any later version.
20 
21 if nargin<2, CHAN=0; end;
22 
23 if isstruct(arg1)
24  HDR=arg1;
25  FILENAME=HDR.FileName;
26 elseif ischar(arg1);
27  HDR.FileName=arg1;
28  fprintf(2,'Warning SCPOPEN: the use of SCPOPEN is discouraged; please use SOPEN instead.\n');
29 end;
30 
31 VER = version;
32 
33 fid = fopen(HDR.FileName,HDR.FILE.PERMISSION,'ieee-le');
34 HDR.FILE.FID = fid;
35 if ~isempty(findstr(HDR.FILE.PERMISSION,'r')), %%%%% READ
36  tmpbytes = fread(fid,inf,'uchar');
37  tmpcrc = crc16eval(tmpbytes(3:end));
38  fseek(fid, 0, 'bof');
39  HDR.FILE.CRC = fread(fid,1,'uint16');
40  if (HDR.FILE.CRC ~= tmpcrc);
41  fprintf(HDR.FILE.stderr,'Warning: CRC check failed (%x vs %x)\n',tmpcrc,HDR.FILE.CRC);
42  end;
43 
44  HDR.FILE.Length = fread(fid,1,'uint32');
45  if HDR.FILE.Length~=HDR.FILE.size,
46  fprintf(HDR.FILE.stderr,'Warning SCPOPEN: header information contains incorrect file size %i %i \n',HDR.FILE.Length,HDR.FILE.size);
47  end;
48  HDR.data = [];
49 
50  DHT = [0,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,7,-7,8,-8,9,-9;0,1,5,3,11,7,23,15,47,31,95,63,191,127,383,255,767,511,1023]';
51  prefix = [1,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,10,10];
52  PrefixLength = prefix;
53 
54  %PREFIX = [0,4,5,12,13,28,29,60,61,124,125,252,253,508,509,1020,1021,1022,1023];
55  PREFIX = [0,4,5,12,13,28,29,60,61,124,125,252,253,508,509,1020,1021,1022,1023]'.*2.^[32-prefix]';
56  codelength = [1,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,18,26];
57  mask = [1,7,7,15,15,31,31,63,63,127,127,255,255,511,511,1023,1023,1023,1023]'.*2.^[32-prefix]';
58  %MASK = dec2bin(mask);
59  %mask = [1,7,7,15,15,31,31,63,63,127,127,255,255,511,511,1023,1023,1023,1023]';
60 
61  mask2 = [1,7,7,15,15,31,31,63,63,127,127,255,255,511,511,1023,1023,1023,1023]';
62  PREFIX2 = DHT(:,2);
63 
64  HT19999 = [prefix',codelength',ones(length(prefix),1),DHT];
65  HT = [prefix',codelength',ones(length(prefix),1),DHT];
66 
67  dd = [0:255]';
68  ACC = zeros(size(dd));
69  c = 0;
70  for k2 = 1:8,
71  ACC = ACC + (dd>127).*(2^c);
72  dd = mod(dd*2, 256);
73  c = c + 1;
74  end;
75 
76  section.CRC = fread(fid,1,'uint16');
77  section.ID = fread(fid,1,'uint16');
78  section.Length = fread(fid,1,'uint32');
79  section.Version = fread(fid,[1,2],'uint8');
80  section.tmp = fread(fid,[1,6],'uint8');
81 
82  NSections = min(11,(section.Length-16)/10);
83  for k = 1:NSections,
84  HDR.Block(k).id = k;
85  HDR.Block(k).length = 0;
86  HDR.Block(k).startpos = -1;
87  end;
88  for K = 1:NSections,
89  k = fread(fid,1,'uint16');
90  len = fread(fid,1,'uint32');
91  pos = fread(fid,1,'uint32');
92  if ((k > 0) && (k < NSections))
93  HDR.Block(k).id = k;
94  HDR.Block(k).length = len;
95  HDR.Block(k).startpos = pos-1;
96 
97 %% [HDR.Block(k).id ,length(tmpbytes), HDR.Block(k).length, HDR.Block(k).length+HDR.Block(k).startpos]
98  %% FIXME: instead of min(...,FileSize) a warning or error message should be reported
99  tmpcrc = crc16eval(tmpbytes(HDR.Block(k).startpos+3:min(HDR.Block(k).startpos+HDR.Block(k).length,HDR.FILE.size)));
100 
101  if (HDR.Block(k).length>0)
102  if (tmpcrc~=(tmpbytes(HDR.Block(k).startpos+(1:2))'*[1;256]))
103  fprintf(HDR.FILE.stderr,'Warning SCPOPEN: faulty CRC %04x in section %i\n',tmpcrc,k-1);
104  end;
105  end;
106  end;
107  end;
108 
109 %%[[HDR.Block.id];[HDR.Block.length];[HDR.Block.startpos]]'
110 
111  % default values - in case Section 6 is missing
112  HDR.NS = 0; HDR.SPR = 0; HDR.NRec = 0; HDR.Calib = zeros(1,0);
113  secList = find([HDR.Block.length]);
114  for K = secList(1:end),
115  if fseek(fid,HDR.Block(K).startpos,'bof');
116  fprintf(HDR.FILE.stderr,'Warning SCPOPEN: section %i not available, although it is listed in Section 0\n',secList(K+1));
117  end;
118  section.CRC = fread(fid,1,'uint16');
119  section.ID = fread(fid,1,'uint16');
120  section.Length = fread(fid,1,'uint32');
121  section.Version = fread(fid,[1,2],'uint8');
122  section.tmp = fread(fid,[1,6],'uint8');
123 
124  HDR.SCP.Section{find(K==secList)} = section;
125  if (section.Length==0),
126  elseif section.ID==0,
127  NSections = (section.Length-16)/10;
128  for k = 1:NSections,
129  HDR.Block(k).id = fread(fid,1,'uint16');
130  HDR.Block(k).length = fread(fid,1,'uint32');
131  HDR.Block(k).startpos = fread(fid,1,'uint32')-1;
132  end;
133 
134  elseif section.ID==1,
135  tag = 0;
136  k1 = 0;
137  Sect1Len = section.Length-16;
138  ListOfRequiredTags = [2,14,25,26];
139  ListOfRecommendedTags = [0,1,5,8,15,34];
140  while (tag~=255) && (Sect1Len>2),
141  tag = fread(fid,1,'uint8');
142  len = fread(fid,1,'uint16');
143  Sect1Len = Sect1Len - 3 - len;
144 %% [tag,len,Sect1Len], %% DEBUGGING information
145  field = fread(fid,[1,len],'uchar');
146  if tag == 0,
147  ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
148  HDR.Patient.Name = char(field); %% LastName
149  elseif tag == 1,
150  ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
151  HDR.Patient.FirstName = char(field);
152  elseif tag == 2,
153  ListOfRequiredTags(find(ListOfRequiredTags==2))=[];
154  HDR.Patient.Id = char(field);
155  elseif tag == 3,
156  HDR.Patient.LastName2 = char(field);
157  elseif tag == 4,
158  HDR.Patient.Age = field(1:2)*[1;256];
159  tmp = field(3);
160  if tmp==1, HDR.Patient.Age = HDR.Patient.Age; % unit='Y';
161  elseif tmp==2, HDR.Patient.Age = HDR.Patient.Age/12; % unit='M';
162  elseif tmp==3, HDR.Patient.Age = HDR.Patient.Age/52; % unit='W';
163  elseif tmp==4, HDR.Patient.Age = HDR.Patient.Age/365.25; % unit='d';
164  elseif tmp==5, HDR.Patient.Age = HDR.Patient.Age/(365.25*24); %unit='h';
165  else warning('units of age not specified');
166  end;
167  elseif (tag == 5)
168  ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
169  if any(field(1:4)~=0)
170  HDR.Patient.Birthday = [field(1:2)*[1;256],field(3:4),12,0,0];
171  end;
172  elseif (tag == 6)
173  if any(field(1:3)),
174  HDR.Patient.Height = field(1:2)*[1;256];
175  tmp = field(3);
176  if tmp==1, % unit='cm';
177  elseif tmp==2, HDR.Patient.Height = HDR.Patient.Height*2.54; %unit='inches';
178  elseif tmp==3, HDR.Patient.Height = HDR.Patient.Height*0.1; %unit='mm';
179  else warning('units of height not specified');
180  end;
181  end;
182  elseif (tag == 7)
183  if any(field(1:3)),
184  HDR.Patient.Weight = field(1:2)*[1;256];
185  tmp = field(3);
186  if tmp==1, % unit='kg';
187  elseif tmp==2, HDR.Patient.Weight = HDR.Patient.Weight/1000; %unit='g';
188  elseif tmp==3, HDR.Patient.Weight = HDR.Patient.Weight/2.2; %unit='pound';
189  elseif tmp==4, HDR.Patient.Weight = HDR.Patient.Weight*0.0284; %unit='ounce';
190  else warning('units of weight not specified');
191  end;
192  end;
193  elseif tag == 8,
194  ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
195  HDR.Patient.Sex = field;
196  elseif tag == 9,
197  HDR.Patient.Race = field;
198  elseif tag == 10,
199 field,
200  if (field(1)~=0)
201  HDR.Patient.Medication = field;
202  else
203  HDR.Patient.Medication.Code = {field(2:3)};
204  HDR.Patient.Medication.txt = field(4:end);
205  end;
206  elseif tag == 11,
207  HDR.Patient.BloodPressure.Systolic = field*[1;256];
208  elseif tag == 12,
209  HDR.Patient.BloodPressure.Diastolic = field*[1;256];
210  elseif tag == 13,
211  HDR.Patient.Diagnosis = char(field);
212  elseif tag == 14,
213  ListOfRequiredTags(ListOfRequiredTags==tag)=[];
214  HDR.SCP1.AcquiringDeviceID = char(field);
215  HDR.VERSION = field(15)/10;
216  elseif tag == 15,
217  ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
218  HDR.SCP1.AnalyisingDeviceID = char(field);
219  elseif tag == 16,
220  HDR.SCP1.AcquiringInstitution = char(field);
221  elseif tag == 17,
222  HDR.SCP1.AnalyzingInstitution = char(field);
223  elseif tag == 18,
224  HDR.SCP1.AcquiringDepartment = char(field);
225  elseif tag == 19,
226  HDR.SCP1.AnalyisingDepartment = char(field);
227  elseif tag == 20,
228  HDR.SCP1.Physician = char(field);
229  elseif tag == 21,
230  HDR.SCP1.LatestComfirmingPhysician = char(field);
231  elseif tag == 22,
232  HDR.SCP1.Technician = char(field);
233  elseif tag == 23,
234  HDR.SCP1.Room = char(field);
235  elseif tag == 24,
236  HDR.SCP1.Emergency = field;
237  elseif tag == 25,
238  ListOfRequiredTags(ListOfRequiredTags==tag)=[];
239  HDR.T0(1,1:3) = [field(1:2)*[1;256],field(3:4)];
240  elseif tag == 26,
241  ListOfRequiredTags(ListOfRequiredTags==tag)=[];
242  HDR.T0(1,4:6) = field(1:3);
243  elseif tag == 27,
244  HDR.Filter.HighPass = field(1:2)*[1;256]/100;
245  elseif tag == 28,
246  HDR.Filter.LowPass = field(1:2)*[1;256]/100;
247  elseif tag == 29,
248  if (field==0)
249  HDR.FILTER.Notch = NaN;
250  elseif bitand(field,1)
251  HDR.FILTER.Notch = 60; % 60Hz Notch
252  elseif bitand(field,2)
253  HDR.FILTER.Notch = 50; % 50Hz Notch
254  elseif bitand(field,3)==0
255  HDR.FILTER.Notch = -1; % Notch Off
256  end;
257  HDR.SCP1.Filter.BitMap = field;
258  elseif tag == 30,
259  HDR.SCP1.FreeText = char(field);
260  elseif tag == 31,
261  HDR.SCP1.ECGSequenceNumber = char(field);
262  elseif tag == 32,
263  HDR.SCP1.MedicalHistoryCodes = char(field);
264  elseif tag == 33,
265  HDR.SCP1.ElectrodeConfigurationCodes = field;
266  elseif tag == 34,
267  ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
268  HDR.SCP1.Timezone = field;
269  elseif tag == 35,
270  HDR.SCP1.MedicalHistory = char(field);
271  elseif tag == 255,
272  % section terminator
273  elseif tag >= 200,
274  % manufacturer specific - not standardized
275  else
276  fprintf(HDR.FILE.stderr,'Warning SCOPEN: unknown tag %i (section 1)\n',tag);
277  end;
278  end;
279  if ~isempty(ListOfRequiredTags)
280  fprintf(HDR.FILE.stderr,'Warning SCPOPEN: the following tags are required but missing in file %s\n',HDR.FileName);
281  disp(ListOfRequiredTags);
282  end;
283  if ~isempty(ListOfRecommendedTags)
284  fprintf(HDR.FILE.stderr,'Warning SCPOPEN: the following tags are recommended but missing in file %s\n',HDR.FileName);
285  disp(ListOfRecommendedTags);
286  end;
287 
288  elseif section.ID==2, % Huffman tables
289  HDR.SCP2.NHT = fread(fid,1,'uint16');
290  HDR.SCP2.NCT = fread(fid,1,'uint16');
291  if HDR.SCP2.NHT~=19999,
292  NHT = HDR.SCP2.NHT;
293  else
294  NHT = 0;
295  end;
296  k3 = 0;
297  for k1 = 1:NHT,
298  HT1 = zeros(HDR.SCP2.NCT,5);
299  for k2 = 1:HDR.SCP2.NCT,
300  tmp = fread(fid,3,'uint8') ;
301  HDR.SCP2.prefix = tmp(1); % PrefixLength
302  HDR.SCP2.codelength = tmp(2); % CodeLength
303  HDR.SCP2.TableModeSwitch = tmp(3);
304  tmp(4) = fread(fid,1,'int16'); % BaseValue
305  tmp(5) = fread(fid,1,'uint32'); % BaseCode
306  k3 = k3 + 1;
307  HT (k3,:) = [tmp'];
308  HT1(k2,:) = [tmp'];
309  end;
310  HDR.SCP2.HTree{k1} = makeTree(HT1);
311  HDR.SCP2.HTs{k1} = HT1;
312  end;
313  if HDR.SCP2.NHT~=19999,
314  HDR.SCP2.HT = HT;
315  else
316  tmp = size(HT19999,1);
317  HDR.SCP2.HT = [ones(tmp,1),[1:tmp]',HT19999];
318  HDR.SCP2.HTree{1} = makeTree(HT19999);
319  HDR.SCP2.HTs{1} = HT19999;
320  end;
321 
322  elseif section.ID==3,
323  HDR.NS = fread(fid,1,'uint8');
324  HDR.FLAG.Byte = fread(fid,1,'uint8');
325  if ~bitand(HDR.FLAG.Byte,4)
326  fprintf(HDR.FILE.stdout,'Warning SCPOPEN: not all leads simultaneously recorded - this mode is not supported.\n');
327  end;
328 
329  HDR.FLAG.ReferenceBeat = mod(HDR.FLAG.Byte,2);
330  %HDR.NS = floor(mod(HDR.FLAG.Byte,128)/8);
331  for k = 1:HDR.NS,
332  HDR.LeadPos(k,1:2) = fread(fid,[1,2],'uint32');
333  HDR.LeadIdCode(k,1) = fread(fid,1,'uint8');
334  end;
335  HDR.N = max(HDR.LeadPos(:))-min(HDR.LeadPos(:))+1;
336  HDR.AS.SPR = HDR.LeadPos(:,2)-HDR.LeadPos(:,1)+1;
337  HDR.SPR = HDR.AS.SPR(1);
338  for k = 2:HDR.NS
339  HDR.SPR = lcm(HDR.SPR,HDR.AS.SPR(k));
340  end;
341 
342  HDR = leadidcodexyz(HDR);
343  for k = 1:HDR.NS,
344  if 0,
345  elseif (HDR.LeadIdCode(k)==0),
346  HDR.Label{k} = 'unspecified lead';
347  elseif (HDR.VERSION <= 1.3) && (HDR.LeadIdCode(k) < 86),
348  % HDR.Label{k} = H.Label(H.LeadIdCode==HDR.LeadIdCode(k));
349  elseif (HDR.VERSION <= 1.3) && (HDR.LeadIdCode(k) > 99),
350  HDR.Label{k} = 'manufacturer specific';
351  elseif (HDR.VERSION >= 2.0) && (HDR.LeadIdCode(k) < 151),
352  % HDR.Label{k} = H.Label(H.LeadIdCode==HDR.LeadIdCode(k));
353  elseif (HDR.VERSION >= 2.0) && (HDR.LeadIdCode(k) > 199),
354  HDR.Label{k} = 'manufacturer specific';
355  else
356  HDR.Label{k} = 'reserved';
357  end;
358  end;
359  HDR.Label = strvcat(HDR.Label);
360 
361  elseif section.ID==4,
362  HDR.SCP4.L = fread(fid,1,'int16');
363  HDR.SCP4.fc0 = fread(fid,1,'int16');
364  HDR.SCP4.N = fread(fid,1,'int16');
365  HDR.SCP4.type = fread(fid,[7,HDR.SCP4.N],'uint16')'*[1,0,0,0; 0,1,0,0;0,2^16,0,0; 0,0,1,0;0,0,2^16,0; 0,0,0,1;0,0,0,2^16];
366 
367  tmp = fread(fid,[2*HDR.SCP4.N],'uint32');
368  HDR.SCP4.PA = reshape(tmp,2,HDR.SCP4.N)';
369  HDR.SCP4.pa = [0;tmp;HDR.N];
370 
371  elseif any(section.ID==[5,6]),
372 
373  SCP = [];
374  SCP.Cal = fread(fid,1,'int16')/1e6; % quant in nV, converted into mV
375  SCP.PhysDim = 'mV';
376  SCP.Dur = fread(fid,1,'int16');
377  SCP.SampleRate = 1e6/SCP.Dur;
378  SCP.FLAG.DIFF = fread(fid,1,'uint8');
379  SCP.FLAG.bimodal_compression = fread(fid,1,'uint8');
380 
381  if isnan(HDR.NS),
382  HDR.ERROR.status = -1;
383  HDR.ERROR.message = sprintf('Error SCPOPEN: could not read %s\n',HDR.FileName);
384  fprintf(HDR.FILE.stderr,'Error SCPOPEN: could not read %s\n',HDR.FileName);
385  return;
386  end;
387 
388  if CHAN==0, CHAN = 1:HDR.NS; end;
389  SCP.SPR = fread(fid,HDR.NS,'uint16');
390  HDR.InChanSelect = CHAN;
391 
392  if section.ID==6,
393  HDR.HeadLen = ftell(fid);
394  HDR.FLAG.DIFF = SCP.FLAG.DIFF;
395  HDR.FLAG.bimodal_compression = SCP.FLAG.bimodal_compression;
396  HDR.data = [];
397  outlen = HDR.SPR;
398  HDR.Calib = sparse(2:HDR.NS+1, 1:HDR.NS, SCP.Cal);
399  elseif isfield(HDR,'SCP4') %% HACK: do no know whether it is correct
400  outlen = floor(1000*HDR.SCP4.L/SCP.Dur);
401  else
402  outlen = inf;
403  end;
404 
405  if ~isfield(HDR,'SCP2'),
406  if any(SCP.SPR(1)~=SCP.SPR),
407  error('SCPOPEN: SPR do not fit');
408  else
409  S2 = fread(fid,[SCP.SPR(1)/2,HDR.NS],'int16');
410  end;
411  %S2 = S2(:,CHAN);
412 
413  elseif (HDR.SCP2.NHT==1) && (HDR.SCP2.NCT==1) && (HDR.SCP2.prefix==0),
414  codelength = HDR.SCP2.HT(1,4);
415  if (codelength==16)
416  S2 = fread(fid,[HDR.N,HDR.NS],'int16');
417  elseif (codelength==8)
418  S2 = fread(fid,[HDR.N,HDR.NS],'int8');
419  else
420  fprintf(HDR.FILE.stderr,'Warning SCPOPEN: codelength %i is not supported yet.',codelength);
421  fprintf(HDR.FILE.stderr,' Contact <a.schloegl@ieee.org>\n');
422  return;
423  end;
424  %S2 = S2(:,CHAN);
425 
426  elseif 1, HDR.SCP2.NHT~=19999;
427  %% User specific Huffman table
428  %% a more elegant Huffman decoder is used here %%
429  for k = 1:HDR.NS,
430  SCP.data{k} = fread(fid,SCP.SPR(k),'uint8');
431  end;
432 % S2 = repmat(NaN,outlen,length(HDR.InChanSelect));
433  clear S2;
434  sz = inf;
435  for k3 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k3); %HDR.NS,
436  outdata{k3} = DecodeHuffman(HDR.SCP2.HTree,HDR.SCP2.HTs,SCP.data{k},outlen);
437  sz = min(sz,length(outdata{k3}));
438  end;
439 
440  for k3 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k3); %HDR.NS,
441  S2(:,k) = outdata{k3}(1:sz);
442  end;
443  accu=0;
444 
445  elseif HDR.SCP2.NHT==19999,
446  HuffTab = DHT;
447  for k = 1:HDR.NS,
448  SCP.data{k} = fread(fid,SCP.SPR(k),'uint8');
449  end;
450  %for k = 1:HDR.NS,
451  for k3 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k3); %HDR.NS,
452  %for k = CHAN(:)',
453  s2 = SCP.data{k};
454  s2 = [s2; repmat(0,ceil(max(HDR.SCP2.HT(:,4))/8),1)];
455  k1 = 0;
456  l2 = 0;
457  accu = 0;
458  c = 0;
459  x = [];
460  HT = HDR.SCP2.HT(find(HDR.SCP2.HT(:,1)==1),3:7);
461  while (l2 < HDR.LeadPos(k,2)),
462  while ((c < max(HT(:,2))) && (k1<length(s2)-1));
463  k1 = k1 + 1;
464  dd = s2(k1);
465  accu = accu + ACC(dd+1)*(2^c);
466  c = c + 8;
467 
468  if 0, %for k2 = 1:8,
469  accu = accu + (dd>127)*(2^c);
470  dd = mod(dd*2,256);
471  c = c + 1;
472  end;
473  end;
474 
475  ixx = 1;
476  %acc = mod(accu,2^32); % bitand returns NaN if accu >= 2^32
477  acc = accu - 2^32*fix(accu*(2^(-32))); % bitand returns NaN if accu >= 2^32
478  while (bitand(acc,2^HT(ixx,1)-1) ~= HT(ixx,5)),
479  ixx = ixx + 1;
480  end;
481 
482  dd = HT(ixx,2) - HT(ixx,1);
483  if HT(ixx,3)==0,
484  HT = HDR.SCP2.HT(find(HDR.SCP2.HT(:,1)==HT(ixx,5)),3:7);
485  fprintf(HDR.FILE.stderr,'Warning SCPOPEN: Switching Huffman Tables is not tested yet.\n');
486  elseif (dd==0),
487  l2 = l2 + 1;
488  x(l2) = HT(ixx,4);
489  else %if (HT(ixx,3)>0),
490  l2 = l2 + 1;
491  %acc2 = fix(accu*(2^(-HT(ixx,1))));
492  %tmp = mod(fix(accu*(2^(-HT(ixx,1)))),2^dd);
493 
494  tmp = fix(accu*(2^(-HT(ixx,1)))); % bitshift(accu,-HT(ixx,1))
495  tmp = tmp - (2^dd)*fix(tmp*(2^(-dd))); % bitand(...,2^dd)
496 
497  %tmp = bitand(accu,(2^dd-1)*(2^HT(ixx,1)))*(2^-HT(ixx,1));
498  % reverse bit-pattern
499  if dd==8,
500  tmp = ACC(tmp+1);
501  else
502  tmp = dec2bin(tmp);
503  tmp = [char(repmat('0',1,dd-length(tmp))),tmp];
504  tmp = bin2dec(tmp(length(tmp):-1:1));
505  end
506  x(l2) = tmp-(tmp>=(2^(dd-1)))*(2^dd);
507  end;
508  accu = fix(accu*2^(-HT(ixx,2)));
509  c = c - HT(ixx,2);
510  end;
511  x = x(:);
512  if k3==1,
513  S2=x(:,ones(1,k));
514  elseif size(x,1)==size(S2,1),
515  S2(:,k) = x;
516  else
517  fprintf(HDR.FILE.stderr,'Error SCPOPEN: Huffman decoding failed (%i) \n',size(x,1));
518  HDR.data = S2;
519  return;
520  end;
521  end;
522 
523 
524  elseif (HDR.SCP2.NHT==19999), % alternative decoding algorithm.
525  warning('this branch is experimental - it might be broken')
526  HuffTab = DHT;
527  for k = 1:HDR.NS,
528  SCP.data{k} = fread(fid,SCP.SPR(k),'uint8');
529  end;
530  %for k = 1:HDR.NS,
531  for k3 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k3); %HDR.NS,
532  %for k = CHAN(:)',
533  tmp = SCP.data{k};
534  accu = [tmp(4)+256*tmp(3)+65536*tmp(2)+2^24*tmp(1)];
535  %accu = bitshift(accu,HDR.SCP2.prefix,32);
536  c = 0; %HDR.SCP2.prefix;
537  l = 4;
538  l2 = 0;
539  clear x;
540  Ntmp = length(tmp);
541  tmp = [tmp; zeros(4,1)];
542  while c <= 32, %1:HDR.SPR(k),
543  ixx = 1;
544  while (bitand(accu,mask(ixx)) ~= PREFIX(ixx)),
545  ixx = ixx + 1;
546  end;
547 
548  if ixx < 18,
549  c = c + prefix(ixx);
550  %accu = bitshift(accu, prefix(ixx),32);
551  accu = mod(accu.*(2^prefix(ixx)),2^32);
552  l2 = l2 + 1;
553  x(l2) = HuffTab(ixx,1);
554 
555  elseif ixx == 18,
556  c = c + prefix(ixx) + 8;
557  %accu = bitshift(accu, prefix(ixx),32);
558  accu = mod(accu.*(2^prefix(ixx)),2^32);
559  l2 = l2 + 1;
560 
561  acc1 = mod(floor(accu*2^(-24)),256);
562  %accu = bitshift(accu, 8, 32);
563  accu = mod(accu*256, 2^32);
564 
565  x(l2) = acc1-(acc1>=2^7)*2^8;
566  acc2 = 0;
567  for kk = 1:8,
568  acc2 = acc2*2 + mod(acc1,2);
569  acc1 = floor(acc1/2);
570  end;
571 
572  elseif ixx == 19,
573  c = c + prefix(ixx);
574  %accu = bitshift(accu, prefix(ixx),32);
575  accu = mod(accu.*(2^prefix(ixx)),2^32);
576  l2 = l2 + 1;
577  while (c > 7) && (l < Ntmp),
578  l = l+1;
579  c = c-8;
580  accu = accu + tmp(l)*2^c;
581  end;
582 
583  acc1 = mod(floor(accu*2^(-16)),2^16);
584  %accu = bitshift(accu, 16, 32);
585  accu = mod(accu.*(2^16), 2^32);
586 
587  x(l2) = acc1-(acc1>=2^15)*2^16;
588  acc2 = 0;
589  for kk=1:16,
590  acc2 = acc2*2+mod(acc1,2);
591  acc1 = floor(acc1/2);
592  end;
593  %x(l2) = acc2;
594  c = c + 16;
595  end;
596 
597  while (c > 7) && (l < Ntmp),
598  l = l+1;
599  c = c-8;
600  accu = accu + tmp(l)*(2^c);
601  end;
602  end;
603 
604  x = x(1:end-1)';
605  if k3==1,
606  S2=x(:,ones(1,k));
607  elseif size(x,1)==size(S2,1),
608  S2(:,k) = x;
609  elseif 1,
610  fprintf(HDR.FILE.stderr,'Error SCPOPEN: length=%i of channel %i different to length=%i of channel 1 \n',size(x,1),k,size(S2,1));
611  return;
612  else
613  fprintf(HDR.FILE.stderr,'Error SCPOPEN: Huffman decoding failed (%i) \n',size(x,1));
614  HDR.data=S2;
615  return;
616  end;
617  end;
618 
619  elseif HDR.SCP2.NHT~=19999,
620  %% OBSOLETE %%
621  fprintf(HDR.FILE.stderr,'Error SOPEN SCP-ECG: user specified Huffman Table not supported\n');
622  HDR.SCP = SCP;
623  return;
624 
625  else
626  HDR.SCP2,
627  end;
628 
629  % Decoding of Difference encoding
630  if SCP.FLAG.DIFF==2,
631  for k1 = 3:size(S2,1);
632  S2(k1,:) = S2(k1,:) + [2,-1] * S2(k1-(1:2),:);
633  end;
634  elseif SCP.FLAG.DIFF==1,
635  S2 = cumsum(S2);
636  end;
637 
638  if section.ID==5,
639  HDR.SCP5 = SCP;
640  HDR.SCP5.data = S2;
641  HDR.SampleRate = SCP.SampleRate;
642 
643  elseif section.ID==6,
644  HDR.SCP6 = SCP;
645  HDR.SampleRate = SCP.SampleRate;
646  HDR.PhysDim = repmat({HDR.SCP6.PhysDim},HDR.NS,1);
647  HDR.data = S2;
648 
649  if HDR.FLAG.bimodal_compression,
650  %% FIXME: THIS IS A HACK - DO NOT KNOW WHETHER IT IS CORRECT.
651  % HDR.FLAG.bimodal_compression = isfield(HDR,'SCP5') & isfield(HDR,'SCP4');
652  HDR.FLAG.bimodal_compression = isfield(HDR,'SCP4');
653  end;
654  if HDR.FLAG.bimodal_compression,
655  if isfield(HDR,'SCP5')
656  F = HDR.SCP5.SampleRate/HDR.SCP6.SampleRate;
657  HDR.SampleRate = HDR.SCP5.SampleRate;
658  HDR.FLAG.F = F;
659  else
660  HDR.FLAG.F = 1;
661  end;
662 
663  tmp=[HDR.SCP4.PA(:,1);HDR.LeadPos(1,2)]-[1;HDR.SCP4.PA(:,2)+1];
664  if ~all(tmp==floor(tmp))
665  tmp,
666  end;
667  t = (1:HDR.N) / HDR.SampleRate;
668  S1 = zeros(HDR.N, HDR.NS);
669 
670  p = 1;
671  k2 = 1;
672  pa = [HDR.SCP4.PA;NaN,NaN];
673  flag = 1;
674  %% FIXME: accu undefined ##
675  accu = 0;
676  for k1 = 1:HDR.N,
677  if k1 == pa(p,2)+1,
678  flag = 1;
679  p = p+1;
680  accu = S2(k2,:);
681  elseif k1 == pa(p,1),
682  flag = 0;
683  k2 = ceil(k2);
684  end;
685 
686  if flag,
687  S1(k1,:) = ((F-1)*accu + S2(fix(k2),:)) / F;
688  k2 = k2 + 1/F;
689  else
690  S1(k1,:) = S2(k2,:);
691  k2 = k2 + 1;
692  end;
693  end;
694 
695  HDR.SCP.S2 = S2;
696  HDR.SCP.S1 = S1;
697  S2 = S1;
698  end;
699 
700  if HDR.FLAG.ReferenceBeat && ~isfield(HDR,'SCP5')
701  fprintf(HDR.FILE.stderr,'Warning SOPEN SCP-ECG: Flag ReferenceBeat set, but no section 5 (containing the reference beat) is available\n');
702  elseif HDR.FLAG.ReferenceBeat,
703 
704  tmp_data = HDR.SCP5.data*(HDR.SCP5.Cal/HDR.SCP6.Cal);
705  for k = find(~HDR.SCP4.type(:,1)'),
706  t1 = (HDR.SCP4.type(k,2):HDR.SCP4.type(k,4));
707  t0 = t1 - HDR.SCP4.type(k,3) + HDR.SCP4.fc0;
708  S2(t1,:) = S2(t1,:) + tmp_data(t0,:);
709  end;
710  end;
711  HDR.data = S2;
712  end;
713 
714  elseif section.ID==7,
715  HDR.SCP7.byte1 = fread(fid,1,'uint8');
716  HDR.SCP7.Nspikes = fread(fid,1,'uint8');
717  HDR.SCP7.meanPPI = fread(fid,1,'uint16');
718  HDR.SCP7.avePPI = fread(fid,1,'uint16');
719 
720  for k=1:HDR.SCP7.byte1,
721  HDR.SCP7.RefBeat{k} = fread(fid,16,'uint8');
722  %HDR.SCP7.RefBeat1 = fread(fid,16,'uint8');
723  end;
724 
725  for k=1:HDR.SCP7.Nspikes,
726  tmp = fread(fid,16,'uint16');
727  tmp(1,2) = fread(fid,16,'int16');
728  tmp(1,3) = fread(fid,16,'uint16');
729  tmp(1,4) = fread(fid,16,'int16');
730  HDR.SCP7.ST(k,:) = tmp;
731  end;
732  for k=1:HDR.SCP7.Nspikes,
733  tmp = fread(fid,6,'uint8');
734  HDR.SCP7.ST2(k,:) = tmp;
735  end;
736  HDR.SCP7.Nqrs = fread(fid,1,'uint16');
737  HDR.SCP7.beattype = fread(fid,HDR.SCP7.Nqrs,'uint8');
738 
739  HDR.SCP7.VentricularRate = fread(fid,1,'uint16');
740  HDR.SCP7.AterialRate = fread(fid,1,'uint16');
741  HDR.SCP7.QTcorrected = fread(fid,1,'uint16');
742  HDR.SCP7.TypeHRcorr = fread(fid,1,'uint8');
743 
744  len = fread(fid,1,'uint16');
745  tag = 255*(len==0);
746  k1 = 0;
747  while tag~=255,
748  tag = fread(fid,1,'uchar');
749  len = fread(fid,1,'uint16');
750  field = fread(fid,[1,len],'uchar');
751 
752  if tag == 0,
753  HDR.Patient.LastName = char(field);
754  elseif tag == 1,
755 
756  end;
757  end;
758  HDR.SCP7.P_onset = fread(fid,1,'uint16');
759  HDR.SCP7.P_offset = fread(fid,1,'uint16');
760  HDR.SCP7.QRS_onset = fread(fid,1,'uint16');
761  HDR.SCP7.QRS_offset = fread(fid,1,'uint16');
762  HDR.SCP7.T_offset = fread(fid,1,'uint16');
763  HDR.SCP7.P_axis = fread(fid,1,'uint16');
764  HDR.SCP7.QRS_axis = fread(fid,1,'uint16');
765  HDR.SCP7.T_axis = fread(fid,1,'uint16');
766 
767  elseif section.ID==8,
768  tmp = fread(fid,9,'uint8');
769  HDR.SCP8.Report = tmp(1);
770  HDR.SCP8.Time = [[1,256]*tmp(2:3),tmp(4:8)'];
771  HDR.SCP8.N = tmp(9);
772  for k = 1:HDR.SCP8.N,
773  ix = fread(fid,1,'uint8');
774  len = fread(fid,1,'uint16');
775  tmp = fread(fid,[1,len],'uchar');
776  HDR.SCP8.Statement{k,1} = char(tmp);
777  end
778 
779  %elseif section.ID==9,
780  % HDR.SCP9.byte1 = fread(fid,1,'uint8');
781 
782  elseif section.ID==10,
783  tmp = fread(fid,2,'uint16');
784  HDR.SCP10.NumberOfLeads = tmp(1);
785  HDR.SCP10.ManufacturerCode = tmp(2);
786  for k = []; 1:HDR.SCP10.NumberOfLeads,
787  tmp = fread(fid,2,'uint16')
788  LeadId = tmp(1);
789  LeadLen = tmp(2);
790  tmp = fread(fid,LeadLen/2,'uint16');
791  HDR.SCP10.LeadId(k)=LeadId;
792  HDR.SCP10.LeadLen(k)=LeadLen;
793  HDR.SCP10.Measurements{k}=tmp;
794  end;
795 
796  elseif section.ID==11,
797  bytes = fread(fid,11,'uint8');
798  HDR.SCP11.T0 = [bytes(2)*256+bytes(3), bytes(4:8)];
799  HDR.SCP11.Confirmed = bytes(1);
800  HDR.SCP11.NumberOfStatements = bytes(9);
801  for k = 1:HDR.SCP11.NumberOfStatements,
802  SeqNo = fread(fid,1,'uint8');
803  len11 = fread(fid,1,'uint16');
804  typeID = fread(fid,1,'uint8');
805  Statement = fread(fid,len11-1,'uint8');
806  HDR.SCP11.Statement.SeqNo(k) = SeqNo;
807  HDR.SCP11.Statement.len11(k) = len11;
808  HDR.SCP11.Statement.typeID(k) = typeID;
809  HDR.SCP11.Statement.Statement{k} = Statement;
810  end;
811  end;
812 
813  if ~section.Length,
814  HDR.ERROR.status = -1;
815  HDR.ERROR.message = 'Error SCPOPEN: \n';
816  return;
817  end;
818  end;
819 
820  HDR.SPR = size(HDR.data,1);
821  HDR.NRec = 1;
822  HDR.AS.endpos = HDR.SPR;
823 
824  HDR.FILE.OPEN = 0;
825  HDR.FILE.POS = 0;
826  HDR.TYPE = 'native';
827  fclose(HDR.FILE.FID);
828 
829 else % writing SCP file
830 
831  NSections = 12;
832  SectIdHdr = zeros(1,16);
833  VERSION = round(HDR.VERSION*10);
834  if ~any(VERSION==[10,13,20])
835  fprintf(HDR.FILE.stderr,'Warning SCPOPEN(WRITE): unknown Version number %4.2f\n',HDR.VERSION);
836  VERSION = 20;
837  end;
838  SectIdHdr(9:10) = VERSION; % Section and Protocol version number
839 
840  POS = 6; B = zeros(1,POS);
841  for K = 0:NSections-1;
842  b = [];
843  if K==0,
844  % SECTION 0
845  b = [SectIdHdr(1:10),'SCPECG', zeros(1,NSections*10)];
846  b(16+(7:10)) = s4b(POS+1);
847 
848  elseif K==1,
849  % SECTION 1
850  b = SectIdHdr;
851  % tag(1),len(1:2),field(1:len)
852  if isfield(HDR.Patient,'Name'),
853  b = [b, 0, s2b(length(HDR.Patient.Name)), HDR.Patient.Name];
854  end;
855  if isfield(HDR.Patient,'Id'),
856  b = [b, 2, s2b(length(HDR.Patient.Id)), HDR.Patient.Id];
857  end;
858  if isfield(HDR.Patient,'Age'),
859  %b = [b, 4, s2b(3), s2b(HDR.Patient.Age),1]; %% use birthday instead
860  end;
861  if isfield(HDR.Patient,'Birthday'),
862  b = [b, 5, s2b(4), s2b(HDR.Patient.Birthday(1)),HDR.Patient.Birthday(2:3)];
863  end;
864  if isfield(HDR.Patient,'Height'),
865  if ~isnan(HDR.Patient.Height),
866  b = [b, 6, s2b(3), s2b(HDR.Patient.Height),1];
867  end;
868  end;
869  if isfield(HDR.Patient,'Weight'),
870  if ~isnan(HDR.Patient.Weight),
871  b = [b, 7, s2b(3), s2b(HDR.Patient.Weight),1];
872  end;
873  end;
874  if isfield(HDR.Patient,'Sex'),
875  if ~isempty(HDR.Patient.Sex)
876  sex = HDR.Patient.Sex;
877  if 0,
878  elseif isnumeric(sex), sex = sex(1);
879  elseif strncmpi(sex,'male',1); sex = 1;
880  elseif strncmpi(sex,'female',1); sex = 2;
881  else sex = 9; % unspecified
882  end;
883  b = [b, 8, s2b(1), sex];
884  end;
885  end;
886  if isfield(HDR.Patient,'Race'),
887  b = [b, 9, s2b(1), HDR.Patient.Race(1)];
888  end;
889  if isfield(HDR.Patient,'BloodPressure'),
890  b = [b, 11, s2b(2), s2b(HDR.Patient.BloodPressure.Systolic)];
891  b = [b, 12, s2b(2), s2b(HDR.Patient.BloodPressure.Diastolic)];
892  end;
893 
894  %% Tag 14
895  tag14.AnalyzingProgramRevisionNumber = ['',char(0)];
896  tag14.SerialNumberAcqDevice = ['',char(0)];
897  tag14.AcqDeviceSystemSoftware = ['',char(0)];
898  tag14.SCPImplementationSoftware = ['BioSig4OctMat v 1.76+',char(0)];
899  tag14.ManufactureAcqDevice = ['',char(0)];
900  t14 = [zeros(1,35), length(tag14.AnalyzingProgramRevisionNumber),tag14.AnalyzingProgramRevisionNumber,tag14.SerialNumberAcqDevice,tag14.AcqDeviceSystemSoftware,tag14.SCPImplementationSoftware,tag14.ManufactureAcqDevice];
901  t14(8) = 255; % Manufacturer
902  %%% ### FIXME ### t14(9:14) = % cardiograph model
903  t14(15) = VERSION; % Version
904  t14(16) = hex2dec('A0'); % Demographics and ECG rhythm data" (if we had also the reference beats we should change it in 0xC0).
905  t14(18) = hex2dec('D0'); % Capabilities of the ECG Device: 0xD0 (acquire, print and store).
906  b = [b, 14, s2b(length(t14)), t14];
907 
908  b = [b, 25, s2b(4), s2b(HDR.T0(1)),HDR.T0(2:3)];
909  b = [b, 26, s2b(3), HDR.T0(4:6)];
910  if ~any(isnan(HDR.Filter.HighPass))
911  b = [b, 27, s2b(2), s2b(round(HDR.Filter.HighPass(1)*100))];
912  end;
913  if ~any(isnan(HDR.Filter.LowPass))
914  b = [b, 28, s2b(2), s2b(round(HDR.Filter.LowPass(1)))];
915  end;
916  b = [b, 255, 0, 0]; % terminator
917  b = b + (b<0)*256;
918 
919  elseif K==3,
920  % SECTION 3
921  b = [SectIdHdr,HDR.NS,4+HDR.NS*8];
922  if ~isfield(HDR,'LeadIdCode'), HDR.LeadIdCode = zeros(1,HDR.NS); end;
923  if (numel(HDR.LeadIdCode)~=HDR.NS); warning('HDR.LeadIdCode does not have HDR.NS elements'); end;
924  if any(HDR.LeadIdCode>255), warning('invalid LeadIdCode'); end;
925  for k = 1:HDR.NS,
926  b = [b, s4b(1), s4b(HDR.SPR*HDR.NRec), mod(HDR.LeadIdCode(k),256)];
927  end;
928 
929  elseif K==6,
930  % SECTION 6
931  Cal = full(HDR.Calib(2:end,:)); Cal = Cal - diag(diag(Cal));
932  if any(Cal(:))
933  fprintf(HDR.FILE.stderr,'Calibration is not a diagonal matrix.\n\tThis can result in incorrect scalings.\n');
934  end;
935  Cal = full(diag(HDR.Calib(2:end,:)));
936  if any(Cal~=Cal(1)),
937  fprintf(HDR.FILE.stderr,'scaling information is not equal for all channels; \n\tThis is not supported by SCP and can result in incorrect scalings.\n');
938  end;
939  [tmp,scale1] = physicalunits(HDR.PhysDim{1});
940  [tmp,scale2] = physicalunits('nV');
941  b = [SectIdHdr, s2b(round(Cal(1)*scale1/scale2)), s2b(round(1e6/HDR.SampleRate)), 0, 0];
942  for k = 1:HDR.NS,
943  b = [b, s2b(HDR.SPR*HDR.NRec*2)];
944  end;
945  data = HDR.data(:);
946  data = data + (data<0)*2^16;
947  tmp = s2b(round(data))';
948  b = [b,tmp(:)'];
949  else
950  b = [];
951  end;
952  if (length(b)>0)
953  if mod(length(b),2), % align to multiple of 2-byte blocks
954  b = [b,0];
955  end;
956  if (length(b)<16), fprintf(HDR.FILE.stderr,'section header %i less then 16 bytes %i', K,length(b)); end;
957  b(3:4) = s2b(K);
958  b(5:8) = s4b(length(b));
959  %b(1:2)= s4b(crc);
960  b(1:2) = s2b(crc16eval(b(3:end)));
961  % section 0: startpos in pointer field
962  B(22+K*10+(7:10)) = s4b(POS+1);
963  end;
964  B = [B(1:POS),b];
965  % section 0 pointer field
966  B(22+K*10+(1:2)) = s2b(K);
967  B(22+K*10+(3:6)) = s4b(length(b)); % length
968  POS = POS + length(b);
969  end
970  B(3:6) = s4b(POS); % length of file
971  B(7:8) = s2b(crc16eval(B(9:22+NSections*10))); % CRC of Section 0
972 
973  B(1:2) = s2b(crc16eval(B(3:end)));
974 
975  % fwrite(fid,crc,'int16');
976  count = fwrite(fid,B,'uchar');
977  fclose(fid);
978 end
979 end %% scpopen
980 
981 
982 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
983 %
984 % Auxillary functions
985 %
986 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
987 
988 function b2 = s2b(i)
989  % converts 16bit into 2 bytes
990  b2 = [bitand(i,255),bitand(bitshift(i,-8),255)];
991  return;
992 end %%%%% s2b %%%%%
993 
994 
995 function b4 = s4b(i)
996  % converts 32 bit into 4 bytes
997  b4 = [s2b(bitand(i,2^16-1)),s2b(bitand(bitshift(i,-16),2^16-1)) ];
998  return;
999 end %%%%% s4b %%%%%
1000 
1001 
1002 function T = makeSubTree(T,bc,len,val)
1003  if (len==0)
1004  T.idxTable = val;
1005  return
1006  end;
1007 if 1,
1008  b = bitand(bc,1)+1;
1009  if ~isfield(T,'branch') T.branch = {[],[]}; end;
1010  T.branch{b} = makeSubTree(T.branch{b},bitshift(bc,-1),len-1,val);
1011 else
1012  if bitand(bc,1)
1013  if ~isfield(T,'node1'),T.node1 = []; end;
1014  T.node1 = makeSubTree(T.node1,bitshift(bc,-1),len-1,val);
1015  else
1016  if ~isfield(T,'node0'),T.node0 = []; end;
1017  T.node0 = makeSubTree(T.node0,bitshift(bc,-1),len-1,val);
1018  end;
1019 end,
1020  return;
1021 end %%%%% makeSubTree %%%%%
1022 
1023 function T = makeTree(HT)
1024  T = [];
1025  for k1 = 1:size(HT,1)
1026  for k2 = 1:HT(k1,1) % CodeLength
1027  T = makeSubTree(T,HT(k1,5),HT(k1,1),k1);
1028  end;
1029  end;
1030 %save matlab T,pause
1031  return;
1032 end %%%%% makeTree %%%%%
1033 
1034 function outdata = DecodeHuffman(HTrees,HTs,indata,outlen)
1035  ActualTable = 1;
1036  k1 = 1; r=0;
1037  k2 = 0;
1038 
1039  if ((outlen>0) && isfinite(outlen))
1040  outdata = repmat(NaN,outlen,1);
1041  else
1042  outdata = [0;0]; %% make it a column vector
1043  end;
1044  Node = HTrees{ActualTable};
1045  while ((k1*8+r <= 8*length(indata)) && (k2<outlen))
1046  if ~isfield(Node,'idxTable')
1047  r = r+1; if (r>8), k1=k1+1; r=1; end;
1048 if 1,
1049  b = bitand(bitshift(indata(k1),r-8),1)+1;
1050  if ~isempty(Node.branch{b})
1051  Node = Node.branch{b};
1052  else
1053  fprintf(2,'Warning SCPOPEN: empty node in Huffman table\n');
1054  end;
1055 else
1056  if bitand(bitshift(indata(k1),r-8),1)
1057  if isfield(Node,'node1')
1058  Node = Node.node1;
1059  else
1060  fprintf(2,'Warning SCPOPEN: empty node in Huffman table\n');
1061  end;
1062  else
1063  if isfield(Node,'node0')
1064  Node = Node.node0;
1065  else
1066  fprintf(2,'Warning SCPOPEN: empty node in Huffman table\n');
1067  end;
1068  end;
1069 end;
1070  end;
1071 
1072  if isfield(Node,'idxTable')
1073  TableEntry = HTs{ActualTable}(Node.idxTable,:);
1074  dlen = TableEntry(2)-TableEntry(1);
1075  if (~TableEntry(3))
1076  ActualTable = TableEntry(4);
1077  elseif (dlen~=0)
1078  acc = 0;
1079  for k3 = 1:dlen,
1080  r = r+1; if (r>8), k1=k1+1; r=1; end;
1081  acc = 2*acc + bitand(bitshift(indata(k1),r-8), 1);
1082  end;
1083  if (acc>=bitshift(1,dlen-1))
1084  acc = acc - bitshift(1,dlen);
1085  end;
1086  k2 = k2+1;
1087  outdata(k2) = acc;
1088  else
1089  k2 = k2+1;
1090  outdata(k2)=TableEntry(4);
1091  end;
1092  Node = HTrees{ActualTable};
1093  end;
1094  end;
1095  return;
1096 end %%%%%%%% DecodeHuffman %%%%%%%%%
1097 
1098 function crc16 = crc16eval(D)
1099 % CRC16EVAL cyclic redundancy check with the polynomiaL x^16+x^12+x^5+1
1100 % i.e. CRC-CCITT http://en.wikipedia.org/wiki/Crc16
1101 
1102  D = uint16(D);
1103 
1104  crchi = 255;
1105  crclo = 255;
1106 
1107  t = '00102030405060708191a1b1c1d1e1f112023222524272629383b3a3d3c3f3e32434041464744454a5b58595e5f5c5d53626160676665646b7a79787f7e7d7c74858687808182838c9d9e9f98999a9b95a4a7a6a1a0a3a2adbcbfbeb9b8bbbab6c7c4c5c2c3c0c1cedfdcdddadbd8d9d7e6e5e4e3e2e1e0effefdfcfbfaf9f8f9181b1a1d1c1f1e110003020504070608393a3b3c3d3e3f30212223242526272b5a59585f5e5d5c53424140474645444a7b78797e7f7c7d72636061666764656d9c9f9e99989b9a95848786818083828cbdbebfb8b9babbb4a5a6a7a0a1a2a3afdedddcdbdad9d8d7c6c5c4c3c2c1c0cefffcfdfafbf8f9f6e7e4e5e2e3e0e1e';
1108  crc16htab = hex2dec(reshape(t,2,length(t)/2)');
1109 
1110  t = '0021426384a5c6e708294a6b8cadceef31107352b594f7d639187b5abd9cffde62432001e6c7a4856a4b2809eecfac8d53721130d7f695b45b7a1938dffe9dbcc4e586a740610223cced8eaf48690a2bf5d4b79671503312fddcbf9e79583b1aa687e4c522036041ae8feccd2a0b684997b6d5f4133251709fbeddfc1b3a597888a9caeb0c2d4e6f80a1c2e304254667b998fbda3d1c7f5eb190f3d235147756eacba8896e4f2c0de2c3a08166472405dbfa99b85f7e1d3cd3f291b0577615344c6d0e2fc8e98aab44650627c0e182a37d5c3f1ef9d8bb9a75543716f1d0b3922e0f6c4daa8be8c926076445a283e0c11f3e5d7c9bbad9f81736557493b2d1f0';
1111  crc16ltab = hex2dec(reshape(t,2,length(t)/2)');
1112 
1113  for k = 1:length(D),
1114  ix = double(bitxor(crchi,D(k)))+1;
1115  crchi = bitxor(crclo,crc16htab(ix));
1116  crclo = crc16ltab(ix);
1117  end;
1118  crc16 = crchi*256+crclo;
1119 
1120 end %%%%% crc16eval %%%%%
1121 
1122 
s2b
function s2b(in i)
crc16eval
function crc16eval(in D)
makeSubTree
function makeSubTree(in T, in bc, in len, in val)
DecodeHuffman
function DecodeHuffman(in HTrees, in HTs, in indata, in outlen)
scpopen
function scpopen(in arg1, in CHAN, in arg4, in arg5, in arg6)
makeTree
function makeTree(in HT)
s4b
function s4b(in i)