Uživatel:Kychot/emg osc.m

Z Wikiverzity
Skočit na navigaci Skočit na vyhledávání
# VERZE PRO OCTAVE 3.0.1
# posouvání podle kláves

% EMG processing
% EMG device: KeyPoint
% Structure of data files OSCILLO
% avg - průměry
% ref = přepočítávají se reference
% rost(rální) = ch2 - ch1
% med(iální) = ch2 - ch4 - ch3
% lat(erální) = ch2 - ch4
% kaud(ální) = ch2
%
%
%    sampling frequency =  5000 Hz
%    one sample = 2 bytes, i.e. signed short integer, Little-endian
%    amplitude: 1 bit = 1 uV, A/D conversion 10 bit
%    max. amplitude = 0xfe0c = 65536 - 65036 = -500 uV ... +500 uV

% Header for 6 channels:
% 0000000      6    200     50     50     50     50     50    500
% 0000020   -206   -202   -208   -205   -204    -21   -205   -204

global emgdir;
#emgdir = "/sda/7/huge/EMG/";	# @INA - šuplík
#emgdir = "/sdb/7/huge/EMG/";	# @HRC - šuplík
#emgdir = "/home/petr/EMG/CD/"; # @mi
#emgdir = "/media/data-60G/petr/EMG/CD/";   # @mi
emgdir = "./";	# current


global datasubdir;

global outdir;
outdir = "./osc-out/";
#outdir = "osc-out/";

global rowsize; %samples
#    rowsize = 8000;  % pro EM5
# Kalibrovano na brumu 50Hz, soubor CD01/A/A0003154.EXP:
#    rowsize = 10000;  %* 100ms/div = 500ms/page @ 20kHz pro OSC BOOST 
#    rowsize = 5000;  %  200ms/div =    1s/page @ 5kHz pro OSC 
#    rowsize = 4000;  %  100ms/div = 800ms/page @ 5kHz pro OSC 
    rowsize = 2500;  %* 100ms/div = 500ms/page @ 5kHz pro OSC 
#    rowsize = 2000;  %  100ms/div = 400ms/page @ 5kHz pro OSC 
#    rowsize = 1000;  %   40ms/div = 200ms/page @ 5kHz pro OSC
#    rowsize =  500;  %   20ms/div = 100ms/page @ 5kHz pro OSC
#    rowsize =  250;  %   10ms/div =  50ms/page @ 5kHz pro OSC
#    rowsize =  200;  %   10ms/div =  40ms/page @ 5kHz pro OSC
#    rowsize =  100;  %    4ms/div =  20ms/page @ 5kHz pro OSC
#    rowsize =   50;  %    2ms/div =  10ms/page @ 5kHz pro OSC
# Kalibrace: 200Hz square dle manualu, ve skutecnosti vypada na 100Hz

global sampling_frequency = 5000; # Hz

global nchans;

function retval = timestr (nseconds_float)
%%%%%%%%          %%%%%%%%%%%%%%%%%%%%%%%%
    nseconds_total = floor(nseconds_float);
    nmilliseconds = (nseconds_float - nseconds_total) * 1000;
    nminutes_float = nseconds_total / 60;
    nminutes = floor(nminutes_float);
    nseconds = nseconds_total - 60*nminutes;
    retval = sprintf("%d:%02d.%03d", nminutes, nseconds, nmilliseconds);
endfunction


function emgplot (page)
%%%%%%%%%%%%%%%%%%%%%%%%
    global headerstr signal dfnamestr rowsize nchans outdir sample_interval_us sensitivity;

    sampling_freq_kHz = 1000 / sample_interval_us;


#    figfilename = strcat (outdir, dfnamestr, "-", sprintf("%04d",page), ".fig")


#    __gnuplot_set__  terminal fig
%    __gnuplot_set__  terminal fig color size 6 7.2 solid textnormal fontsize 16
%    __gnuplot_set__  terminal fig color size 11.69 8.27 solid textnormal fontsize 16
%    __gnuplot_set__  output "'/home/petr/EMG/octave/out/test7.fig'"
%    __gnuplot_raw__  figfilenm = figfilename
%    __gnuplot_set__  output "'/home/petr/EMG/octave/out/tmp.fig'"
%    __gnuplot_set__  output "'/huge/EMG/octave/osc-out/tmp.fig'"
#    __gnuplot_set__  terminal fig color size 10.91 7.48 solid textnormal fontsize 16
%    __gnuplot_set__  terminal fig monochrome size 10.91 7.48 solid textnormal fontsize 16
# na HRC:
%    __gnuplot_set__  output "'./osc-out/tmp.fig'"
# na INA - kyž bývala stará verze 2.1:
%    __gnuplot_set__  output "./osc-out/tmp.fig"
#    __gnuplot_set__  output "/tmp/tmp.fig"	

