📄 scpopen.m
字号:
tmp(5) = fread(fid,1,'uint32');
k3 = k3 + 1;
HT(k3,:) = [k1,k2,tmp'];
end;
end;
if HDR.SCP2.NHT~=19999,
HDR.SCP2.HT = HT;
else
tmp = size(HT19999,1);
HDR.SCP2.HT = [ones(tmp,1),[1:tmp]',HT19999];
end;
elseif section.ID==3,
HDR.NS = fread(fid,1,'char');
HDR.FLAG.Byte = fread(fid,1,'char');
HDR.FLAG.ReferenceBeat = mod(HDR.FLAG.Byte,2);
%HDR.NS = floor(mod(HDR.FLAG.Byte,128)/8);
for k = 1:HDR.NS,
HDR.LeadPos(k,1:2) = fread(fid,[1,2],'uint32');
HDR.Lead(k,1) = fread(fid,1,'uint8');
end;
HDR.N = max(HDR.LeadPos(:))-min(HDR.LeadPos(:))+1;
LeadIdTable = {'I';'II';'V1';'V2';'V3';'V4';'V5';'V6';'V7';'V2R';'V3R';'V4R';'V5R';'V6R';'V7R';'X';'Y';'Z';'CC5';'CM5';'left arm';'right arm';'left leg';'I';'E';'C';'A';'M';'F';'H';'I-cal'};
if max(HDR.Lead)<length(LeadIdTable),
HDR.Label = LeadIdTable(HDR.Lead);
end;
elseif section.ID==4,
HDR.SCP4.L = fread(fid,1,'int16');
HDR.SCP4.fc0 = fread(fid,1,'int16');
HDR.SCP4.N = fread(fid,1,'int16');
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];
tmp = fread(fid,[2*HDR.SCP4.N],'uint32');
HDR.SCP4.PA = reshape(tmp,2,HDR.SCP4.N)';
HDR.SCP4.pa = [0;tmp;HDR.N];
elseif any(section.ID==[5,6]),
SCP = [];
SCP.Cal = fread(fid,1,'int16')/1e6;
SCP.PhysDim = 'mV';
SCP.Dur = fread(fid,1,'int16');
SCP.SampleRate = 1e6/SCP.Dur;
SCP.FLAG.DIFF = fread(fid,1,'uint8');
SCP.FLAG.bimodal_compression = fread(fid,1,'uint8');
SCP.SPR = fread(fid,HDR.NS,'uint16');
if section.ID==6,
HDR.HeadLen = ftell(fid);
HDR.FLAG.DIFF = SCP.FLAG.DIFF;
HDR.FLAG.bimodal_compression = SCP.FLAG.bimodal_compression;
HDR.data = [];
end;
if ~isfield(HDR,'SCP2'),
if any(SCP.SPR(1)~=SCP.SPR),
error('SCPOPEN: SPR do not fit');
else
S2 = fread(fid,[SCP.SPR(1)/2,HDR.NS],'int16');
end;
elseif HDR.SCP2.NHT==19999,
HuffTab = DHT;
S2 = []; %repmat(NaN,HDR.N,HDR.NS);
for k = 1:HDR.NS,
SCP.data{k} = fread(fid,SCP.SPR(k),'uint8');
s2 = SCP.data{k};
s2 = [s2; repmat(0,ceil(max(HDR.SCP2.HT(:,4))/8),1)];
k1 = 0;
l2 = 0;
accu = 0;
c = 0;
x = [];
HT = HDR.SCP2.HT(find(HDR.SCP2.HT(:,1)==1),3:7);
while (l2 < HDR.LeadPos(k,2)),
while c < max(HT(:,2));
k1 = k1 + 1;
% [k1,length(s2),k,HDR.NS],ret
dd = s2(k1);
accu = accu + ACC(dd+1)*(2^c);
c = c + 8;
if 0, %for k2 = 1:8,
accu = accu + (dd>127)*(2^c);
dd = mod(dd*2,256);
c = c + 1;
end;
end;
ixx = 1;
%acc = mod(accu,2^32); % bitand returns NaN if accu >= 2^32
acc = accu - 2^32*fix(accu*(2^(-32))); % bitand returns NaN if accu >= 2^32
while (bitand(acc,2^HT(ixx,1)-1) ~= HT(ixx,5)),
ixx = ixx + 1;
end;
dd = HT(ixx,2) - HT(ixx,1);
if HT(ixx,3)==0,
HT = HDR.SCP2.HT(find(HDR.SCP2.HT(:,1)==HT(ixx,5)),3:7);
fprintf(HDR.FILE.stderr,'Warning SCPOPEN: Switching Huffman Tables is not tested yet.\n');
elseif (dd==0),
l2 = l2 + 1;
x(l2) = HT(ixx,4);
else %if (HT(ixx,3)>0),
l2 = l2 + 1;
%acc2 = fix(accu*(2^(-HT(ixx,1))));
%tmp = mod(fix(accu*(2^(-HT(ixx,1)))),2^dd);
tmp = fix(accu*(2^(-HT(ixx,1)))); % bitshift(accu,-HT(ixx,1))
tmp = tmp - (2^dd)*fix(tmp*(2^(-dd))); % bitand(...,2^dd)
%tmp = bitand(accu,(2^dd-1)*(2^HT(ixx,1)))*(2^-HT(ixx,1));
% reverse bit-pattern
if dd==8,
tmp = ACC(tmp+1);
else
tmp = dec2bin(tmp);
tmp = [char(repmat('0',1,dd-length(tmp))),tmp];
tmp = bin2dec(tmp(length(tmp):-1:1));
end
x(l2) = tmp-(tmp>=(2^(dd-1)))*(2^dd);
end;
accu = fix(accu*2^(-HT(ixx,2)));
c = c - HT(ixx,2);
end;
if k==1,
S2 = x';
elseif size(x,2)==size(S2,1),
S2(:,k) = x';
else
fprintf(HDR.FILE.stderr,'Error SCPOPEN: Huffman decoding failed (%i) \n',size(x,1));
HDR.data = S2;
return;
end;
end;
elseif (HDR.SCP2.NHT==19999), % alternative decoding algorithm.
HuffTab = DHT;
for k = 1:HDR.NS,
SCP.data{k} = fread(fid,SCP.SPR(k),'uint8');
end;
clear S2;
for k = 1:HDR.NS,
tmp = SCP.data{k};
accu = [tmp(4)+256*tmp(3)+65536*tmp(2)+2^24*tmp(1)];
%accu = bitshift(accu,HDR.SCP2.prefix,32);
c = 0; %HDR.SCP2.prefix;
l = 4;
l2 = 0;
clear x;
Ntmp = length(tmp);
tmp = [tmp; zeros(4,1)];
while c <= 32, %1:HDR.SPR(k),
ixx = 1;
while (bitand(accu,mask(ixx)) ~= PREFIX(ixx)),
ixx = ixx + 1;
end;
if ixx < 18,
c = c + prefix(ixx);
%accu = bitshift(accu, prefix(ixx),32);
accu = mod(accu.*(2^prefix(ixx)),2^32);
l2 = l2 + 1;
x(l2) = HuffTab(ixx,1);
elseif ixx == 18,
c = c + prefix(ixx) + 8;
%accu = bitshift(accu, prefix(ixx),32);
accu = mod(accu.*(2^prefix(ixx)),2^32);
l2 = l2 + 1;
acc1 = mod(floor(accu*2^(-24)),256);
%accu = bitshift(accu, 8, 32);
accu = mod(accu*256, 2^32);
x(l2) = acc1-(acc1>=2^7)*2^8;
acc2 = 0;
for kk = 1:8,
acc2 = acc2*2 + mod(acc1,2);
acc1 = floor(acc1/2);
end;
elseif ixx == 19,
c = c + prefix(ixx);
%accu = bitshift(accu, prefix(ixx),32);
accu = mod(accu.*(2^prefix(ixx)),2^32);
l2 = l2 + 1;
while (c > 7) & (l < Ntmp),
l = l+1;
c = c-8;
accu = accu + tmp(l)*2^c;
end;
acc1 = mod(floor(accu*2^(-16)),2^16);
%accu = bitshift(accu, 16, 32);
accu = mod(accu.*(2^16), 2^32);
x(l2) = acc1-(acc1>=2^15)*2^16;
acc2 = 0;
for kk=1:16,
acc2 = acc2*2+mod(acc1,2);
acc1 = floor(acc1/2);
end;
%x(l2) = acc2;
c = c + 16;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -