bleeck@3: % This external file is included as part of the 'aim-mat' distribution package
bleeck@3: % (c) 2011, University of Southampton
bleeck@3: % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@3: % download of current version is on the soundsoftware site: 
bleeck@3: % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@3: % documentation and everything is on http://www.acousticscale.org
tomwalters@0: 
tomwalters@0: function [lowpasssig,maxinfo,mininfo]=peak_picker(sig,options);
tomwalters@0: % sophisticated peak picker
tomwalters@0: % the signal is first lowpassfilterd with the frequency given in options
tomwalters@0: % options comes with:
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % frequency with witch the psth is filtered
tomwalters@0: % options.lowpassfrequency=500;
tomwalters@0: % search for peaks in the range from to:
tomwalters@0: % options.search_peak_start_search=0.001; % seach also before the offical latency
tomwalters@0: % options.search_peak_stop_search=0.025; % seach this time for peaks
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % every distinct peak must comply the following conditions:
tomwalters@0: % a certain activity during the duration of the maximum (a spike count per
tomwalters@0: % presentation)
tomwalters@0: if ~isfield(options,'min_count_per_peak')
tomwalters@0:     options.min_count_per_peak=-inf;% default: no certain activity: all peaks are significant
tomwalters@0: end
tomwalters@0: % a certain height of the maximum peak
tomwalters@0: if ~isfield(options,'min_height_of_heighest_peak')
tomwalters@0:     options.min_height_of_heighest_peak=0; % default: no certain height: all peaks are significant
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % return values are the filtered signal and information about the found
tomwalters@0: % peaks:
tomwalters@0: % maxinfo.values  % the hights of the found maxima
tomwalters@0: % maxinfo.times   % where are the found maxima
tomwalters@0: % maxinfo.contrast % the contrast is the relation between the maximum and the surrounding minima
tomwalters@0: % maxinfo.qvalue  % how wide the maximum is at half height.
tomwalters@0: % maxinfo.activity % sum of activity between the adjacent minima
tomwalters@0: 
tomwalters@0: % every found maxima (except the first and the last) is surrounded by two
tomwalters@0: % minima. Same is true for the found minimas
tomwalters@0: 
tomwalters@0: % do we want some grafical output (for debugging)
tomwalters@0: grafix=options.grafix;
tomwalters@0: % define the return values
tomwalters@0: maxinfo=[];
tomwalters@0: mininfo=[];
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %do the lowpassfiltering    
tomwalters@0: lowpasssig=lowpass(sig,options.lowpassfrequency);
tomwalters@0: [atmax,ahmax]=getlocalmaxima(lowpasssig);
tomwalters@0: [atmin,ahmin]=getlocalminima(lowpasssig);
tomwalters@0: 
tomwalters@0: % restrict to the requiered range
tomwalters@0: indeces1=find(atmax>options.search_peak_start_search);
tomwalters@0: indeces2=find(atmax<options.search_peak_stop_search);
tomwalters@0: indeces=intersect(indeces1,indeces2);
tomwalters@0: tmax=atmax(indeces);
tomwalters@0: hmax=ahmax(indeces);
tomwalters@0: indeces1=find(atmin>options.search_peak_start_search);
tomwalters@0: indeces2=find(atmin<options.search_peak_stop_search);
tomwalters@0: indeces=intersect(indeces1,indeces2);
tomwalters@0: tmin=atmin(indeces);
tomwalters@0: hmin=ahmin(indeces);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % make an iterated search through all maxima and decide which ones to keep
tomwalters@0: finished=0;
tomwalters@0: while ~finished
tomwalters@0:     % throw out exactly one maximum
tomwalters@0:     [tmaxnew,tminnew,hmaxnew,hminnew]=try_reduce_maxima(tmax,tmin,hmax,hmin,options);
tomwalters@0:     if length(tmaxnew)==length(tmax) || length(tmaxnew)==1 || length(tminnew)==1
tomwalters@0:         finished=1;
tomwalters@0:     end
tomwalters@0:     tmax=tmaxnew;    tmin=tminnew;    hmax=hmaxnew;    hmin=hminnew; % new values
tomwalters@0: %     if grafix
tomwalters@0: %         plotall(fignum,lowpasssig,tmax,tmin,hmax,hmin);
tomwalters@0: %         p=0;
tomwalters@0: %     end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if grafix
tomwalters@0:     fignum=figure;
tomwalters@0:     set(gcf,'Number','off');
tomwalters@0: %     set(gcf,'name',sprintf('peak picker (An:%s Un:%s Ex:%s) ',data.unitinfo.an_num,data.unitinfo.un_num,data.unitinfo.ex_num));
tomwalters@0:     set(gcf,'name',('peak picker'));
tomwalters@0:     plotall(fignum,lowpasssig,tmax,tmin,hmax,hmin);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % if grafix
tomwalters@0: %     plotall(lowpasssig,tmax,tmin,hmax,hmin);
tomwalters@0: % end
tomwalters@0: % now only significant values are left!
tomwalters@0: % put them in the return structure:
tomwalters@0: 
tomwalters@0: nr_max=length(tmax);
tomwalters@0: maxinfo.values=hmax;  % the hights of the found maxima
tomwalters@0: maxinfo.times=tmax;   % where are the found maxima
tomwalters@0: [highest_peak_height,highest_peak_index]=max(hmax);
tomwalters@0: 
tomwalters@0: for i=1:nr_max
tomwalters@0:     maxinfo.contrast(i)=getcontrast(i,tmax,tmin,hmax,hmin); % the contrast is the relation between the maximum and the surrounding minima
tomwalters@0:     maxinfo.qvalue(i)=getquality(i,tmax,tmin,hmax,hmin);  % how wide the maximum is at half height.
tomwalters@0:     maxinfo.activity(i)=getactivity(i,lowpasssig,tmax,tmin,hmax,hmin,options.latency);  % sum of activity between the adjacent minima
tomwalters@0: end
tomwalters@0: 
tomwalters@0: nr_min=length(tmin);
tomwalters@0: mininfo.values=hmin;  % the hights of the found maxima
tomwalters@0: mininfo.times=tmin;   % where are the found maxima
tomwalters@0: for i=1:nr_min
tomwalters@0:     mininfo.contrast(i)=getmincontrast(i,tmax,tmin,hmax,hmin); % the contrast is the relation between the maximum and the surrounding minima
tomwalters@0:     mininfo.qvalue(i)=getminquality(i,tmax,tmin,hmax,hmin);  % how wide the maximum is at half height.
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % now we have all peaks, find out which ones are significant for us
tomwalters@0: maxinfo.distinct_max=[];
tomwalters@0: count=1;
tomwalters@0: height_criterium=options.min_height_of_heighest_peak*highest_peak_height;
tomwalters@0: count_criterium=options.min_count_per_peak;
tomwalters@0: 
tomwalters@0: for i=1:nr_max
tomwalters@0:     % it must have a certain contrast
tomwalters@0:     if maxinfo.contrast(i)>=options.min_contrast_for_distinct_peak
tomwalters@0:         % and it must have a certain height
tomwalters@0:         if maxinfo.values(i)>=height_criterium
tomwalters@0:             if maxinfo.activity(i)>=count_criterium
tomwalters@0:                 maxinfo.distinct_max(count).contrast=maxinfo.contrast(i);
tomwalters@0:                 maxinfo.distinct_max(count).qvalue=maxinfo.qvalue(i);
tomwalters@0:                 maxinfo.distinct_max(count).activity=maxinfo.activity(i);
tomwalters@0:                 maxinfo.distinct_max(count).hmax=hmax(i);
tomwalters@0:                 maxinfo.distinct_max(count).tmax=tmax(i);
tomwalters@0:                 count=count+1;
tomwalters@0:             end
tomwalters@0:         end
tomwalters@0:     end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if grafix
tomwalters@0:     nr_all_max=length(atmax);
tomwalters@0:     for i=1:nr_all_max
tomwalters@0:         plot(atmax(i)*1000,ahmax(i),'o','markersize',2,'Markerfacecolor','r','Markeredgecolor','r');
tomwalters@0:     end
tomwalters@0:     nr_all_min=length(atmin);
tomwalters@0:     for i=1:nr_all_min
tomwalters@0:         plot(atmin(i)*1000,ahmin(i),'o','markersize',2,'Markerfacecolor','g','Markeredgecolor','g');
tomwalters@0:     end
tomwalters@0: 
tomwalters@0:     % plot a dot for every disticnt maximum found
tomwalters@0:     for i=1:nr_max
tomwalters@0:         testmaxpos=maxinfo.times(i);
tomwalters@0:         testmaxval=maxinfo.values(i);
tomwalters@0:         contrast=maxinfo.contrast(i);
tomwalters@0:         count=maxinfo.activity(i);
tomwalters@0:         plot(testmaxpos*1000,testmaxval,'o','markersize',8,'Markerfacecolor','r','Markeredgecolor','r');
tomwalters@0:     end
tomwalters@0:     % plot a dot for every minimum found
tomwalters@0:     for i=1:length(maxinfo.distinct_max);
tomwalters@0:         testmaxpos=maxinfo.distinct_max(i).tmax;
tomwalters@0:         testmaxval=maxinfo.distinct_max(i).hmax;
tomwalters@0:         contrast=maxinfo.distinct_max(i).contrast;
tomwalters@0:         count=maxinfo.distinct_max(i).activity;
tomwalters@0:         text(testmaxpos*1000+1,testmaxval,sprintf('contr %3.3f',contrast),'fontsize',6,'ver','bottom');
tomwalters@0:         text(testmaxpos*1000+1,testmaxval,sprintf('count %3.3f',count),'fontsize',6,'ver','top');
tomwalters@0:     end
tomwalters@0:     for i=1:nr_min
tomwalters@0:         testminpos=mininfo.times(i);
tomwalters@0:         testminval=mininfo.values(i);
tomwalters@0:         plot(testminpos*1000,testminval,'o','markersize',4,'Markerfacecolor','g','Markeredgecolor','g');
tomwalters@0:     end
tomwalters@0: end
tomwalters@0: return
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: function plotall(fignum,lowpasssig,tmax,tmin,hmax,hmin);
tomwalters@0: figure(fignum);
tomwalters@0: clf
tomwalters@0: hold on
tomwalters@0: fill(lowpasssig,'b');
tomwalters@0: set(gca,'xlim',[0,60]);
tomwalters@0: for i=1:length(tmax)
tomwalters@0:     plot(tmax(i)*1000,hmax(i),'Marker','o','MarkerSize',6,'MarkerFaceColor','r','MarkerEdgeColor','r')
tomwalters@0: end
tomwalters@0: for i=1:length(tmin)
tomwalters@0: %     plot(time2bin(lowpasssig,tmin(i)),hmin(i),'Marker','o','MarkerSize',6,'MarkerFaceColor','g','MarkerEdgeColor','g')
tomwalters@0: end
tomwalters@0: xlabel('time (ms)');
tomwalters@0: ylabel('spikes / sweep / bin');
tomwalters@0: title(' PSTH plus found maxima');
tomwalters@0: % movegui(nfig,'northwest');
tomwalters@0: return
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: function [tmaxnew,tminnew,hmaxnew,hminnew]=try_reduce_maxima(tmax,tmin,hmax,hmin,options)
tomwalters@0: % throw out the peak with the smallest contrast 
tomwalters@0: % if removed, the minimum to the right of it is removed too
tomwalters@0: 
tomwalters@0: nr_max=length(tmax);
tomwalters@0: 
tomwalters@0: % find the maximum with the neighbour that is closest in height and throw
tomwalters@0: % it out!
tomwalters@0: for i=1:nr_max
tomwalters@0:     contrastr(i)=getrightcontrast(i,tmax,tmin,hmax,hmin);
tomwalters@0: end
tomwalters@0: for i=1:nr_max
tomwalters@0:     contrastl(i)=getleftcontrast(i,tmax,tmin,hmax,hmin);
tomwalters@0: end
tomwalters@0: [mincontrastr,minconindexr]=min(abs(contrastr));    % thats the smallest right 
tomwalters@0: [mincontrastl,minconindexl]=min(abs(contrastl));    % thats the smallest left
tomwalters@0: 
tomwalters@0: if mincontrastr<mincontrastl
tomwalters@0:     side='right';
tomwalters@0:     [mincontrast,minconindex]=min(contrastr);
tomwalters@0: else
tomwalters@0:     side='left';
tomwalters@0:     [mincontrast,minconindex]=min(contrastl);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: testmaxpos=tmax(minconindex);
tomwalters@0: % if mincontrast<options.min_contrast_for_distinct_peak
tomwalters@0: if mincontrast<options.min_contrast_for_peak
tomwalters@0:     switch side
tomwalters@0:         case 'right'
tomwalters@0:             [testminpos,testminval,indexmin]=getminimumrightof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0:         case 'left'
tomwalters@0:             [testminpos,testminval,indexmin]=getminimumleftof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0:     end
tomwalters@0:     
tomwalters@0:     % remove the minimum right of the maximum
tomwalters@0:     tmaxnew=mysetdiff(tmax,tmax(minconindex));
tomwalters@0:     hmaxnew=mysetdiff(hmax,hmax(minconindex));
tomwalters@0:     if indexmin>0 && indexmin <length(tmin)
tomwalters@0:         tminnew=mysetdiff(tmin,tmin(indexmin));
tomwalters@0:         hminnew=mysetdiff(hmin,hmin(indexmin));
tomwalters@0:     else
tomwalters@0:         tminnew=tmin;
tomwalters@0:         hminnew=hmin;
tomwalters@0:     end
tomwalters@0:     return % only remove one!
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % if still here, then nothing happend
tomwalters@0: tmaxnew=tmax;    tminnew=tmin;     hmaxnew=hmax;    hminnew=hmin;
tomwalters@0: 
tomwalters@0: return
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function ret=mysetdiff(values,singleval)
tomwalters@0: % the same as setdiff, but without the sorting
tomwalters@0: ret=[];
tomwalters@0: for i=1:length(values)
tomwalters@0:     if singleval~=values(i)
tomwalters@0:         ret=[ret values(i)];
tomwalters@0:     end
tomwalters@0: end
tomwalters@0: return
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function contrast=getcontrast(peaknumber,tmax,tmin,hmax,hmin)
tomwalters@0: % return the contrast of the peak with the number
tomwalters@0: testmaxpos=tmax(peaknumber);
tomwalters@0: testmaxval=hmax(peaknumber);
tomwalters@0: [leftminpos,leftminval]=getminimumleftof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: [rightminpos,rightminval]=getminimumrightof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: if isempty(leftminpos) && isempty(rightminpos)  % if both are empty, its difficult
tomwalters@0:     contrast=0;
tomwalters@0: elseif isempty(leftminpos)  % if only the left is empty, take the right instead
tomwalters@0:     contrast=(testmaxval-rightminval)/(testmaxval+rightminval);
tomwalters@0: elseif isempty(rightminpos)  % if only the right is empty, take the left instead
tomwalters@0:     contrast=(testmaxval-leftminval)/(testmaxval+leftminval);
tomwalters@0: else
tomwalters@0:     mean_min_val=mean([leftminval rightminval]);
tomwalters@0:     contrast=(testmaxval-mean_min_val)/(testmaxval+mean_min_val);
tomwalters@0: end
tomwalters@0: return
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function contrast=getrightcontrast(peaknumber,tmax,tmin,hmax,hmin)
tomwalters@0: % return the contrast with the right minimum of the peak with the number
tomwalters@0: testmaxpos=tmax(peaknumber);
tomwalters@0: testmaxval=hmax(peaknumber);
tomwalters@0: % [leftminpos,leftminval]=getminimumleftof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: [rightminpos,rightminval]=getminimumrightof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: % if isempty(leftminpos) && isempty(rightminpos)  % if both are empty, its difficult
tomwalters@0: %     contrast=0;
tomwalters@0: % elseif isempty(leftminpos)  % if only the left is empty, take the right instead
tomwalters@0: %     contrast=(testmaxval-rightminval)/(testmaxval+rightminval);
tomwalters@0: if isempty(rightminpos)  % if only the right is empty, take the left instead
tomwalters@0:     contrast=inf;
tomwalters@0: else
tomwalters@0:     %     mean_min_val=mean([leftminval rightminval]);
tomwalters@0:     contrast=(testmaxval-rightminval)/(testmaxval+rightminval);
tomwalters@0: end
tomwalters@0: return
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function contrast=getleftcontrast(peaknumber,tmax,tmin,hmax,hmin)
tomwalters@0: % return the contrast with the right minimum of the peak with the number
tomwalters@0: testmaxpos=tmax(peaknumber);
tomwalters@0: testmaxval=hmax(peaknumber);
tomwalters@0: [leftminpos,leftminval]=getminimumleftof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: % [rightminpos,rightminval]=getminimumrightof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: % if isempty(leftminpos) && isempty(rightminpos)  % if both are empty, its difficult
tomwalters@0: %     contrast=0;
tomwalters@0: % elseif isempty(leftminpos)  % if only the left is empty, take the right instead
tomwalters@0: %     contrast=(testmaxval-rightminval)/(testmaxval+rightminval);
tomwalters@0: if isempty(leftminpos)  % if only the right is empty, take the left instead
tomwalters@0:     contrast=inf;
tomwalters@0: else
tomwalters@0:     %     mean_min_val=mean([leftminval rightminval]);
tomwalters@0:     contrast=(testmaxval-leftminval)/(testmaxval+leftminval);
tomwalters@0: end
tomwalters@0: return
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function contrast=getmincontrast(troughnumber,tmax,tmin,hmax,hmin)
tomwalters@0: % return the contrast of the peak with the number
tomwalters@0: testminpos=tmin(troughnumber);
tomwalters@0: testminval=hmin(troughnumber);
tomwalters@0: [leftmaxpos,leftmaxval]=getmaximumleftof(testminpos,tmax,tmin,hmax,hmin);
tomwalters@0: [rightmaxpos,rightmaxval]=getmaximumrightof(testminpos,tmax,tmin,hmax,hmin);
tomwalters@0: if isempty(leftmaxpos) && isempty(rightmaxpos)  % if both are empty, its difficult
tomwalters@0:     contrast=0;
tomwalters@0: elseif isempty(leftmaxpos)  % if only the left is empty, take the right instead
tomwalters@0:     contrast=(testminval-rightmaxval)/(testminval+rightmaxval);
tomwalters@0: elseif isempty(rightmaxpos)  % if only the right is empty, take the left instead
tomwalters@0:     contrast=(testminval-leftmaxval)/(testminval+leftmaxval);
tomwalters@0: else
tomwalters@0:     mean_max_val=mean([leftmaxval rightmaxval]);
tomwalters@0:     contrast=(testminval-mean_max_val)/(testminval+mean_max_val);
tomwalters@0: end
tomwalters@0: return
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function peakheight=getpeakheight(peaknumber,tmax,tmin,hmax,hmin)
tomwalters@0: testmaxpos=tmax(peaknumber);
tomwalters@0: testmaxval=hmax(peaknumber);
tomwalters@0: [leftminpos,leftminval]=getminimumleftof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: [rightminpos,rightminval]=getminimumrightof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: if isempty(leftminpos) && isempty(rightminpos)  % if both are empty, its difficult
tomwalters@0:     peakheight=testmaxval;
tomwalters@0: elseif isempty(leftminpos)  % if only the left is empty, take the right instead
tomwalters@0:     peakheight=testmaxval-rightminval;
tomwalters@0: elseif isempty(rightminpos)  % if only the right is empty, take the left instead
tomwalters@0:     peakheight=testmaxval-leftminval;
tomwalters@0: else
tomwalters@0:     mean_min_val=mean([leftminval rightminval]);
tomwalters@0:     peakheight=testmaxval-mean_min_val;
tomwalters@0: end
tomwalters@0: return
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function troughheight=gettroughheight(troughnumber,tmax,tmin,hmax,hmin)
tomwalters@0: testminpos=tmin(troughnumber);
tomwalters@0: testminval=hmin(troughnumber);
tomwalters@0: [leftmaxpos,leftmaxval]=getminimumleftof(testminpos,tmax,tmin,hmax,hmin);
tomwalters@0: [rightmaxpos,rightmaxval]=getminimumrightof(testminpos,tmax,tmin,hmax,hmin);
tomwalters@0: if isempty(leftmaxpos) && isempty(rightmaxpos)  % if both are empty, its difficult
tomwalters@0:     troughheight=abs(testminval);
tomwalters@0: elseif isempty(leftmaxpos)  % if only the left is empty, take the right instead
tomwalters@0:     troughheight=abs(testminval-rightmaxval);
tomwalters@0: elseif isempty(rightmaxpos)  % if only the right is empty, take the left instead
tomwalters@0:     troughheight=abs(testminval-leftmaxval);
tomwalters@0: else
tomwalters@0:     mean_min_val=mean([leftmaxval rightmaxval]);
tomwalters@0:     troughheight=abs(testminval-mean_min_val);
tomwalters@0: end
tomwalters@0: return
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: function quality=getquality(peaknumber,tmax,tmin,hmax,hmin);  % how wide the maximum is at the height of its surrounding minima
tomwalters@0: testmaxpos=tmax(peaknumber);
tomwalters@0: testmaxval=hmax(peaknumber);
tomwalters@0: [leftminpos,leftminval]=getminimumleftof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: [rightminpos,rightminval]=getminimumrightof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: if isempty(leftminpos) && isempty(rightminpos)  % if both are empty, its difficult
tomwalters@0:     quality=0;
tomwalters@0: elseif isempty(leftminpos)  % if only the left is empty, take the right instead
tomwalters@0:     quality=0;
tomwalters@0: elseif isempty(rightminpos)  % if only the right is empty, take the left instead
tomwalters@0:     quality=0;
tomwalters@0: else
tomwalters@0:     diff_min_val=abs(leftminpos-rightminpos);
tomwalters@0:     quality=testmaxval/diff_min_val;
tomwalters@0: end
tomwalters@0: return
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function quality=getminquality(troughknumber,tmax,tmin,hmax,hmin);  % how wide the maximum is at the height of its surrounding maxima
tomwalters@0: testminpos=tmin(troughknumber);
tomwalters@0: testminval=hmin(troughknumber);
tomwalters@0: [leftmaxpos,leftmaxval]=getmaximumleftof(testminpos,tmax,tmin,hmax,hmin);
tomwalters@0: [rightmaxpos,rightmaxval]=getmaximumrightof(testminpos,tmax,tmin,hmax,hmin);
tomwalters@0: if isempty(leftmaxpos) && isempty(rightmaxpos)  % if both are empty, its difficult
tomwalters@0:     quality=0;
tomwalters@0: elseif isempty(leftmaxpos)  % if only the left is empty, take the right instead
tomwalters@0:     quality=0;
tomwalters@0: elseif isempty(rightmaxpos)  % if only the right is empty, take the left instead
tomwalters@0:     quality=0;
tomwalters@0: else
tomwalters@0:     diff_min_val=abs(leftmaxval-rightmaxval);
tomwalters@0:     quality=testminval/diff_min_val;
tomwalters@0: end
tomwalters@0: return
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: function activity=getactivity(peaknumber,sig,tmax,tmin,hmax,hmin,latency);  % sum of activity between the adjacent minima
tomwalters@0: testmaxpos=tmax(peaknumber);
tomwalters@0: testmaxval=hmax(peaknumber);
tomwalters@0: [leftminpos,leftminval]=getminimumleftof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: [rightminpos,rightminval]=getminimumrightof(testmaxpos,tmax,tmin,hmax,hmin);
tomwalters@0: if isempty(leftminpos) && isempty(rightminpos)  % if both are empty, its difficult
tomwalters@0:     activity=0;
tomwalters@0: elseif isempty(leftminpos)  % if only the left is empty, take the right instead
tomwalters@0:     activity=0; % take the activity from the start of the signal instead
tomwalters@0:     sr=1000/getsr(sig);
tomwalters@0:     activity=sum(sig,latency,rightminpos)/sr;
tomwalters@0: 
tomwalters@0: elseif isempty(rightminpos)  % if only the right is empty, take the left instead
tomwalters@0:     activity=0;
tomwalters@0: else
tomwalters@0:     sr=1000/getsr(sig);
tomwalters@0:     activity=sum(sig,leftminpos,rightminpos)/sr;
tomwalters@0: end
tomwalters@0: return