#    shift = 20;    __gnuplot_set__ ytics -20,20,240;
#    shift = 40;    __gnuplot_set__ ytics -40,40,480;
#    shift = 50;    __gnuplot_set__ ytics -50,50,600;
#    shift = 100;  __gnuplot_set__ ytics -100,100,1200;
#    shift = 200;  __gnuplot_set__ ytics -200,200,2400;
%    shift = 500;   __gnuplot_set__ ytics -500,500,6000;
    shift = 500;
#    shift = 1000;  __gnuplot_set__ ytics -1000,1000,12000;
#    shift = 2000;  __gnuplot_set__ ytics -2000,2000,24000;
#    shift = 5000;  __gnuplot_set__ ytics -5000,5000,60000;
#    shift = 5000;  __gnuplot_set__ ytics -5000,5000,70000; # 2 sety

%    __gnuplot_set__ ytics autofreq;


    titlestr = strcat (sprintf("%d chan at %d kHz, ", nchans, sampling_freq_kHz), sprintf("%d ",sensitivity), "uV/Div=409.6Lev, ", "500samp/div, 100ms/div   ",headerstr, "-", sprintf("%04d",page), "    ", timestr((page-1)*0.5));
    title (titlestr);
%    xlabel (strcat("1 chan,  U [", int2str(shift), "level/div = _____ uV/div]   -> t [50 ms/div]"));

    axis([0, rowsize, -shift, 12*shift], "nolabel");
    grid ("on");
%    plot_border ("blank");
    hold on;

    for s = 1:rowsize
	for ch = 1:nchans
	    channel(ch,s) = signal(nchans*(s-1)+ch);
	endfor
    endfor

% transformace referencí pro 4 kanály:

#    for s = 1:rowsize
#	  chanref(1,s) = channel(2,s) - channel(1,s);	% rost(rální) = ch2 - ch1
#	  chanref(2,s) = channel(2,s) - channel(4,s) - channel(3,s);	% med(iální) = ch2 - ch4 - ch3
#	  chanref(3,s) = channel(2,s) - channel(4,s);	% lat(erální) = ch2 - ch4
#	  chanref(4,s) = channel(2,s);			% kaud(ální) = ch2
#    endfor

% průměry referencí pro 4 kanály:

#    for s = 1:rowsize
#	  chanavg(1,s) = (chanref(1,s) + chanref(4,s))/2;	% rost+kaud
#	  chanavg(2,s) = (chanref(2,s) + chanref(3,s))/2;	% med+lat
#	  chanavg(3,s) = (chanavg(1,s) + chanavg(2,s))/2;	% all
#	  chanavg(4,s) = chanref(4,s);				% stejný=kaud
#    endfor


% průměry referencí pro 2+2 kanály
% ipsilat: 1=plošná, 2=monopol; kontralat:3=plošná, 4=monopolár.; společná ref. mezi:
% příklad, když je TrP vpravo - záda ze shora, * je plošná, + je monopolární:
% l.sin.     *3    +4     *ref.   +2    *1      l.dx. 
#    for s = 1:rowsize
#	  chan2avg(1,s) = channel(1,s);			% plošná ipsi
#	  chan2avg(2,s) = channel(2,s)-channel(1,s)/2;	% monopol. ipsim, prům.ref.
#	  chan2avg(3,s) = channel(3,s);			% plošná kontra
#	  chan2avg(4,s) = channel(4,s)-channel(3,s)/2;	% monopol. kontra, prům. ref.
#    endfor

    for ch = 1:nchans
#        offset = (12 - 2*ch)*shift;   # for 6 channels
        offset = (13 - 3*ch)*shift;   # for 4 channels
        plot ([1:rowsize], channel(ch,:)+offset, "k-;;");	% původní
#        plot ([1:rowsize], chanref(ch,:)+offset, "k-;;");	% přepočítaná reference
#        plot ([1:rowsize], chanavg(ch,:)+offset, "k-;;");	% průměry ze čtyř
#        plot ([1:rowsize], chan2avg(ch,:)+offset, "k-;;");	% průměry ze 2+2
    endfor

    hold off;
#    __gnuplot_set__  terminal x11
#    __gnuplot_set__  output "'/dev/null'"

