Commit 69553c45 authored by Enrico Glerean's avatar Enrico Glerean

adding files

parent cda6f9b0
function [hdr, record] = edfreadUntilDone(fname, varargin)
% Read European Data Format file into MATLAB
%
% [hdr, record] = edfread(fname)
% Reads data from ALL RECORDS of file fname ('*.edf'). Header
% information is returned in structure hdr, and the signals
% (waveforms) are returned in structure record, with waveforms
% associated with the records returned as fields titled 'data' of
% structure record.
%
% [...] = edfread(fname, 'assignToVariables', assignToVariables)
% If assignToVariables is true, triggers writing of individual
% output variables, as defined by field 'labels', into the caller
% workspace.
%
% [...] = edfread(...,'desiredSignals',desiredSignals)
% Allows user to specify the names (or position numbers) of the
% subset of signals to be read. |desiredSignals| may be either a
% string, a cell array of comma-separated strings, or a vector of
% numbers. (Default behavior is to read all signals.)
% E.g.:
% data = edfread(mydata.edf,'desiredSignals','Thoracic');
% data = edfread(mydata.edf,'desiredSignals',{'Thoracic1','Abdominal'});
% or
% data = edfread(mydata.edf,'desiredSignals',[2,4,6:13]);
%
% FORMAT SPEC: Source: http://www.edfplus.info/specs/edf.html SEE ALSO:
% http://www.dpmi.tu-graz.ac.at/~schloegl/matlab/eeg/edf_spec.htm
%
% The first 256 bytes of the header record specify the version number of
% this format, local patient and recording identification, time information
% about the recording, the number of data records and finally the number of
% signals (ns) in each data record. Then for each signal another 256 bytes
% follow in the header record, each specifying the type of signal (e.g.
% EEG, body temperature, etc.), amplitude calibration and the number of
% samples in each data record (from which the sampling frequency can be
% derived since the duration of a data record is also known). In this way,
% the format allows for different gains and sampling frequencies for each
% signal. The header record contains 256 + (ns * 256) bytes.
%
% Following the header record, each of the subsequent data records contains
% 'duration' seconds of 'ns' signals, with each signal being represented by
% the specified (in the header) number of samples. In order to reduce data
% size and adapt to commonly used software for acquisition, processing and
% graphical display of polygraphic signals, each sample value is
% represented as a 2-byte integer in 2's complement format. Figure 1 shows
% the detailed format of each data record.
%
% DATA SOURCE: Signals of various types (including the sample signal used
% below) are available from PHYSIONET: http://www.physionet.org/
%
%
% % EXAMPLE 1:
% % Read all waveforms/data associated with file 'ecgca998.edf':
%
% [header, recorddata] = edfread('ecgca998.edf');
%
% % EXAMPLE 2:
% % Read records 3 and 5, associated with file 'ecgca998.edf':
%
% header = edfread('ecgca998.edf','AssignToVariables',true);
% % Header file specifies data labels 'label_1'...'label_n'; these are
% % created as variables in the caller workspace.
%
% Coded 8/27/09 by Brett Shoelson, PhD
% brett.shoelson@mathworks.com
% Copyright 2009 - 2012 MathWorks, Inc.
%
% Modifications:
% 5/6/13 Fixed a problem with a poorly subscripted variable. (Under certain
% conditions, data were being improperly written to the 'records' variable.
% Thanks to Hisham El Moaqet for reporting the problem and for sharing a
% file that helped me track it down.)
%
% 5/22/13 Enabled import of a user-selected subset of signals. Thanks to
% Farid and Cindy for pointing out the deficiency. Also fixed the import of
% signals that had "bad" characters (spaces, etc) in their names.
%
% 10/30/14 Now outputs frequency field directly, and (thanks to Olga Imas
% for the suggestion.)
% HEADER RECORD
% 8 ascii : version of this data format (0)
% 80 ascii : local patient identification
% 80 ascii : local recording identification
% 8 ascii : startdate of recording (dd.mm.yy)
% 8 ascii : starttime of recording (hh.mm.ss)
% 8 ascii : number of bytes in header record
% 44 ascii : reserved
% 8 ascii : number of data records (-1 if unknown)
% 8 ascii : duration of a data record, in seconds
% 4 ascii : number of signals (ns) in data record
% ns * 16 ascii : ns * label (e.g. EEG FpzCz or Body temp)
% ns * 80 ascii : ns * transducer type (e.g. AgAgCl electrode)
% ns * 8 ascii : ns * physical dimension (e.g. uV or degreeC)
% ns * 8 ascii : ns * physical minimum (e.g. -500 or 34)
% ns * 8 ascii : ns * physical maximum (e.g. 500 or 40)
% ns * 8 ascii : ns * digital minimum (e.g. -2048)
% ns * 8 ascii : ns * digital maximum (e.g. 2047)
% ns * 80 ascii : ns * prefiltering (e.g. HP:0.1Hz LP:75Hz)
% ns * 8 ascii : ns * nr of samples in each data record
% ns * 32 ascii : ns * reserved
% DATA RECORD
% nr of samples[1] * integer : first signal in the data record
% nr of samples[2] * integer : second signal
% ..
% ..
% nr of samples[ns] * integer : last signal
if nargin > 5
error('EDFREAD: Too many input arguments.');
end
if ~nargin
error('EDFREAD: Requires at least one input argument (filename to read).');
end
[fid,msg] = fopen(fname,'r');
if fid == -1
error(msg)
end
assignToVariables = false; %Default
targetSignals = []; %Default
for ii = 1:2:numel(varargin)
switch lower(varargin{ii})
case 'assigntovariables'
assignToVariables = varargin{ii+1};
case 'targetsignals'
targetSignals = varargin{ii+1};
otherwise
error('EDFREAD: Unrecognized parameter-value pair specified. Valid values are ''assignToVariables'' and ''targetSignals''.')
end
end
% HEADER
hdr.ver = str2double(char(fread(fid,8)'));
hdr.patientID = fread(fid,80,'*char')';
hdr.recordID = fread(fid,80,'*char')';
hdr.startdate = fread(fid,8,'*char')';% (dd.mm.yy)
% hdr.startdate = datestr(datenum(fread(fid,8,'*char')','dd.mm.yy'), 29); %'yyyy-mm-dd' (ISO 8601)
hdr.starttime = fread(fid,8,'*char')';% (hh.mm.ss)
% hdr.starttime = datestr(datenum(fread(fid,8,'*char')','hh.mm.ss'), 13); %'HH:MM:SS' (ISO 8601)
hdr.bytes = str2double(fread(fid,8,'*char')');
reserved = fread(fid,44);
hdr.records = str2double(fread(fid,8,'*char')');
hdr.duration = str2double(fread(fid,8,'*char')');
% Number of signals
hdr.ns = str2double(fread(fid,4,'*char')');
for ii = 1:hdr.ns
hdr.label{ii} = regexprep(fread(fid,16,'*char')','\W','');
end
if isempty(targetSignals)
targetSignals = 1:numel(hdr.label);
elseif iscell(targetSignals)||ischar(targetSignals)
targetSignals = find(ismember(hdr.label,regexprep(targetSignals,'\W','')));
end
if isempty(targetSignals)
error('EDFREAD: The signal(s) you requested were not detected.')
end
for ii = 1:hdr.ns
hdr.transducer{ii} = fread(fid,80,'*char')';
end
% Physical dimension
for ii = 1:hdr.ns
hdr.units{ii} = fread(fid,8,'*char')';
end
% Physical minimum
for ii = 1:hdr.ns
hdr.physicalMin(ii) = str2double(fread(fid,8,'*char')');
end
% Physical maximum
for ii = 1:hdr.ns
hdr.physicalMax(ii) = str2double(fread(fid,8,'*char')');
end
% Digital minimum
for ii = 1:hdr.ns
hdr.digitalMin(ii) = str2double(fread(fid,8,'*char')');
end
% Digital maximum
for ii = 1:hdr.ns
hdr.digitalMax(ii) = str2double(fread(fid,8,'*char')');
end
for ii = 1:hdr.ns
hdr.prefilter{ii} = fread(fid,80,'*char')';
end
for ii = 1:hdr.ns
hdr.samples(ii) = str2double(fread(fid,8,'*char')');
end
for ii = 1:hdr.ns
reserved = fread(fid,32,'*char')';
end
hdr.label = hdr.label(targetSignals);
hdr.label = regexprep(hdr.label,'\W','');
hdr.units = regexprep(hdr.units,'\W','');
hdr.frequency = hdr.samples./hdr.duration;
disp('Step 1 of 2: Reading requested records. (This may take a few minutes.)...');
if nargout > 1 || assignToVariables
% Scale data (linear scaling)
scalefac = (hdr.physicalMax - hdr.physicalMin)./(hdr.digitalMax - hdr.digitalMin);
dc = hdr.physicalMax - scalefac .* hdr.digitalMax;
% RECORD DATA REQUESTED
tmpdata = struct;
hdr.records = 1e10; %Read a maximum of 1e10 records
dataExists = true;
for recnum = 1:hdr.records
if ~dataExists,break,end
for ii = 1:hdr.ns
% Read or skip the appropriate number of data points
if ismember(ii,targetSignals)
% Use a cell array for DATA because number of samples may vary
% from sample to sample
tempvar = fread(fid,hdr.samples(ii),'int16') * scalefac(ii) + dc(ii);
if isempty(tempvar)
dataExists = false;
hdr.records = recnum-1;
break
else
tmpdata(recnum).data{ii} = tempvar;
end
else
fseek(fid,hdr.samples(ii)*2,0);
end
end
end
hdr.units = hdr.units(targetSignals);
hdr.physicalMin = hdr.physicalMin(targetSignals);
hdr.physicalMax = hdr.physicalMax(targetSignals);
hdr.digitalMin = hdr.digitalMin(targetSignals);
hdr.digitalMax = hdr.digitalMax(targetSignals);
hdr.prefilter = hdr.prefilter(targetSignals);
hdr.transducer = hdr.transducer(targetSignals);
record = zeros(numel(hdr.label), hdr.samples(1)*hdr.records);
% NOTE: 5/6/13 Modified for loop below to change instances of hdr.samples to
% hdr.samples(ii). I think this underscored a problem with the reader.
disp('Step 2 of 2: Parsing data...');
recnum = 1;
for ii = 1:hdr.ns
if ismember(ii,targetSignals)
ctr = 1;
for jj = 1:hdr.records
try
record(recnum, ctr : ctr + hdr.samples(ii) - 1) = tmpdata(jj).data{ii};
end
ctr = ctr + hdr.samples(ii);
end
recnum = recnum + 1;
end
end
hdr.ns = numel(hdr.label);
hdr.samples = hdr.samples(targetSignals);
if assignToVariables
for ii = 1:numel(hdr.label)
try
eval(['assignin(''caller'',''',hdr.label{ii},''',record(ii,:))'])
end
end
% Uncomment this line to duplicate output in a single matrix
% ''record''. (Could be memory intensive if dataset is large.)
record = [];% No need to output record data as a matrix?
end
end
fclose(fid);
\ No newline at end of file
clear all
close all
rng(0)
a=randn(100,10);
%% matrix or struct
a=randn(100,10);
s='hello planet';
x=[];
x.time=1:100;
x.amplitude=sin(1:100);
x.str=s;
\ No newline at end of file
clear all
close all
rng(0)
a=randn(100,10);
%% matrix or struct
a=randn(100,10);
s='hello planet';
x=[];
x.time=1:100;
x.amplitude=sin(1:100);
x.str=s;
%% the commands I have typed the most in my life
whos
class(a)
class(x) size(s)
size(a) size(x)
x
a(1)
disp("We are here at the code")
%error('stop everything')
disp(num2str(a(1),2))
sprintf('%05d',10)
\ No newline at end of file
function out=arrProd(s,v);
% Make a function that takes as input parameters:
% A scalar
% A vector of size N
% Checks that the first is a scalar
if(length(s)~=1)
error('first input is not a scalar');
end
% Checks that the second is a vector
if(~(length(v)>1 && any(size(v)==1) && length(size(v)==2)))
error('second input is not a vector');
end
% Multiplies each element of the vector with a for loop
N=length(v);
out=zeros(N,1);
for n=1:N
out(n)=s*v(n);
end
\ No newline at end of file
/*==========================================================
* arrayProduct.c - example in MATLAB External Interfaces
*
* Multiplies an input scalar (multiplier)
* times a 1xN matrix (inMatrix)
* and outputs a 1xN matrix (outMatrix)
*
* The calling syntax is:
*
* outMatrix = arrayProduct(multiplier, inMatrix)
*
* This is a MEX-file for MATLAB.
* Copyright 2007-2012 The MathWorks, Inc.
*
*========================================================*/
#include "mex.h"
/* The computational routine */
void arrayProduct(double x, double *y, double *z, mwSize n)
{
mwSize i;
/* multiply each element y by x */
for (i=0; i<n; i++) {
z[i] = x * y[i];
}
}
/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
double multiplier; /* input scalar */
double *inMatrix; /* 1xN input matrix */
size_t ncols; /* size of matrix */
double *outMatrix; /* output matrix */
/* check for proper number of arguments */
if(nrhs!=2) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs","Two inputs required.");
}
if(nlhs!=1) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs","One output required.");
}
/* make sure the first input argument is scalar */
if( !mxIsDouble(prhs[0]) ||
mxIsComplex(prhs[0]) ||
mxGetNumberOfElements(prhs[0])!=1 ) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notScalar","Input multiplier must be a scalar.");
}
/* make sure the second input argument is type double */
if( !mxIsDouble(prhs[1]) ||
mxIsComplex(prhs[1])) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notDouble","Input matrix must be type double.");
}
/* check that number of rows in second input argument is 1 */
if(mxGetM(prhs[1])!=1) {
mexErrMsgIdAndTxt("MyToolbox:arrayProduct:notRowVector","Input must be a row vector.");
}
/* get the value of the scalar input */
multiplier = mxGetScalar(prhs[0]);
/* create a pointer to the real data in the input matrix */
inMatrix = mxGetPr(prhs[1]);
/* get dimensions of the input matrix */
ncols = mxGetN(prhs[1]);
/* create the output matrix */
plhs[0] = mxCreateDoubleMatrix(1,(mwSize)ncols,mxREAL);
/* get a pointer to the real data in the output matrix */
outMatrix = mxGetPr(plhs[0]);
/* call the computational routine */
arrayProduct(multiplier,inMatrix,outMatrix,(mwSize)ncols);
}
# include "mex.h"
void mexFunction ( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] );
/**********************************************************************/
void mexFunction ( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
/**********************************************************************/
/*
Purpose:
MEXFUNCTION is a MATLAB/C interface for the factorial function.
Discussion:
This file should be called "fact.c". It should be placed in the
MATLAB user's path. It can either be compiled externally, with
a command like
mex fact.c
creating a compiled MEX file, or, inside of MATLAB, the command
mex fact.c
accomplishes the same task. Once the file has been compiled,
The MATLAB user can invoke the function by typing:
y = fact ( x )
The routine carries out the computation of the factorial of X,
whose value is returned in Y.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
17 July 2006
Author:
Duane Hanselman and Bruce Littlefield
Reference:
Duane Hanselman, Bruce Littlefield,
Mastering MATLAB 7,
Pearson Prentice Hall, 2005,
ISBN: 0-13-143018-1.
The Mathworks,
MATLAB External Interfaces,
The Mathworks, 2000.
Parameters:
Input, int NLHS, the number of output arguments.
Input, mxArray *PLHS[NLHS], pointers to the output arguments.
Input, int NRHS, the number of input arguments.
Input, const mxArray *PRHS[NRHS], pointers to the input arguments.
*/
{
int i;
double x;
double y;
double *y_pointer;
/*
INPUT:
Retrieve the (first) (scalar) input argument from the line
y = fact ( x )
*/
x = mxGetScalar ( prhs[0] );
/*
COMPUTATION:
Here is where the computation is done.
In many cases, it would make sense to have this computation be in
separate C function, so that the mexFunction is really just an
interface.
*/
y = 1.0;
for ( i = 1; i <= ( int ) x; i++ )
{
y = y * ( double ) i;
}
/*
OUTPUT:
Make space for the output argument,
retrieve the value of the pointer to that space,
and copy the value of Y into the address for the output.
*/
plhs[0] = mxCreateDoubleMatrix ( 1, 1, mxREAL );
y_pointer = mxGetPr ( plhs[0] );
*y_pointer = y;
return;
}
\ No newline at end of file
clear all
close all
nii=load_nii('ch2better.nii');
%%
if(0)
for z=1:size(nii.img,3)
imshow(nii.img(:,:,z))
imwrite(nii.img(:,:,z),['slices/' num2str(z) '.png'])
pause(0.01)
colormap(gray)
end
end
%%
X=size(nii.img,1);
Y=size(nii.img,2);
Z=size(nii.img,3);
X=301;
Y=370;
Z=316;
fileID = fopen('mybinary.bin','w');
for z=1:Z
disp(num2str(z))
data=imread(['slices/' num2str(z) '.png']);
fwrite(fileID,data,'uint8');
end
fclose(fileID)
error('stop')
%%
m = memmapfile('mybinary.bin','Format',{'uint8',[X Y Z],'x'})
save memMap m
%% and then
clear all
load memMap
% see that the size it uses is small
imshow(squeeze(m.Data.x(170,:,:)))
clear all
close all
N=100;
s=2;
v=1:10e6;
timeElapsed=zeros(N,1);
timeElapsedMEX=zeros(N,1);
for n=1:N
disp(num2str(n))
tic;
%out=arrProd(s,v);
out=s*v;
timeElapsed(n)=toc;
tic
out=arrayProduct(s,v);
timeElapsedMEX(n)=toc;
end
%%
figure
histogram(timeElapsed,20)
hold on
histogram(timeElapsedMEX,20)
\ No newline at end of file
clear all
close all
Z=316;
for z =1:Z
vol(:,:,z)=double(imread(['slices/' num2str(z) '.png']));
end
\ No newline at end of file
clear all
close all
rng(0)
t=(0:0.01:1000)';
x=(20*randn(length(t),1));
y=(20*randn(length(t),1));
figure(1)
plot(x)
%%
figure
plot(y)
%%
figure(3)
plot(x)
hold on
plot(y)
%%
figure
plot(x,y,'o')
%%
x=x+100*cos(t);
y=y+100*sin(t);
X=x-min(x);
Y=y-min(y);
M=ceil(max([X(:);Y(:)]));
hmap=zeros(M+1);
for t=1:length(X)
hmap(round(X(t))+1,round(Y(t))+1)=1+hmap(round(X(t))+1,round(Y(t))+1);
end
%%
figure
imagesc(hmap)
colormap(hot)
error('stop')
%%
close all
rng(0)
x=zscore(cumsum(randn(100,1)))
y=100*zscore(cumsum(randn(100,1)));
yyaxis left
plot(x)
yyaxis right
plot(y)
\ No newline at end of file
close all
clear all
rng(0);
T=1e5;
IF=10;
x=IF*randn(T,1);
y=IF*randn(T,1);
plot(x,y,'o')
x=x+100*cos((1:T)'/1e4);
y=y+100*sin((1:T)'/1e4);
x=x-min(x);
y=y-min(y);
M=max([x;y]);
data=zeros(ceil(M)+1);
for t=1:T
data(round(x(t))+1,round(y(t))+1)=data(round(x(t))+1,round(y(t))+1)+1;
end
imagesc(data)
map=parula(5);
map=hot
colormap(map)
colorbar
%%
clear all
close all
T=1e5;
x=zscore(cumsum(randn(T,1)));
X=fft(x,1024);
M=abs(X);
M(513:end)=[];
plot(10*log10(M))
ax=gca;
ax.XScale='log'