Contents
- Gain from axon binary file
- Load data from axon binary file
- Controller for abfload function
- Filter waveform data
- Cropp waveform
- Peak detection
- Event Detection
- Flag Bad Events
- Generate Descriptive Stats
- Exponential decay fit
- Expoential rise fit
- Resample data version1
- Resample data version1
- Get relative data paths to data files
Gain from axon binary file
back to top%--------------------------------- function datagain = abfgain(abf_file)
% gain = abfgain(abf_filename) % % input % abf_file = name of abf file on path % % output % datagain = amplitude gain %Search the file for the Gain information myfileN = abf_file; %enter file name fid = fopen(myfileN,'r','ieee-le'); %open the abf file sz=1; %# of data entries to read numType = 'float'; %precision format myGainfileoffset = 268; %find the correct offset where data entry for gain begins fseek(fid,myGainfileoffset,'bof'); %move field position datagain = fread(fid,sz,numType); %read formatted data fclose(fid); %close file %---------------------------------
Load data from axon binary file
back to top%--------------------------------- function [d,si]=abfload(fn,varargin)
% ** function [d,si]=abfload(fn,varargin) % loads and returns data in the Axon abf format PRIOR TO VERSION 2.0. % Data may have been acquired in the following modes: % (1) event-driven variable-length % (2) event-driven fixed-length % (3) gap-free % Information about scaling, the time base and the number of channels and % episodes is extracted from the header of the abf file (see also abfinfo.m). % All optional input parameters listed below (= all except the file name) % must be specified as parameter/value pairs, e.g. as in % d=abfload('d:\data01.abf','start',100,'stop','e'); % % >>> INPUT VARIABLES >>> % % NAME TYPE, DEFAULT DESCRIPTION % fn char array abf data file name % start scalar, 0 only gap-free-data: start of cutout to be read (unit: sec) % stop scalar or char, only gap-free-data: end of cutout to be read (unit: sec). % 'e' May be set to 'e' (end of file). % sweeps 1d-array or char, only episodic data: sweep numbers to be read. By default, % 'a' all sweeps will be read ('a') % channels cell array names of channels to be read, like {'IN 0','IN 8'}; % or char, 'a' ** make sure spelling is 100% correct (including blanks) ** % if set to 'a', all channels will be read % chunk scalar, 0.05 only gap-free-data: the elementary chunk size (Megabytes) % to be used for the 'discontinuous' mode of reading data % (when less channels shall be read than exist in file) % machineF char array, the 'machineformat' input parameter of the % 'ieee-le' matlab fopen function. 'ieee-le' is the correct % option for windows; depending on the % platform the data were recorded/shall be read % by abfload 'ieee-be' is the alternative. % % <<< OUTPUT VARIABLES <<< % % NAME TYPE DESCRIPTION % d the data read, the format depending on the recording mode % 1. GAP-FREE: % 2d array 2d array of size % '<data pts>' by '<number of chans>' % Examples of access: % d(:,2) data from channel 2 at full length % d(1:100,:) first 100 data points from all channels % 2. EPISODIC FIXED-LENGTH: % 3d array 3d array of size % '<data pts per sweep>' by '<number of chans>' by '<number of sweeps>' % Examples of access: % d(:,2,:) is a matrix which contains all events (at full length) % of channel 2 in its columns % d(1:200,:,[1 11]) contains first 200 data points of events #1 % and #11 of all channels % 3. EPISODIC VARIABLE-LENGTH: % cell array cell array whose elements correspond to single sweeps. Each element is % a (regular) array of size % '<data pts per sweep>' by '<number of chans>' % Examples of access: % d{1} a 2d-array which contains sweep #1 (all of it, all channels) % d{2}(1:100,2) a 1d-array containing first 100 data points of of channel 2 in sweep #1 % % si scalar the sampling interval in usec % % Date of this version: June 17, 2011 % future improvements: % - handle expansion of header in transition to file version 1.65 better % --- defaults % gap-free start=0.0; stop='e'; % episodic sweeps='a'; % general channels='a'; % the size of data chunks (see above) in Mb. 0.05 Mb is an empirical value % which works well for abf with 6-16 channels and recording durations of % 5-30 min chunk=0.05; machineF='ieee-le'; verbose=1; % assign values of optional input parameters, if any were given pvpmod(varargin); if verbose, disp(['** ' mfilename]); end d=[]; si=[]; if ischar(stop) if ~strcmpi(stop,'e') error('input parameter ''stop'' must be specified as ''e'' (=end of recording) or as a scalar'); end end % --- obtain vital header parameters and initialize them with -1 % temporary initializing var tmp=repmat(-1,1,16); % for the sake of typing economy set up a cell array (and convert it to a struct below) % column order is % name, position in header in bytes, type, value) headPar={ 'fFileVersionNumber',4,'float',-1; 'nOperationMode',8,'int16',-1; 'lActualAcqLength',10,'int32',-1; 'nNumPointsIgnored',14,'int16',-1; 'lActualEpisodes',16,'int32',-1; 'lFileStartTime',24,'int32',-1; 'lDataSectionPtr',40,'int32',-1; 'lSynchArrayPtr',92,'int32',-1; 'lSynchArraySize',96,'int32',-1; 'nDataFormat',100,'int16',-1; 'nADCNumChannels', 120, 'int16', -1; 'fADCSampleInterval',122,'float', -1; 'fSynchTimeUnit',130,'float',-1; 'lNumSamplesPerEpisode',138,'int32',-1; 'lPreTriggerSamples',142,'int32',-1; 'lEpisodesPerRun',146,'int32',-1; 'fADCRange', 244, 'float', -1; 'lADCResolution', 252, 'int32', -1; 'nFileStartMillisecs', 366, 'int16', -1; 'nADCPtoLChannelMap', 378, 'int16', tmp; 'nADCSamplingSeq', 410, 'int16', tmp; 'sADCChannelName',442, 'uchar', repmat(tmp,1,10); 'fADCProgrammableGain', 730, 'float', tmp; 'fInstrumentScaleFactor', 922, 'float', tmp; 'fInstrumentOffset', 986, 'float', tmp; 'fSignalGain', 1050, 'float', tmp; 'fSignalOffset', 1114, 'float', tmp; 'nTelegraphEnable',4512,'int16',tmp; 'fTelegraphAdditGain',4576,'float',tmp }; fields={'name','offs','numType','value'}; s=cell2struct(headPar,fields,2); numOfParams=size(s,1); clear tmp headPar; if ~exist(fn,'file'), error(['could not find file ' fn]); end if verbose, disp(['opening ' fn '..']); end [fid,messg]=fopen(fn,'r',machineF); if fid == -1,error(messg);end; % determine absolute file size fseek(fid,0,'eof'); fileSz=ftell(fid); fseek(fid,0,'bof'); % read all vital information in header % convert names in structure to variables and read value from header for g=1:numOfParams, if fseek(fid, s(g).offs,'bof')~=0, fclose(fid); error(['something went wrong locating ' s(g).name]); end; sz=length(s(g).value); eval(['[' s(g).name ',n]=fread(fid,sz,''' s(g).numType ''');']); if n~=sz, fclose(fid); error(['something went wrong reading value(s) for ' s(g).name]); end; end; if lActualAcqLength<nADCNumChannels, fclose(fid); error('less data points than sampled channels in file'); end % the numerical value of all recorded channels (numbers 0..15) recChIdx=nADCSamplingSeq(1:nADCNumChannels); % the corresponding indices into loaded data d recChInd=1:length(recChIdx); % the channel names, e.g. 'IN 8' recChNames=(reshape(char(sADCChannelName),10,16))'; recChNames=recChNames(recChIdx+1,:); chInd=[]; eflag=0; if ischar(channels) if strcmp(channels,'a') chInd=recChInd; else fclose(fid); error('input parameter ''channels'' must either be a cell array holding channel names or the single character ''a'' (=all channels)'); end else for i=1:length(channels) tmpChInd=strmatch(channels{i},recChNames,'exact'); if ~isempty(tmpChInd) chInd=[chInd tmpChInd]; else % set error flag to 1 eflag=1; end end; end if eflag fclose(fid); disp('**** available channels:'); disp(recChNames); disp(' '); disp('**** requested channels:'); disp(strvcat(channels)); error('at least one of the requested channels does not exist in data file (see above)'); end % fFileVersionNumber needs a fix - for whatever reason its value is sometimes % a little less than what it should be (e.g. 1.6499999xxxx instead of 1.65) fFileVersionNumber=.001*round(fFileVersionNumber*1000); % gain of telegraphed instruments, if any if fFileVersionNumber>=1.65 addGain=nTelegraphEnable.*fTelegraphAdditGain; addGain(addGain==0)=1; else addGain=ones(size(fTelegraphAdditGain)); end % tell me where the data start blockSz=512; switch nDataFormat case 0 dataSz=2; % bytes/point precision='int16'; case 1 dataSz=4; % bytes/point precision='float32'; otherwise fclose(fid); error('invalid number format'); end; headOffset=lDataSectionPtr*blockSz+nNumPointsIgnored*dataSz; % fADCSampleInterval is the TOTAL sampling interval si=fADCSampleInterval*nADCNumChannels; if ischar(sweeps) && sweeps=='a' nSweeps=lActualEpisodes; sweeps=1:lActualEpisodes; else nSweeps=length(sweeps); end; switch nOperationMode case 1 if verbose, disp('data were acquired in event-driven variable-length mode'); end warndlg('function abfload has not yet been thorougly tested for data in event-driven variable-length mode - please double-check that the data loaded is correct','Just a second, please'); if (lSynchArrayPtr<=0 || lSynchArraySize<=0), fclose(fid); error('internal variables ''lSynchArraynnn'' are zero or negative'); end switch fSynchTimeUnit case 0 % time information in synch array section is in terms of ticks synchArrTimeBase=1; otherwise % time information in synch array section is in terms of usec synchArrTimeBase=fSynchTimeUnit; end; % the byte offset at which the SynchArraySection starts lSynchArrayPtrByte=blockSz*lSynchArrayPtr; % before reading Synch Arr parameters check if file is big enough to hold them % 4 bytes/long, 2 values per episode (start and length) if lSynchArrayPtrByte+2*4*lSynchArraySize<fileSz, fclose(fid); error('file seems not to contain complete Synch Array Section'); end if fseek(fid,lSynchArrayPtrByte,'bof')~=0, fclose(fid); error('something went wrong positioning file pointer to Synch Array Section'); end [synchArr,n]=fread(fid,lSynchArraySize*2,'int32'); if n~=lSynchArraySize*2, fclose(fid); error('something went wrong reading synch array section'); end % make synchArr a lSynchArraySize x 2 matrix synchArr=permute(reshape(synchArr',2,lSynchArraySize),[2 1]); % the length of episodes in sample points segLengthInPts=synchArr(:,2)/synchArrTimeBase; % the starting ticks of episodes in sample points WITHIN THE DATA FILE segStartInPts=cumsum([0 (segLengthInPts(1:end-1))']*dataSz)+headOffset; % start time (synchArr(:,1)) has to be divided by nADCNumChannels to get true value % go to data portion if fseek(fid,headOffset,'bof')~=0, fclose(fid); error('something went wrong positioning file pointer (too few data points ?)'); end for i=1:nSweeps, % if selected sweeps are to be read, seek correct position if ~isequal(nSweeps,lActualEpisodes), fseek(fid,segStartInPts(sweeps(i)),'bof'); end; [tmpd,n]=fread(fid,segLengthInPts(sweeps(i)),precision); if n~=segLengthInPts(sweeps(i)), warning(['something went wrong reading episode ' int2str(sweeps(i)) ': ' segLengthInPts(sweeps(i)) ' points should have been read, ' int2str(n) ' points actually read']); end; dataPtsPerChan=n/nADCNumChannels; if rem(n,nADCNumChannels)>0, fclose(fid); error('number of data points in episode not OK'); end % separate channels.. tmpd=reshape(tmpd,nADCNumChannels,dataPtsPerChan); % retain only requested channels tmpd=tmpd(chInd,:); tmpd=tmpd'; % if data format is integer, scale appropriately; if it's float, tmpd is fine if ~nDataFormat for j=1:length(chInd), ch=recChIdx(chInd(j))+1; tmpd(:,j)=tmpd(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))... *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch); end; end % now place in cell array, an element consisting of one sweep with channels in columns d{i}=tmpd; end; case {2,5} if nOperationMode==2 if verbose disp('data were acquired in event-driven fixed-length mode'); end else if verbose disp('data were acquired in waveform fixed-length mode (clampex only)'); end end % determine first point and number of points to be read startPt=0; dataPts=lActualAcqLength; dataPtsPerChan=dataPts/nADCNumChannels; if rem(dataPts,nADCNumChannels)>0, fclose(fid); error('number of data points not OK'); end dataPtsPerChanPerSweep=dataPtsPerChan/lActualEpisodes; if rem(dataPtsPerChan,lActualEpisodes)>0 fclose(fid); error('number of data points not OK'); end dataPtsPerSweep=dataPtsPerChanPerSweep*nADCNumChannels; if fseek(fid,startPt*dataSz+headOffset,'bof')~=0 fclose(fid); error('something went wrong positioning file pointer (too few data points ?)'); end d=zeros(dataPtsPerChanPerSweep,length(chInd),nSweeps); % the starting ticks of episodes in sample points WITHIN THE DATA FILE selectedSegStartInPts=((sweeps-1)*dataPtsPerSweep)*dataSz+headOffset; for i=1:nSweeps, fseek(fid,selectedSegStartInPts(i),'bof'); [tmpd,n]=fread(fid,dataPtsPerSweep,precision); if n~=dataPtsPerSweep, fclose(fid); error(['something went wrong reading episode ' int2str(sweeps(i)) ': ' dataPtsPerSweep ' points should have been read, ' int2str(n) ' points actually read']); end; dataPtsPerChan=n/nADCNumChannels; if rem(n,nADCNumChannels)>0 fclose(fid); error('number of data points in episode not OK'); end; % separate channels.. tmpd=reshape(tmpd,nADCNumChannels,dataPtsPerChan); % retain only requested channels tmpd=tmpd(chInd,:); tmpd=tmpd'; % if data format is integer, scale appropriately; if it's float, d is fine if ~nDataFormat for j=1:length(chInd), ch=recChIdx(chInd(j))+1; tmpd(:,j)=tmpd(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))... *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch); end; end % now fill 3d array d(:,:,i)=tmpd; end; case 3 if verbose, disp('data were acquired in gap-free mode'); end % from start, stop, headOffset and fADCSampleInterval calculate first point to be read % and - unless stop is given as 'e' - number of points startPt=floor(1e6*start*(1/fADCSampleInterval)); % this corrects undesired shifts in the reading frame due to rounding errors in the previous calculation startPt=floor(startPt/nADCNumChannels)*nADCNumChannels; % if stop is a char array, it can only be 'e' at this point (other values would have % been caught above) if ischar(stop), dataPtsPerChan=lActualAcqLength/nADCNumChannels-floor(1e6*start/si); dataPts=dataPtsPerChan*nADCNumChannels; else dataPtsPerChan=floor(1e6*(stop-start)*(1/si)); dataPts=dataPtsPerChan*nADCNumChannels; if dataPts<=0 fclose(fid); error('start is larger than or equal to stop'); end end if rem(dataPts,nADCNumChannels)>0 fclose(fid); error('number of data points not OK'); end if fseek(fid,startPt*dataSz+headOffset,'bof')~=0, fclose(fid); error('something went wrong positioning file pointer (too few data points ?)'); end % *** decide on the most efficient way to read data: % (i) all (of one or several) channels requested: read, done % (ii) one (of several) channels requested: use the 'skip' feature of % fread % (iii) more than one but not all (of several) channels requested: % 'discontinuous' mode of reading data. Read a reasonable chunk of data % (all channels), separate channels, discard non-requested ones (if % any), place data in preallocated array, repeat until done. This is % faster than reading the data in one big lump, separating channels and % discarding the ones not requested if length(chInd)==1 && nADCNumChannels>1 % --- situation (ii) % jump to proper reading frame position in file if fseek(fid,(chInd-1)*dataSz,'cof')~=0 fclose(fid); error('something went wrong positioning file pointer (too few data points ?)'); end % read, skipping nADCNumChannels-1 data points after each read dataPtsPerChan=dataPts/nADCNumChannels; [d,n]=fread(fid,dataPtsPerChan,precision,dataSz*(nADCNumChannels-1)); if n~=dataPtsPerChan, fclose(fid); error(['something went wrong reading file (' int2str(dataPtsPerChan) ' points should have been read, ' int2str(n) ' points actually read']); end elseif length(chInd)/nADCNumChannels<1 % --- situation (iii) % prepare chunkwise upload: % preallocate d d=repmat(nan,dataPtsPerChan,length(chInd)); % the number of data points corresponding to the maximal chunk size, % rounded off such that from each channel the same number of points is % read (do not forget that each data point will by default be made a % double of 8 bytes, no matter what the original data format is) chunkPtsPerChan=floor(chunk*2^20/8/nADCNumChannels); chunkPts=chunkPtsPerChan*nADCNumChannels; % the number of those chunks.. nChunk=floor(dataPts/chunkPts); % ..and the remainder restPts=dataPts-nChunk*chunkPts; restPtsPerChan=restPts/nADCNumChannels; % chunkwise row indices into d dix=(1:chunkPtsPerChan:dataPtsPerChan)'; dix(:,2)=dix(:,1)+chunkPtsPerChan-1; dix(end,2)=dataPtsPerChan; if verbose && nChunk disp(['reading file in ' int2str(nChunk) ' chunks of ~' num2str(chunk) ' Mb']); end % do it: if no remainder exists loop through all rows of dix, % otherwise spare last row for the lines below (starting with % 'if restPts') for ci=1:size(dix,1)-(restPts>0) [tmpd,n]=fread(fid,chunkPts,precision); if n~=chunkPts fclose(fid); error(['something went wrong reading chunk #' int2str(ci) ' (' ... int2str(chunkPts) ' points should have been read, ' int2str(n) ' points actually read']); end % separate channels.. tmpd=reshape(tmpd,nADCNumChannels,chunkPtsPerChan); d(dix(ci,1):dix(ci,2),:)=tmpd(chInd,:)'; end % collect the rest [tmpd,n]=fread(fid,restPts,precision); if n~=restPts fclose(fid); error(['something went wrong reading last chunk (' ... int2str(restPts) ' points should have been read, ' int2str(n) ' points actually read']); end % separate channels.. tmpd=reshape(tmpd,nADCNumChannels,restPtsPerChan); d(dix(end,1):dix(end,2),:)=tmpd(chInd,:)'; else % --- situation (i) [d,n]=fread(fid,dataPts,precision); if n~=dataPts, fclose(fid); error(['something went wrong reading file (' int2str(dataPts) ' points should have been read, ' int2str(n) ' points actually read']); end; % separate channels.. d=reshape(d,nADCNumChannels,dataPtsPerChan); d=d'; end % if data format is integer, scale appropriately; if it's float, d is fine if ~nDataFormat for j=1:length(chInd), ch=recChIdx(chInd(j))+1; d(:,j)=d(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))... *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch); end; end otherwise disp('recording mode of data must be event-driven variable-length (1), event-driven fixed-length (2) or gap-free (3) -- returning empty matrix'); d=[]; si=[]; end; fclose(fid); %---------------------------------
Controller for abfload function
back to top%--------------------------------- function pvpmod(x)
% PVPMOD - evaluate parameter/value pairs % pvpmod(x) assigns the value x(i+1) to the parameter defined by the % string x(i) in the calling workspace. This is useful to evaluate % <varargin> contents in an mfile, e.g. to change default settings % of any variable initialized before pvpmod(x) is called. % % (c) U. Egert 1998 % this loop is assigns the parameter/value pairs in x to the calling % workspace. if ~isempty(x) for i = 1:2:size(x,2) assignin('caller', x{i}, x{i+1}); end; end; %---------------------------------
Filter waveform data
back to top%--------------------------------- function yfilt = myFilter(y)
%load data that need to be filtered %a zerophase filter will be used to process data %the zero phase filter parameters will filter out large noise data %like drift. %input % unfiltered signal "y" %output % filtered signal "yfilt" % Display start of filtering disp('Filtering Signal...'); disp(' '); % filter cell1 and assign to workspace windowSize = round(0.025*length(y))*2 + 1; %get an odd value for window, approx 5% wavefrom length tempfilt = filtfilt(ones(windowSize,1)/windowSize,1,y); %set zero phase filter parameters y = y - tempfilt; %remove any large noise fluxuations, like drift y = abs(-y); %make all values positive setzero = min(y); %get min value y = y-setzero; %%bring min value to zero yfilt = y; % Display end of filtering disp('...Filtering Complete'); disp(' '); %---------------------------------
Cropp waveform
back to top%--------------------------------- function [croppedSignal, limitsCrop, hfigOriginal, hfigCrop] = signalCrop(Signal)
% % [croppedSignal, limitsCrop] = signalCrop(Signal); % % Crop a signal by selecting points with mouse % % % % Input % % (Signal) - input waveform that needs to be cropped % % % % Output % % (croppedSignal) - cropped waveform % % (limitsCrop) - original indices or selected waveform to be cropped % % (hfigOriginal) - handle to original figure window % % (hfigCrop) - handle to cropped figure window mySignal = Signal; % initialize titles on figure hfigOriginal = figure('name','Original Signal'); plot(mySignal); %plot signal title({'CLICK TO CROP THE SIGNAL';... 'Use LEFT Click to add Vertical lmit';... 'Use RIGHT CLICK to add LAST limit'}); disp('CROP LIMITS'); disp('Left mouse button picks points.'); disp('Right mouse button picks last point.'); disp(' '); % use ginput to get user interaction. Continue until terminated with right % click xyCropp = []; %initialize clicked coordinates variable n = 0; %counts number of times clicked lastpoint = 1; %tracks left click, right click and middle click != 1 pointsY = axis; %assigns axis dimensions pointsY = pointsY(3:4); %keeps only mySignal dimensions while lastpoint == 1 [xi,yi,lastpoint] = ginput(1); %ginput calls input routine line([xi,xi],pointsY,'color','r') %uses last input and draws vertical line in red n = n+1; %increment clicks xyCropp(:,n) = [xi;yi]; %save click coordinates, use later to reference new cropped graph end % select subregion of signal using user input starty = round(min(xyCropp(1,:))); endy = round(max(xyCropp(1,:))); croppedSignal = mySignal(starty:endy); limitsCrop = [starty,endy]; % save the cropped limits for future use % plot the cropped signal hfigCrop = figure('name','Cropped Signal'); plot(croppedSignal); title('Cropped Signal'); %---------------------------------
Peak detection
back to top%--------------------------------- function [peakSet peakThreshold hfigPeak] = getPeaks(Signal)
% Interactive Peak Detection by Threshold % peakSet = getPeaks(Signal); % % Input % Signal - give a waveform signal, should be continuous % % Output % peakSet - gives the indices of all peaks above threshold level % peakThreshold - peak threshold value % hfigPeak - handle to peakSet figure window % % Interective Threshold method of peak detection. Left click allows to % refresh peak detection. Any click other than left ends peak search % and outputs peaks. % Search relies on user input threshold. All values above threshold are % isolated and clusted by checking for discontinuites with respect to % indices. The max value (y(x)) of each peak set of values is determined. % Noise, classified as signals that oscillate through the threshold on the % rising or falling legs of a peak, generates false positive results. % False positives are removed by comparing two lists, threshold vs % 50%threshold. Values that are found in both lists are preserved. % % E. Maravilla, 2014/03/05 mySignal = Signal; disp('PEAK THRESHOLD'); disp('Left mouse button sets threshold.') disp('Right mouse button sets threshold and terminates interaction.') disp(' '); xyThresh = []; %initialize threshold variable n = 0; %counts number of times clicked lastpoint = 1; %tracks left click, =1. Right click and middle click =2 & =3 % Initialize figure titles hfigPeak = figure('name','Cropped Signal for Peak Detection'); plot(mySignal); %plot signal title({'Click to ADD threshold value';... 'use LEFT CLICK to add a new threshold';... 'use RIGHT CLICK to add final threshold'}); pointsX = axis; %assigns axis dimensions pointsX = pointsX(1:2); %keeps only X dimensions % use ginput function to get user input. iterate until terminated with % right click while lastpoint == 1 [xi,yi,lastpoint] = ginput(1); %ginput calls input routine plot(mySignal); %redraw figure title({'Click to ADD threshold value';... 'use LEFT CLICK to add a new threshold';... 'use RIGHT CLICK to add final threshold'}); xlim(pointsX); %maintain x-axis limits line(pointsX,[yi,yi],'color','r'); %uses last input and draws vertical line in red n = n+1; %increment clicks xyThresh = [xi;yi]; %save click coordinates peaksThresh = yi; %threshold %%Geat max point from set of values above a threshold xThresh = find(mySignal > peaksThresh); %get index (w.r.t mySignal) above threshold value xThreshdiff = diff(xThresh); %get the neighbor distance between successive values that satisfied threshold screen %length is 1 less than xThresh length xCluster = find(xThreshdiff > 1); %get index that marks end of peak xCluster = xCluster+1; %shift, now index marks start of peaks %xCutoff = xThresh(xCluster); %get the original inices w.r.t ycropped of the cutoffs %yValues = mySignal(xThresh); %get y values of peaks %iterate through each set of peak values and find local max Apeak = []; for k = 1:length(xCluster)+1 %need one more iteration to find last peak if k == 1 %case when finding first peak peakSet = xThresh(1:xCluster(k)-1); [tempY,tempX] = max(mySignal(peakSet)); Apeak(1) = peakSet(tempX); elseif k == length(xCluster)+1 %case when finding last peak peakSet = xThresh(xCluster(k-1):end); %cannot access k when index = length+1 [tempY,tempX] = max(mySignal(peakSet)); Apeak(k) = peakSet(tempX); else %case when finding peaks between first and last peakSet = xThresh(xCluster(k-1):xCluster(k)-1); %never reached k becasue have separate case [tempY,tempX] = max(mySignal(peakSet)); Apeak(k) = peakSet(tempX); end end hold on; plot(Apeak,mySignal(Apeak),'go'); hold off; % Eliminate noise on rising/falling legs of peaks that cross thresh % level. This will eliminate false positive peaks if lastpoint ~= 1 plot(mySignal); %redraw figure title({'Click to ADD threshold value';... 'use LEFT CLICK to add a new threshold';... 'use RIGHT CLICK to add final threshold'}); xlim(pointsX); %maintina x-axis limits line(pointsX,[yi,yi],'color','r'); %uses last input and draws vertical line in red n = n+1; %increment clicks xyThresh = [xi;yi]; %save click coordinates peaksThresh = 0.25*yi; %--==Use lower thresh, Keep peaks that are present under both threshold conditions==-- %%Geat max point from set of values above a threshold xThresh = find(mySignal > peaksThresh); %get index (w.r.t mySignal) above threshold value xThreshdiff = diff(xThresh); %get the neighbor distance between successive values that satisfied threshold screen %length is 1 less than xThresh length xCluster = find(xThreshdiff > 1); %get index that marks end of peak xCluster = xCluster+1; %shift, now index marks start of peaks %xCutoff = xThresh(xCluster); %get the original inices w.r.t ycropped of the cutoffs %yValues = mySignal(xThresh); %get y values of peaks %iterate through each set of peak values and find local max Apeak2 = []; for k = 1:length(xCluster)+1 %need one more iteration to find last peak if k == 1 %case when finding first peak peakSet = xThresh(1:xCluster(k)-1); [tempY,tempX] = max(mySignal(peakSet)); Apeak2(1) = peakSet(tempX); elseif k == length(xCluster)+1 %case when finding last peak peakSet = xThresh(xCluster(k-1):end); %cannot access k when index = length+1 [tempY,tempX] = max(mySignal(peakSet)); Apeak2(k) = peakSet(tempX); else %case when finding peaks between first and last peakSet = xThresh(xCluster(k-1):xCluster(k)-1); %never reached k becasue have separate case [tempY,tempX] = max(mySignal(peakSet)); Apeak2(k) = peakSet(tempX); end end %%compare both peak lists, comparing x-coordinate, if found twice then keep [comboApeak,choose1] = max([length(Apeak),length(Apeak2)]); %find longest list switch choose1 case 1 temp = zeros(1,length(Apeak)); temp(1:length(Apeak2)) = Apeak2; Apeak2 = temp; matchApeak = []; countPeaks = 1; for k = 1:length(Apeak) checkPeak = find(Apeak == Apeak2(k)); %will give index of match or empty matrix if ~isempty(checkPeak) matchApeak(countPeaks) = checkPeak; countPeaks = countPeaks+1; else continue end end Apeak = Apeak(matchApeak); case 2 temp = zeros(1,length(Apeak2)); temp(1:length(Apeak)) = Apeak; Apeak = temp; matchApeak = []; countPeaks = 1; for k = 1:length(Apeak2) checkPeak = find(Apeak2 == Apeak(k)); if ~isempty(checkPeak) matchApeak(countPeaks) = checkPeak; countPeaks = countPeaks+1; else continue end end Apeak = Apeak2(matchApeak); end hold on; plot(Apeak,mySignal(Apeak),'go'); hold off; end end % close gcf peakThreshold = yi; peakSet = Apeak; %---------------------------------
Event Detection
back to top%--------------------------------- function [myEvents, myEventsFlag, hfigEvents] = getEvents(peakSet,Signal)
%Identify, Plot, and Flag Events as "Good" and "Bad" % % [myEvents, myEventsFlag] = getEvents(peakSet,Signal) % % -=INPUT=- % peakSet - index of list of all peaks w.r.t. Signal % Signal - data signal % % -=OUTPUT=- % myEvents - [#x3] list of #events with [start,peak,end] of event % myEventsFlag - [#x1] array of 1's and 0's indentifying "good" & "bad" events, 1=bad % hfigEvents = handle to events figure window y3 = Signal; %reassign variable peakcount = 1; %initialize peakcount, will be used to keep count when peak from peakSet is skipped for i=1:length(peakSet); % initial conditions for event % need to find a & c, the start and end index of event peaktemp = peakSet(i); %get index of peak w.r.t event % get index of peak b = peaktemp; % Left leg of peak if peaktemp < 101 % skip peak if the left of waveform from peak is less than 101 points continue end % compute a running average of slope values for j = 1:100 %iterate an arbitrary number of times and collect average slope values (span 5 points) ytemp = y3(peaktemp-5 +1 -(j-1):peaktemp -(j-1)); %-5+1 = span 5 to the left, -(j-1) = shift left by amount indicated counter slopetemp = diff(ytemp); %the 1st order discrete derivative using a curve that spans 5 points slopemean(j) = mean(slopetemp); %average of slope values end % find the point where the average slope falls below a given % of min slope value mycutoff = 0.0001; %arbitrary value for threshold cutoff slopethresh = []; %reset slpethresh while isempty(slopethresh) slopethresh = find(slopemean < mycutoff *min(slopemean)); mycutoff = mycutoff + 0.001; %increment threshold until successful if mycutoff > 0.5 break %stop loop if threshold value exceeds 50% end end slopethresh = slopethresh(1); %extract first instance along entire event that satisfied cutoff criteria % get start index of event a = peaktemp-2 -(slopethresh-1); % Right leg of peak if length(y3)-peaktemp < 101 % skip peak if the right of waveform from peak is less than 101 points continue end for j = 1:100 ytemp = y3(peaktemp +(j-1):peaktemp+5 -1 +(j-1)); %+5-1 = span 5, +(j-1) = shift right with counter slopetemp = diff(ytemp); slopemean(j) = mean(slopetemp); %average of slope values end % find the point where the average slope falls below % of min slope value windowSize = 13; Aveslopemean = filter(ones(1,windowSize)/windowSize,1,slopemean); %compute average of slopemean mycutoff = 0.0001; %establish threshold criteria slopethresh = []; %reset slopethresh while isempty(slopethresh) slopethresh = find(Aveslopemean > mycutoff *min(Aveslopemean)); mycutoff = mycutoff + 0.001; %increment threshold until successful if mycutoff > 0.5 break %stop loop if threshold value exceeds 50% end end slopethresh = slopethresh(1); %extract first instance along entire event that satisfied cutoff criteria % get end index of event c = peaktemp+2 +(slopethresh-1); % create an array storing the start, peak, and end of event peakEvents(peakcount,1:3) = [a, b, c]; %store curves to array peakcount = peakcount+1; %increment counter for event end myEvents = peakEvents; %copy to variable used by function % Plot the Signal with "bad" events highlighted in red hfigEvents = figure; plot(y3,'k'); hold on; for k = 1:length(peakEvents) a = peakEvents(k,1); c = peakEvents(k,3); plot(a:c,y3(a:c),'r'); end % Get flagged events myEventsFlag = flagEvents(myEvents,y3); %call flagEvents function %---------------------------------
Flag Bad Events
back to top%--------------------------------- function myEventsFlag = flagEvents(myEvents,Signal)
%Flag Bad Events % myEventsFilter = filterEvents(myEvents,Signal) % % Input % myEvents - output from myEvents = getEvents % Signal - waveform % % Output % myEventsFilter - array of 1 or 0 indicating "bad" events % % Modify vector shape % Combine the events into one variable for k = 1:length(myEvents) E{k,1} = Signal(myEvents(k,1):myEvents(k,3)); %y-values E{k,2} = myEvents(k,1):myEvents(k,3); %x-indices end % Split the event(s) in half, at peak value for k = 1:length(myEvents) Esplit{k,1} = Signal(myEvents(k,1):myEvents(k,2)); %left leg Esplit{k,2} = Signal(myEvents(k,2):myEvents(k,3)); %right leg end % First Filter Condition % Get the min y-value of each half and the mean y-value of the two combined for k = 1:length(Esplit) ypeakval = max(E{k,1}); %get peak value ytemp1 = Esplit{k,1}; %make temp waveform of rising leg [min1temp(k) min1idx(k)] = min(ytemp1); %get min y-value of risign leg and index ytemp2 = Esplit{k,2}; %make temp waveform of falling leg [min2temp(k) min2idx(k)] = min(ytemp2); %get min y-value of falling leg and index mean12temp(k) = mean([min1temp(k), min2temp(k)]); %get mean of both waveform mins yAmp = ypeakval - mean12temp(k); %get relative amplitude yshift = abs(min1temp(k) - min2temp(k)); %get difference between start and end conditions % Establish a flag that characterizes/identiies problematic event waveforms if yshift > 0.2 * yAmp %characteristic of "bad" waveforms myflag1(k) = 1; %bad waveform else myflag1(k) = 0; %good waveform end if length(ytemp2) < 10 myflag1(k) = 1; %bad waveform end end % Second Filter Condition % Establish 2nd criteria for identifying problematic events mystd = std(mean12temp); %standard deviation meanall = mean(mean12temp); %mean value of all events under consideration % 2nd criteria identifies false positives, removes events that have a mean 1 standard deviation away from group mean for k = 1:length(Esplit) if mean12temp(k) > meanall + 2*mystd %2 standard deviations from mean is general rule for exclusion myflag2(k) = 1; %bad waveform else myflag2(k) = 0; %good waveform end end % Combine both flags logically. If flag1==flag2, pass flag1 to flag. If % flag1 ~= flag2, pass 1 to flag. for k = 1:length(Esplit) if myflag1(k) == myflag2(k) myflag(k) = myflag1(k); elseif myflag1(k) ~= myflag2(k) myflag(k) = 1; end end % Third Filter Condition % Isolate events that were not filtered % use myflag to remove filtered events myEventsFlagKeep = find(myflag == 0); if length(myEventsFlagKeep) >= 1 for k = 1:length(myEventsFlagKeep) myEventsFilt{k,1} = E{myEventsFlagKeep(k),1}; myEventsFilt{k,2} = E{myEventsFlagKeep(k),2}; end % Establish routine that captures bias selection of "bad" events by using % the amplitude and the length of the evnent to establish programmable % criteria that recreates bias user selection with 90% accuracy. for k = 1:length(myEventsFilt) a = myEventsFilt{k,1}(1); b = max(myEventsFilt{k,1}); c = myEventsFilt{k,1}(end); ab = b-a; bc = b-c; EventAmp(k) = mean([ab bc]); EventLength(k) = myEventsFilt{k,2}(end) - myEventsFilt{k,2}(1); end EventAmpAve = mean(EventAmp); EventLengthAve = mean(EventLength); EventAmpComplex = EventAmp/EventAmpAve; EventLengthComplex = EventLength/EventLengthAve; EventComplex = EventAmpComplex.*EventLengthComplex; EventComplexAve = mean(EventComplex); % Loop through all events kept after filtering and filter again based on % the complex value determined from amplitude and length of event. Flag % events that fall below 20% of the EventComplexAve value for k = 1:length(myEventsFilt) if EventComplex(k) < 0.20*EventComplexAve myflag(myEventsFlagKeep(k)) = 1; %assign to index of myFlag end end end % % Plottting Events and identify the "bad" events % for k = 1:length(Esplit) % % ytemp1 = Esplit{k,1}; %make temp waveform of rising leg % ytemp2 = Esplit{k,2}; %make temp waveform of falling leg % % plot rising and falling legs in same figure, but in different axis % figure; % subplot(1,2,1); % if myflag(k) == 1 % plot(Esplit{k,1},'r','LineWidth',3); %paint red if "bad" % else % plot(Esplit{k,1},'g','LineWidth',3); %pain green if "good" % end % hold on; plot(min1idx(k),ytemp1(min1idx(k)),'ok','MarkerSize',8,'LineWidth',3.5); %mark mean of rising leg % plot(1:length(ytemp1),repmat(mean12temp(k),length(ytemp1),1),'--k','LineWidth',3.5); %mark mean of all event legs % subplot(1,2,2); % if myflag(k) == 1 % plot(Esplit{k,2},'r','LineWidth',3); %paint red if "bad" % else % plot(Esplit{k,2},'g','LineWidth',3); %pain green if "good" % end % hold on; plot(min2idx(k),ytemp2(min2idx(k)),'ok','MarkerSize',8,'LineWidth',3.5); %mark mean of rising leg % plot(1:length(ytemp2),repmat(mean12temp(k),length(ytemp2),1),'k--','LineWidth',3.5); %mark mean of all event legs % % end myEventsFlag = myflag'; %---------------------------------
Generate Descriptive Stats
back to top%--------------------------------- function [EventDataStats hfigStats] = getEventStats(myEvents,Signal,myEventsFlag,si)
% EventDataStats = getEventStats(myEvents,Signal,myEventsFlag); % get activaton and decay descriptives % % Input % myEvents = row array of peak waveform limits (start,peak,end) or % (start,end) % % Signal = waveform collecting stats from % % myEventsFlag = array of 1 & 0 indentifying "bad" events, 1 == bad % si = sample interval in seconds % % Output % EventDataStats - % EventDataStats(k,1) = Amptemp; %amplitude % EventDataStats(k,2) = 1/beta1temp; %activation time constant % EventDataStats(k,3) = 1/beta2temp; %decay time constant % EventDataStats(k,4) = AUCtemp; %area under curve % EventDataStats(k,5) = Rsquare1; %goodness of fit for activation % EventDataStats(k,6) = Rsquare2; %goodness of fit for decay % hfigStats - array of figure handles to fit curve figure windows % % external functions called % myfity1 = established to fit an exponential decay curve % used for transpose on activation, but should be modified % % Remove "bad" events from applying fit curves myEventstemp = myEvents; clear myEvents count = 1; %setup counter for indexing into new array for k = 1:length(myEventstemp) if myEventsFlag(k) == 1 continue %do not increment index if skipping a "bad" event else myEvents(count,:) = myEventstemp(k,:); count = count+1; %increment counter if copying "good" event end end % Display information about "good" and "bad" events disp('Analyzing...'); %show in command line disp(' '); disp(['Total Events Found = ', num2str(length(myEventstemp))]); %show total # of events disp(['"Good" Events Found = ', num2str(length(myEvents))]); %show # of "good" events disp(['"Bad" Events Found = ', num2str(length(myEventstemp) - length(myEvents))]); %show #of "Bad" events disp(' '); % Convert myEvents to a triple data array if is double set array if size(myEvents,2) == 2 %convert to triple set (start,peak,end) myEventsTemp = myEvents; clear myEvents for k = 1:length(myEventsTemp) a = myEventsTemp(k,1); c = myEventsTemp(k,2); [unused,bindex] = max(Signal(a:c)); b = a + bindex - 1; %get the peak with respect to start of event, subtract 1 b/c array indexing starts at 1, not 0 myEvents(k,1:3) = [a,c,b]; end end hfigs = []; %initialize array % Loop through each event and apply myfity, external fit curve function for k = 1:length(myEvents) %waveforms have 3 indices, column vector of successive events disp(['Fitting Event ', num2str(k)]); a = myEvents(k,1); b = myEvents(k,2); c = myEvents(k,3); ab = Signal(a:b); bc = Signal(b:c); ac = Signal(a:c); % compute the factor to make si = 0.0001 resF = si/0.0001; % resample data using resample function [yRising yDecay yResampled siNew] = myResample_simple_extend(Signal,a,c,si,resF); ytemp = ab;%activation curve ytemp = yRising; %resampled activation curve mintemp = min(ytemp); ytemp = ytemp - mintemp; %translate curve to origin %xtemp = 0:(length(ytemp)-1); %xtemp = xtemp*si; % plot if fitting an even event makefig1 = 'n'; % variable for makefig decision if makefig1 == 'y' hfigtemp1 = figure('name',['Rising Edge_Event_', num2str(k)]); %enter title as figure name hfigs(k,1) = hfigtemp1; end %[mtemp, stemp] = myfity1(flipud(ytemp),makefig1,si); %curve fit model & stats for activation [mtemp, stemp] = myfity2(ytemp,makefig1,siNew,' '); %curve fit model & stats for resampled activation curve beta1temp = mtemp.beta; %activation 1/tau Rsquare1 = stemp.rsquare; %R square error of fit % eval(['EPC' num2str(k) 'beta1 = beta1temp;']) % eval(['EPC' num2str(k) 'y1 = ytemp;']) % eval(['EPC' num2str(k) 'x1 = xtemp;']) % eval(['EPC' num2str(k) 'y1curve = mtemp;'])%curve fit model for activation % eval(['EPC' num2str(k) 'y1stat = stemp;'])%curve fit stats for activation ytemp = bc; %decay curve ytemp = yDecay; %resampled decay curve mintemp = min(ytemp); ytemp = ytemp - mintemp; %translate curve to origin %xtemp = 0:(length(ytemp)-1); %xtemp = xtemp*si; makefig2 = 'n'; % variable for makefig decision if makefig2 == 'y' hfigtemp2 = figure('name',['Falling Edge_Event_', num2str(k)]); %enter title as figure name hfigs(k,2) = hfigtemp2; end %[mtemp, stemp] = myfity1(ytemp,makefig2,si); %curve fit model & stats for decay [mtemp, stemp] = myfity1(ytemp,makefig2,siNew,' '); %curve fit model & stats for resampled decay curve beta2temp = mtemp.beta; %decay 1/tau Rsquare2 = stemp.rsquare; %R square error of fit % eval(['EPC' num2str(k) 'beta2 = beta2temp;']) % eval(['EPC' num2str(k) 'y2 = ytemp;']) % eval(['EPC' num2str(k) 'x2 = xtemp;']) % eval(['EPC' num2str(k) 'y1curve = mtemp;']) %curve fit model for decay % eval(['EPC' num2str(k) 'y1stat = stemp;']) %curve fit stats for decay [yRising yDecay yResampled siNew] = myResample_simple(Signal,a,c,si,resF); ytemp = ac; %whole waveform ytemp = yResampled; mintemp = min(ytemp); ytemp = ytemp - mintemp; %translate to origin for AUC determination xtemp = 0:(length(ytemp)-1); %xtemp = xtemp*si; xtemp = xtemp*siNew; % eval(['EPC' num2str(k) 'y3 = ytemp;']) % eval(['EPC' num2str(k) 'x3 = xtemp;']) AUCtemp = trapz(xtemp,ytemp); %compute area under curve % eval(['EPC' num2str(k) 'AUC = AUCtemp;']) % base1temp = mean (Signal(a-3:a)); % base2temp = mean (Signal(c:c+3)); % basetemp = mean([base1temp,base2temp]); % peaktemp = Signal(b); % Amptemp = abs(peaktemp - basetemp); %abs amplitude Amptemp = max(ytemp); % eval(['EPC' num2str(k) 'Amp = Amptemp;']) EventDataStats(k,1) = Amptemp; EventDataStats(k,2) = 1/beta1temp; EventDataStats(k,3) = 1/beta2temp; EventDataStats(k,4) = AUCtemp; EventDataStats(k,5) = Rsquare1; EventDataStats(k,6) = Rsquare2; %eval(['EPC' num2str(k) 'y1 = c1( a( (1+(k-1)*3) ):a( (2+(k-1)*3) ) );']) % eval( 'ytemp = eval([ ''min ( eval([''''EPC'''' num2str(k) ''''y1'''']) )'' ])' ); end hfigStats = hfigs; % Dispaly end of analysis disp(' '); disp('...Analysis Complete'); %---------------------------------
Exponential decay fit
back to top%--------------------------------- function [Curvefit, goodness] = myfity1(y1,makeplot,si,unitsY) % [Curvefit, goodness] = myfity1(y1,makeplot,si,yUnit) % Curve Fit % % Input % y1 - curve to be fitted must be exponential decay curve % makeplot - enter plot option PLOT (y), DONT PLOT (n) % si - sample interval in seconds (si) % unitsY - units of y-axis % % Output % Curvefit - curve fit to data % goodness - goodness of fit least square error % Y = y1(:); % make column vector % % Extend curve plateu value % a = ones(1, 10+length(Y) ); %make array length Y + 10 % a = a*Y(end); %make array equal to platue value % Y = [Y a]; %concatenate, extend platue x = 0:(length(Y)-1); %start x-axis from zero x = x*si; %convert to seconds x = x(:); %make column vector %x = x+(5*FrameTreat); %shift the curve to match time stamp %myc = round(min(Y)); %myalpha = round(max(Y)-myc); myalpha = max(Y); %temp = 0.63*myalpha+myc; %find value 63% temp = 0.63*myalpha; %find value 63% temp = abs(Y-temp); %subtract value 66% of max from curve [~, tempx] = min(temp); %get index of min difference = index of 63%myalpha+myc mybeta = 1/(x(tempx)); %r = 5*FrameTreat; %create right shift known constant (seconds) % % For condition where intercept is passed into fitoptions % s = fitoptions('Mehtod','NonlinearLeastSquares','Robust','LAR',... % 'Lower',[myalpha-.01,0,myc-.01],... % 'Upper',[myalpha+.01,inf,myc+.01],... % 'Startpoint',[myalpha,mybeta,myc]); %set the fitting parameters s = fitoptions('Mehtod','NonlinearLeastSquares','Robust','LAR',... 'Lower',[myalpha-1,0],... 'Upper',[myalpha+1,inf],... 'Startpoint',[myalpha,mybeta]); %set the fitting parameters % f = fittype( @(alpha,beta,c,x) alpha*(1 - exp(-beta *(x-r))) +c,... % 'coefficients',{'alpha','beta','c'},'options',s ); %set the curve fit function % f = fittype( @(alpha,beta,x) alpha*(1 - exp(-beta *(x))),... % 'coefficients',{'alpha','beta'},'options',s ); %set the curve fit function f = fittype( @(alpha,beta,x) alpha*(exp(-beta *(x))),... 'coefficients',{'alpha','beta'},'options',s ); %set the curve fit function [Curvefit, goodness] = fit(x,Y,f); %fit data to an eponential rise curve % Plot the original waveform, fit, and fit stats in figure % make plot of curve fit if requested if strncmpi(makeplot,'y',1) %case insensitive comparison of first character %figure; plot(Curvefit,'k-'),hold on %create a plot of the curve fit and original data plot(x, Y,'k.','markersize',12) %plot data points legend('Fit Curve', 'Data','Location','SouthEast') %place legend xlabel('Time(seconds)','FontSize',20) %x-axis label ylabel(['Amplitude (', unitsY, ')'], 'FontSize',20) %y-axis label Txticks = get(gca,'XTick'); %x-axis tick array Pxticks = round(length(get(gca,'XTick'))/4); %1/4 length of x-axis tick array xticks = Txticks(Pxticks); %1/4 length value of x-axis ticks %tempytick = abs(x-xticks); %x-axis data - x-axis 1/4 length tick value %[~, tempytick] = min(tempytick); %find position of min difference tempytick = 1; % make strings of equations using tex syntax strn1 = '$y={\alpha}\exp\left({-\beta}{x}\right)$'; strn2 = sprintf('$y={%0.4f}\\exp\\left({-%0.g}{x}\\right)$',Curvefit.alpha, Curvefit.beta); strn3 = sprintf('$R^2={%0.4f}$',abs(goodness.rsquare)); strn4 = sprintf('$\\tau=\\frac{1}{\\beta}={%0.4g}\\left( seconds \\right)$',(1/Curvefit.beta)); % % For plotting w.r.t. origin text('Interpreter','latex','string',... {strn1;strn2;strn3;strn4},... 'Position',[xticks,Y(tempytick)],... 'VerticalAlignment','Cap','FontSize',16) end
end %---------------------------------
Expoential rise fit
back to top%--------------------------------- function [Curvefit, goodness] = myfity2(y1,makeplot,si,unitsY) % [Curvefit, goodness] = myfity2(y1,makeplot,si,yUnit) % Curve Fit % % Input % y1 - curve to be fitted must be exponential rising curve % makeplot - enter plot option PLOT (y), DONT PLOT (n) % si - sample interval in seconds (si) % unitsY - units of y-axis % % Output % Curvefit - curve fit to data % goodness - goodness of fit least square error % Y = y1(:);%curve, make column vector %Y data index starts at 1 % Translate Y to origin, eliminating constant terms minY = min(Y); Y = Y - minY; % Create time vector x = 0:(length(Y)-1); %start x-axis from zero x = x*si; %convert to seconds x = x(:); %make column vector %x = x+(5*FrameTreat); %shift the curve to match time stamp %myc = round(min(Y)); %myalpha = round(max(Y)-myc); myalpha = max(Y); %temp = 0.63*myalpha+myc; %find value 63% temp = 0.63*myalpha; %find value 63% temp = abs(Y-temp); %subtract value 66% of max from curve [~, tempx] = min(temp); %get index of min difference = index of 63%myalpha+myc mybeta = 1/(x(tempx)); %r = 5*FrameTreat; %create right shift known constant (seconds) % s = fitoptions('Mehtod','NonlinearLeastSquares','Robust','LAR',... % 'Lower',[myalpha-.01,0,myc-.01],... % 'Upper',[myalpha+.01,inf,myc+.01],... % 'Startpoint',[myalpha,mybeta,myc]); %set the fitting parameters s = fitoptions('Mehtod','NonlinearLeastSquares','Robust','LAR',... 'Lower',[myalpha-1,0],... 'Upper',[myalpha+1,inf],... 'Startpoint',[myalpha,mybeta]); %set the fitting parameters % f = fittype( @(alpha,beta,c,x) alpha*(1 - exp(-beta *(x-r))) +c,... % 'coefficients',{'alpha','beta','c'},'options',s ); %set the curve fit function % f = fittype( @(alpha,beta,x) alpha*(exp(-beta *(x))),... % 'coefficients',{'alpha','beta'},'options',s ); %set the curve fit function f = fittype( @(alpha,beta,x) alpha*(1 - exp(-beta *(x))),... 'coefficients',{'alpha','beta'},'options',s ); %set the curve fit function [Curvefit, goodness] = fit(x,Y,f); %fit data to an eponential rise curve % Plot the original waveform, fit, and fit stats in figure % make plot of curve fit if requested if strncmpi(makeplot,'y',1) %case insensitive comparison of first character %figure; figure, plot(Curvefit,'k-'),hold on %create a plot of the curve fit and original data plot(x, Y,'k.','markersize',12) %plot data points legend('Fit Curve', 'Data','Location','SouthEast') %place legend xlabel('Time(seconds)','FontSize',20) %x-axis label ylabel(['Amplitude (', unitsY, ')'], 'FontSize',20) %y-axis label Txticks = get(gca,'XTick'); %x-axis tick array Pxticks = round(length(get(gca,'XTick'))/4); %1/4 length of x-axis tick array xticks = Txticks(Pxticks); %1/4 length value of x-axis ticks tempytick = abs(x-xticks); %x-axis data - x-axis 1/4 length tick value [~, tempytick] = min(tempytick); %find position of min difference % make strings of equations using tex syntax strn1 = '$y={\alpha}\left(1-\exp({-\beta}{x})\right)$'; strn2 = sprintf('$y={%0.4f}\\left(1-\\exp({-%0.g}{x})\\right)$',Curvefit.alpha, Curvefit.beta); strn3 = sprintf('$R^2={%0.4f}$',abs(goodness.rsquare)); strn4 = sprintf('$\\tau=\\frac{1}{\\beta}={%0.4g}\\left( seconds \\right)$',(1/Curvefit.beta)); % % For plotting w.r.t. origin text('Interpreter','latex','string',... {strn1;strn2;strn3;strn4},... 'Position',[xticks,Y(tempytick)],... 'VerticalAlignment','Cap','FontSize',16) end end %---------------------------------
Resample data version1
back to top%--------------------------------- function [yRising yDecay yResampled siNew] = myResample_simple(Signal,idxStart,idxEnd,si,resF)
% [yRising yDecay yResampled siNew] = myResample(Signal,idxStart,idxEnd,si) % implement resampling by interpolation % resample data to increase sample rate by a factor of 10 % example: old sample rate = 500 samples per second % new sample rate = 500*10 = 5000 samples per second % % Output % yRising = new y-values of length [((n-1)*10)+1], where n = old # of samples % yDecay = '' % yResampled = [yRising yDecay] % siNew = sampling interval for resampled data % % Input % Signal = y-values of signal as vector % idxStart = start of event % idxEnd = end of event % si = sample interval % resF = resampling factor % 19-Apr-2014 13:30:42 % % START % % use signal start and end limits to extract event y = Signal(idxStart:idxEnd); mintemp = min(y); y = y - mintemp; %translate curve to origin % setup and perform interpolation x = 1:length(y); %make x index vector of length y xi = 1:(1/resF):length(x); %make interpolation vector yi = interp1(x,y,xi,'spline'); %use spline interpolation method to resample data % get new sample interval time siNew = si/resF; %compute new sample interval based on resF resampling % % NEXT % % after resampling data find x-index for both rising and falling legs of % curve that point to the two y-values that are closest to 75%max(y-values) % divide curve into halfs and find index of y-value that has the minimum % square error away from 75%max(y-value) % % split curve into rising and decaying parts [yMaxval yMaxidx] = max(yi); if yMaxidx > length(yi)-20 %redo if max value is found at end of curve [yMaxval yMaxidx] = max(yi(1:end-20)); end PercMax = 0.75*yMaxval; %calculate 75% of amplitude aStart = 1; bMid = yMaxidx; cEnd = length(yi); % crop left leg from min y-value to end yLeft = yi(aStart:bMid); [notused yRisingtempidx] = min(yLeft); yRisingtemp = yLeft(yRisingtempidx:end); % % find left leg index whith y-value closest to 0.75*amp % yDiff = (yRisingtemp-PercMax).^2; % [notused yLeftidx] = min(yDiff); % eliminate starting tail from rising curve, using max(d/dx) & max(d/dx^2) [~,dx1] = max(diff(yRisingtemp)); [~,dx2] = max(diff(diff(yRisingtemp))); if dx2 < dx1 %&& abs(dx1-dx2) < 0.5*length(yRisingtemp) dxpos = dx2; else dxpos = dx1; end % Reset if selection failed by creeping too far up if yRisingtemp(dxpos)-min(yRisingtemp) > 0.3*max( (yRisingtemp-min(yRisingtemp)) ) dxpos = 1; end % build rising curve yRising = yRisingtemp(dxpos:end); % crop right leg from start to min y-value yRight = yi(bMid:cEnd); [notused yDecaytempidx] = min(yRight); yDecaytemp = yRight(1:yDecaytempidx); % find right leg index whith y-value closest to 0.75*amp yDiff = (yDecaytemp-PercMax).^2; [notused yRightidx] = min(yDiff); % build decay curve yDecay = yDecaytemp(yRightidx:end); % build whole curve using the left and right legs of event yResampled = yi((yRisingtempidx+dxpos-1):bMid+yDecaytempidx-1); %---------------------------------
Resample data version1
back to top%--------------------------------- function [yRising yDecay yResampled siNew] = myResample_simple(Signal,idxStart,idxEnd,si,resF)
% [yRising yDecay yResampled siNew] = myResample(Signal,idxStart,idxEnd,si) % implement resampling by interpolation % resample data to increase sample rate by a factor of 10 % example: old sample rate = 500 samples per second % new sample rate = 500*10 = 5000 samples per second % % Output % yRising = new y-values of length [((n-1)*10)+1], where n = old # of samples % yDecay = '' % yResampled = [yRising yDecay] % siNew = sampling interval for resampled data % % Input % Signal = y-values of signal as vector % idxStart = start of event % idxEnd = end of event % si = sample interval % resF = resampling factor % 19-Apr-2014 13:30:42, 12-May-2014 10:06:54 % % START % % use signal start and end limits to extract event y = Signal(idxStart:idxEnd); mintemp = min(y); y = y - mintemp; %translate curve to origin % setup and perform interpolation x = 1:length(y); %make x index vector of length y xi = 1:(1/resF):length(x); %make interpolation vector yi = interp1(x,y,xi,'spline'); %use spline interpolation method to resample data % get new sample interval time siNew = si/resF; %compute new sample interval based on resF resampling % % NEXT % % after resampling data find x-index for both rising and falling legs of % curve that point to the two y-values that are closest to 75%max(y-values) % divide curve into halfs and find index of y-value that has the minimum % square error away from 75%max(y-value) % % split curve into rising and decaying parts [yMaxval yMaxidx] = max(yi); if yMaxidx > length(yi)-20 %redo if max value is found at end of curve [yMaxval yMaxidx] = max(yi(1:end-20)); end PercMax = 0.75*yMaxval; %calculate 75% of amplitude aStart = 1; bMid = yMaxidx; cEnd = length(yi); % crop left leg from min y-value to end yLeft = yi(aStart:bMid); [notused yRisingtempidx] = min(yLeft); yRisingtemp = yLeft(yRisingtempidx:end); % % find left leg index whith y-value closest to 0.75*amp % yDiff = (yRisingtemp-PercMax).^2; % [notused yLeftidx] = min(yDiff); % eliminate starting tail from rising curve, using max(d/dx) & max(d/dx^2) [~,dx1] = max(diff(yRisingtemp)); [~,dx2] = max(diff(diff(yRisingtemp))); if dx2 < dx1 %&& abs(dx1-dx2) < 0.5*length(yRisingtemp) dxpos = dx2; else dxpos = dx1; end % Reset if selection failed by creeping too far up if yRisingtemp(dxpos)-min(yRisingtemp) > 0.3*max( (yRisingtemp-min(yRisingtemp)) ) dxpos = 1; end % build rising curve yRising = yRisingtemp(dxpos:end); % extend platue by 5% a = ones(1, round(0.05*length(yRising))+length(yRising) ); %make array length Y + 5%(length Y) a = a*yRising(end); %make array equal to platue value yRising = [yRising a]; %concatenate, extend platue % crop right leg from start to min y-value yRight = yi(bMid:cEnd); [notused yDecaytempidx] = min(yRight); yDecaytemp = yRight(1:yDecaytempidx); % find right leg index whith y-value closest to 0.75*amp yDiff = (yDecaytemp-PercMax).^2; [notused yRightidx] = min(yDiff); % build decay curve yDecay = yDecaytemp(yRightidx:end); % build whole curve using the left and right legs of event yResampled = yi((yRisingtempidx+dxpos-1):bMid+yDecaytempidx-1); %---------------------------------
Get relative data paths to data files
back to top
%---------------------------------
function rawDataPaths()
% ADD SPECIFIC FOLDERS TO SEARCH PATH CONTAINING RAW DATA, ABF FILES
% ...\EPCs\Group 1;
% ...\EPCs\Group 2;
% ...\MEPPs\Group 1;
% ...\MEPPs\Group 2;
% along with all subfolders
% FINAL WORKING FORM
% 07-May-2014 00:39:20
% ADD SPECIFIC FOLDERS TO SEARCH PATH CONTAINING RAW DATA, ABF FILES
mystring = pwd; %get full path of current directory
myPath_2 = regexp(mystring, '.+(?=\\.+\\.+$)','match'); %get parsed path name, two folders up
FolderSearch_2 = char(myPath_2); %get path to folder 2 level above where scripts is
% Generate character list for all paths in FolderSearch2
All_paths_2temp = genpath(FolderSearch_2); %
% Using the original fullfile output, Searching for
% ...\EPCs\Group 1;
% ...\EPCs\Group 2;
% ...\MEPPs\Group 1;
% ...\MEPPs\Group 2;
% Scan string with all full file paths
matchStr_2 = regexp(All_paths_2temp, '[A-Z]:\\[^;]+\\(EPCs|MEPPs)\\Group \d;', 'match');
% Replace the semicolon(;) with backslash(\)
matchStr_22 = regexprep(matchStr_2, ';', '\\');
% Add paths of folders
for k = 1:length(matchStr_22)
addpath(genpath(matchStr_22{k}));