%    [err, msg] = rename ("/home/petr/EMG/octave/out/tmp.fig", figfilename)
%    [err, msg] = rename ("./osc-out/tmp.fig", figfilename)
#    err
#    msg

endfunction


function dfopen (dfname)
%%%%%%%%%%%%%%%%%%%%%%%%
    global emgdir datasubdir headerstr dfnamestr signal rowsize nchans sample_interval_us sampling_frequency_kHz sensitivity;

#    emgdir
#    datasubdir
#    dfname
    dfpathname = strcat (emgdir, datasubdir, dfname)
    [fstat,err,msg] = stat(dfpathname);
    filetimestr = strftime("%a %Y-%m-%d %R", localtime(fstat.mtime));
#    filetimestr = 'ZdC G,vys. z 2007-26-09, pocitano 2007-10-17';
    headerstr = strcat (filetimestr, "  ", datasubdir, dfname);
    dfnamestr = dfname;

    datafile = fopen (dfpathname, "r", "ieee-le");

    [nchans, count] = fread (datafile,  1, "int16");
    [sample_interval_us, count] = fread (datafile,  1, "int16");
    [sensitivity, count] = fread (datafile,  nchans, "int16");
    sensitivity;
    sensitivity = sensitivity.';
    headersize = 4 + 2*nchans;   % bytes
    datasize = (fstat.size - headersize)/2; % samples
    nsamples = datasize / nchans
#    nseconds_float = nsamples / sampling_frequency
    nseconds_float = nsamples * sample_interval_us / 1000000
    timestr (nseconds_float)
    pagedatasize = nchans*rowsize;    % samples  pro 6-channel 200ms OSC @ 5kHz 
    npages = datasize / pagedatasize

#   for page=1300:120:npages		# po 1 min při 5 kHz nebo 15s při 20 kHz
#   for page=19000:40:npages		# po 5 sec při 20 kHz
#   for page=2000:240:npages		# po 30s při 20 kHz
#   for page=1:480:npages		# po 1 min při 20 kHz
#   for page=1:120:npages		# po 1 min při 5 kHz nebo 15s při 20 kHz
#   for page=8520:40:npages		# po 1 min 2E
#   for page=1:120:npages		# po 1 min při 5 kHz nebo 15s při 20 kHz
#   for page=1:60:npages		# po 30 s
#   for page=1:30:npages		# po 15 s
#   for page=1:20:npages		# po 10 s
%   for page=1:10:npages		# po 5 s
#   for page=1:30:209 # 000A3234.EXP
#   for page=1:60:1723 # 000B3234.EXP
%   for page=1:1
%    for page=1:nchunks
#   for page=4800:8:npages		# po 1 sec při 20 kHz, od   

     page=1;
     while (page>0 && page<npages)

	fseek (datafile, headersize + (page-1)*2*pagedatasize, SEEK_SET);
        position=ftell (datafile)

	[signal, count] = fread (datafile,  pagedatasize, "int16");

	clf;
	emgplot(page);

	page
#	disp("5m=A 1m=S 10s=D 2s=F 500ms=G 125ms=T <> Y=125ms H=500ms J=2s K=10s L=1m p=5m")
	disp("5m=A 1m=S 10s=D 2s=F 500ms=G  <> H=500ms J=2s K=10s L=1m p=5m")
	fflush(1);
#	input("Strike any key for next page");
	o_kolik = kbhit();
#	o_kolik
	switch(o_kolik)
	    case 'a'
	      page -= 600;
	    case 's'
	      page -= 120;
	    case 'd'
	      page -= 20;
	    case 'f'
	      page -= 4;
	    case 'g'
	      page -= 1;
	    case 'h'
	      page += 1;
	    case 'j'
	      page += 4;
	    case 'k'
	      page += 20;
	    case 'l'
	      page += 120;
	    case 'p'
	      page += 600;
#	    case 'a'
#	      page -= 2400;
#	    case 's'
#	      page -= 480;
#	    case 'd'
#	      page -= 80;
#	    case 'f'
#	      page -= 16;
#	    case 'g'
#	      page -= 4;
#	    case 't'
#	      page -= 1;
#	    case 'y'
#	      page += 1;
#	    case 'h'
#	      page += 4;
#	    case 'j'
#	      page += 16;
#	    case 'k'
#	      page += 80;
#	    case 'l'
#	      page += 480;
#	    case 'p'
#	      page += 2400;
	endswitch
#	page += 480;
    endwhile
    fclose (datafile);
endfunction



#datasubdir = "CD01/A/";
				# 2007-06-06:
#dfopen ("00003154.EXP");
#dfopen ("A0003154.EXP");
#dfopen ("B0003154.EXP");
#dfopen ("C0003154.EXP");
#dfopen ("D0003154.EXP");
				# 2007-06-26:
#dfopen ("00003220.EXP");
#dfopen ("A0003220.EXP");
#dfopen ("B0003220.EXP");
#dfopen ("C0003220.EXP");
#dfopen ("D0003220.EXP");
#dfopen ("F0003220.EXP");
#dfopen ("G0003220.EXP");

#datasubdir = "CD01/B/";  # Wed 2007-06-27 PHe
#dfopen ("2E003220.EXP"); #
#dfopen ("3A003220.EXP"); #
#dfopen ("3B003220.EXP"); #
#dfopen ("3C003220.EXP"); #
#dfopen ("3D003220.EXP"); #
#dfopen ("3E003220.EXP"); #
#dfopen ("3F003220.EXP"); #

# C
#datasubdir = "CD03/C/";  # 2007-07-18 
#dfopen ("000A3234.EXP"); # 
#dfopen ("000B3234.EXP"); # 
#dfopen ("000C3234.EXP"); # 

#dfopen ("0003276A.EXP"); # ZdC

#datasubdir = "CD03/D/";  # 2007-07-25 ZC 4
#dfopen ("400A3248.EXP"); #
#dfopen ("400B3248.EXP"); #
#dfopen ("400C3248.EXP"); #
#dfopen ("400D3248.EXP"); #
#dfopen ("400E3248.EXP"); #
#dfopen ("400F3248.EXP"); #

#datasubdir = "CD03/E/";  # 2007-08-29
#dfopen ("0003276A.EXP"); #

#datasubdir = "CD04/F/";   # 2007-09-19
#dfopen ("600A3336.EXP"); #
#dfopen ("600B3336.EXP"); #
#dfopen ("600C3336.EXP"); #

#datasubdir = "CD03/G/";  # 2007-09-26 ZC
#dfopen ("700A3358.EXP"); #
#dfopen ("700B3358.EXP"); #
#dfopen ("700C3358.EXP"); #

#datasubdir = "CD03/H/";  # 2007-10-17
#datasubdir = "CD04/H/";  # 2007-10-17
#dfopen ("800A3427.EXP"); #
#dfopen ("800B3427.EXP"); #
#dfopen ("800C3427.EXP"); #

# datasubdir = "CD04/I/";  # 2007-10-24
# dfopen ("900A3449.EXP"); #
# dfopen ("900B3449.EXP"); #
# dfopen ("900C3449.EXP"); #
# dfopen ("900D3449.EXP"); #

# datasubdir = "CD04/J/";  # 2007-11-07
# dfopen ("A00A3496.EXP"); #
# dfopen ("A00B3496.EXP"); #

# K
# datasubdir = "CD04/K/";  # 2007-11-21 . , testování
# dfopen ("000A3533.EXP"); #
# dfopen ("000B3533.EXP"); #
# dfopen ("000C3533.EXP"); #
# dfopen ("000D3533.EXP"); #
# dfopen ("000E3533.EXP"); #
# dfopen ("000F3533.EXP"); #
# dfopen ("000G3533.EXP"); #

# L
# datasubdir = "CD04/L/";  # 2007-12-05 PHo.
# dfopen ("C00A3573.EXP"); #
# dfopen ("C00B3573.EXP"); #

# datasubdir = "CD04/M/";  # 2007-12-19 
# dfopen ("D00A3611.EXP"); #
# dfopen ("D00B3611.EXP"); #
# dfopen ("D00C3611.EXP"); #
# dfopen ("D00D3611.EXP"); #

# N
# datasubdir = "CD04/N/";  # 2008-01-09 .
# dfopen ("E00A3635.EXP"); # Začátek 50 uV/div
# dfopen ("E00B3635.EXP"); #  Hlavní soubor
# dfopen ("E00C3635.EXP"); #
# dfopen ("E00D3635.EXP"); #
# dfopen ("E00E3635.EXP"); #


# datasubdir = "CD03/O/";  # 2008-01-06  nebo MV.
# dfopen ("F00A3649.EXP"); # hlavní
# dfopen ("F00B3649.EXP"); # filtry
# dfopen ("F00C3649.EXP"); # vytahování jehly

# P
# datasubdir = "CD05/P/";  # 2008-03-12 AH. 1956
# dfopen ("G00A3790.EXP"); #
# dfopen ("G00B3790.EXP"); #

# R
# datasubdir = "CD05/R/";  # 2008-03-19 ZM. 1985
# dfopen ("H00A3801.EXP"); #
# dfopen ("H00B3801.EXP"); #

# S
# datasubdir = "CD05/S/";  # 2008-03-26 LH. 1960
# dfopen ("I00A3818.EXP"); #
# dfopen ("I00A3818.EXP"); #

# T
# datasubdir = "CD05/T/";  # 2008-04-02 KP. 
# dfopen ("J00A3840.EXP"); #
# dfopen ("J00B3840.EXP"); #

#
# datasubdir = "CD05/U/";  # 2008-04-16
# dfopen ("K00A3877.EXP"); #

#
# datasubdir = "CD06/V/";  # 2008-06-18
# dfopen ("L00A4047.EXP"); #
# dfopen ("L00B4047.EXP"); #
# dfopen ("L00C4047.EXP"); #

# X
# datasubdir = "CD06/X/";  # 2008-06-25 PHe. 1953
# dfopen ("M00A4073.EXP"); #
# dfopen ("M00B4073.EXP"); #
# dfopen ("M00C4073.EXP"); #

# Y
# datasubdir = "CD07/Y/";  # 2008-07-02 .
# dfopen ("N00A4104.EXP"); #
# dfopen ("N00B4104.EXP"); #
# dfopen ("N00C4104.EXP"); #

# Z
# datasubdir = "CD??/Z/";  # 2008-07-16 LM.
# O00A4135 chybí !

#
# datasubdir = "CD07/Z1/";  # 2008-09-24
# dfopen ("P00A4254.EXP"); #
# dfopen ("P00B4254.EXP"); #
# dfopen ("P00C4254.EXP"); #

#
# datasubdir = "CD07/Z2-pokus/";  # 2008-10-15
# dfopen ("XXX00001.EXP"); #
# dfopen ("XXX00002.EXP"); #
# dfopen ("XXX00004.EXP"); #
# dfopen ("XXX00008.EXP"); #

#
# datasubdir = "CD07/";  # 2008-10-22
# dfopen ("XXXX0000.EXP"); #

#
# datasubdir = "CD08/2A/"; # 2008-10-22
# dfopen ("2AA04336.EXP"); #

#
# datasubdir = "CD09/2B/"; # 2008-10-29 ??protokol??
# dfopen ("2BCh4347.EXP"); #          (~9.8 MB)
# dfopen ("2BSi4347.EXP"); #          (~116 MB)

# 2C
# datasubdir = "CD09/2C/"; # 2008-11-05 ŠB. 1983
# dfopen ("2CSi4370.EXP"); # 32:00 min (298.2 MB) celé, . vytahování

#
# datasubdir = "CD10/2D/"; # 2008-11-19 ZS. 1983-01
# dfopen ("2DA04403.EXP"); # 03:00 min (29.2 MB), zač.
# dfopen ("2DB04403.EXP"); # 33:05 min (317.6 MB)  vše

# 2E
# datasubdir = "CD11/2E/"; # 2008-11-26 JL. 1964-07
# dfopen ("2EA04425.EXP"); # 02:15 min (20.7 MB)
# dfopen ("2EB04425.EXP"); # 04:10 min (40.00 MB)
# dfopen ("2EC04425.EXP"); # kousek (0.1 MB)
# dfopen ("2ED04425.EXP"); # 02:10 min (20.8 MB) AP
# dfopen ("2EE04425.EXP"); # 30:00 min (288.0 MB) vše
# dfopen ("2EF04425.EXP"); # 06:19 min (55.5? MB) různá zesílení, vytahování el. 

# 2F
# datasubdir = "CD12/2F/"; # 2008-12-17 LS. 1986-07
# dfopen ("2FA4490.EXP"); # 02:40 min (   MB) začátek s kalibrací  
# dfopen ("2FB4490.EXP"); # 02:45 min (   MB) začátek
# dfopen ("2FC4490.EXP"); # >37 min   (385 MB) hlavní část a konec, vytahování

# 2G
# datasubdir = "CD13/2G/"; # 2009-01-14 JB. 1984-08
 datasubdir = ""; 	   # current dir.
# dfopen ("2GA04529.EXP"); # 30:00 min (288.2 MB) celý záznam, bez vytahování
 dfopen ("2GA04529q.EXP"); # 30:00 min (288.2 MB) celý záznam, 4x redukovaný vzork. kmitočet