bleeck@4: % function [pitchstrength, dominant_time] = analyse_timeinterval_profile(ti_profile, peaks, a_priori, fqp_fq)
bleeck@4: %
bleeck@4: %   To analyse the time interval profile of the auditory image
bleeck@4: %
bleeck@4: %   INPUT VALUES:
bleeck@4: %  
bleeck@4: %   RETURN VALUE:
bleeck@4: %
bleeck@4:  
bleeck@4: % (c) 2011, University of Southampton
bleeck@4: % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@4: % download of current version is on the soundsoftware site: 
bleeck@4: % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@4: % documentation and everything is on http://www.acousticscale.org
bleeck@4: 
bleeck@4: function result = analyse_timeinterval_profile(ti_profile,options)
bleeck@4: 
bleeck@4: if nargin < 2
bleeck@4: 	options=[];
bleeck@4: end
bleeck@4: 
bleeck@4: if isfield(options,'target_frequency')
bleeck@4: end
bleeck@4: 
bleeck@4: if ~isfield(options,'options')
bleeck@4: 	target_frequency=0;
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: % for debug resons
bleeck@4: plot_switch = 1;
bleeck@4: 
bleeck@4: 
bleeck@4: intsum = ti_profile;   
bleeck@4: intsum_vals = getdata(intsum);
bleeck@4: if plot_switch
bleeck@4:     minimum_time_interval=0;  % in ms
bleeck@4:     maximum_time_interval=100;
bleeck@4:     figure(654654)
bleeck@4: 	clf
bleeck@4:     tip = intsum_vals;
bleeck@4:     tip_x = bin2time(ti_profile, 1:length(tip));  % Get the times
bleeck@4:     tip_x = tip_x((tip_x>=(minimum_time_interval/1000)) & tip_x<=(maximum_time_interval/1000));  
bleeck@4:     tip = tip(time2bin(ti_profile, tip_x(1)):time2bin(ti_profile, tip_x(end)));
bleeck@4:     % tip_x is in ms. Change to Hz
bleeck@4:     tip_x = 1./tip_x;
bleeck@4:     plot(tip_x, tip, 'b');
bleeck@4:     set(gca,'XScale','log'); 
bleeck@4: 	set(gca,'Ylim',	[min(intsum)-10 max(intsum)+10]);
bleeck@4: 	set(gca,'Xlim',	[30 1500]);
bleeck@4:     hold on
bleeck@4: end
bleeck@4: if length(peaks)==0
bleeck@4:     % If there is no peak than we just plot the function
bleeck@4:     pitchstrength = 0;
bleeck@4:     found_frequency = 0;
bleeck@4:     return
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: sig=envelope(intsum)
bleeck@4: 
bleeck@4: % ++++++++++++ Part one: Peak finding ++++++++++
bleeck@4: % Calculate the envelope (curve of the peaks)
bleeck@4: 
bleeck@4: % sort peaks from low time interval to a heigh one
bleeck@4: % or from left to right
bleeck@4: peaks_lr = sortstruct(peaks,'x');  % 
bleeck@4: 
bleeck@4: if plot_switch
bleeck@4:     px=[];py=[];
bleeck@4:     for i=1:length(peaks)
bleeck@4:         px = [px tip_x(peaks{i}.x)];
bleeck@4:         py = [py tip(peaks{i}.x)];
bleeck@4:     end
bleeck@4:     plot(px,py,'kx');
bleeck@4: end
bleeck@4: 
bleeck@4: % Create envelope of peaks
bleeck@4: peak_envX = [];
bleeck@4: peak_envY = [];
bleeck@4: for i=1:length(peaks_lr)
bleeck@4:     peak_envX = [peak_envX peaks_lr{i}.x];
bleeck@4:     peak_envY = [peak_envY peaks_lr{i}.y];
bleeck@4: end
bleeck@4: 
bleeck@4: if plot_switch
bleeck@4:      plot(tip_x(peak_envX),  peak_envY, ':k');
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: % Find Maxima of the envelope 
bleeck@4: % create signal
bleeck@4: 
bleeck@4: sig=buildfrompoints(sig,xx,yy)
bleeck@4: 
bleeck@4: 
bleeck@4: peak_envsig = signal(length(peak_envX), 1);
bleeck@4: peak_envsig = setvalues(peak_envsig, peak_envY);
bleeck@4: params = 0; 
bleeck@4: peak_env_peaks = PeakPicker(peak_envsig, params);
bleeck@4: % sort peaks of the envelope from low time interval to a heigh one
bleeck@4: % or from left to right
bleeck@4: peaks_env_peaks_lr = sortstruct(peak_env_peaks,'x');  % 
bleeck@4: 
bleeck@4: if plot_switch
bleeck@4: 	for i=2:length(peaks_env_peaks_lr) % construct all maxima
bleeck@4: 		fre=peaks_env_peaks_lr{i}.t;
bleeck@4: 		plot(fre,gettimevalue(ti_profile,1/fre),'Marker','o','Markerfacecolor','g','MarkeredgeColor','g','MarkerSize',5);
bleeck@4: 	end
bleeck@4: end
bleeck@4: 
bleeck@4: % % Now make sure, that highest peak is not the first peak of envelope
bleeck@4: % peak_oi = peaks{1};
bleeck@4: % if (peak_oi.x == peak_envX(peaks_env_peaks_lr{1}.x))
bleeck@4: %     % The first Peak is the heighest -> take second highest of envelope
bleeck@4: %     if (length(peak_env_peaks)>=2)
bleeck@4: %          poix = peak_envX(peak_env_peaks{2}.x);
bleeck@4: %          poiy = peak_env_peaks{2}.y;
bleeck@4: %          for i=1:length(peaks)
bleeck@4: %              if poix==peaks{i}.x
bleeck@4: %                 peak_oi = peaks{i};  
bleeck@4: %              end
bleeck@4: %          end
bleeck@4: %      end
bleeck@4: % end
bleeck@4: 
bleeck@4: 
bleeck@4: % Stefans new method:
bleeck@4: % Take the one with the highest contrast
bleeck@4: if length(peaks_env_peaks_lr) > 1
bleeck@4: 	for i=2:length(peaks_env_peaks_lr) % construct all maxima
bleeck@4: 		maxpeak=peaks_env_peaks_lr{i}.y;
bleeck@4: 		minpeak1=peaks_env_peaks_lr{i}.left.y;
bleeck@4: 		minpeak2=peaks_env_peaks_lr{i}.right.y;
bleeck@4: 		minpeak=mean([minpeak1 minpeak2]);
bleeck@4: 		
bleeck@4: 		% the definition of pitch strength: Contrast
bleeck@4: 		% 	contrast(i)=maxpeak/(mean([minpeak1 minpeak2]));
bleeck@4: 		
bleeck@4: 		% the difference between both
bleeck@4: 		if maxpeak-minpeak > 0
bleeck@4: 			contrast(i)=(maxpeak-minpeak)/1000;
bleeck@4: 		else
bleeck@4: 			contrast(i)=0;
bleeck@4: 		end
bleeck@4: 		if plot_switch
bleeck@4: 			fre=1/peaks{i}.t;
bleeck@4: 			plot(fre,gettimevalue(ti_profile,1/fre),'Marker','o','Markerfacecolor','g','MarkeredgeColor','g','MarkerSize',5);
bleeck@4: 			text(fre,gettimevalue(ti_profile,1/fre)*1.1,sprintf('%2.2f',contrast(i)));
bleeck@4: 		end
bleeck@4: 	end
bleeck@4: else	% there is only one peak
bleeck@4: 	maxpeak=peaks{1}.y;
bleeck@4: 	minpeak1=peaks{1}.left.y;
bleeck@4: 	minpeak2=peaks{1}.right.y;
bleeck@4: 	minpeak=mean([minpeak1 minpeak2]);
bleeck@4: 	if maxpeak-minpeak > 0
bleeck@4: 		contrast(1)=(maxpeak-minpeak)/1000;
bleeck@4: 	else
bleeck@4: 		contrast(1)=0;
bleeck@4: 	end
bleeck@4: end
bleeck@4: 
bleeck@4: [maxcontrast,wheremax]=max(contrast);
bleeck@4: 
bleeck@4: peak_oi=peaks{wheremax};
bleeck@4: 
bleeck@4: pitchstrength= contrast(wheremax);
bleeck@4: found_frequency = 1/peak_oi.t;
bleeck@4: 
bleeck@4: if plot_switch
bleeck@4: 	fre=found_frequency;
bleeck@4: 	plot(fre,gettimevalue(ti_profile,1/fre),'Marker','o','Markerfacecolor','g','MarkeredgeColor','g','MarkerSize',pitchstrength*10);
bleeck@4: 	text(fre*1.1,gettimevalue(ti_profile,1/fre)+5, ['highest pitchstrength found at ' num2str(round(fre)) 'Hz: ' num2str(fround(pitchstrength,2)) ],'VerticalAlignment','bottom','HorizontalAlignment','center','color','g','Fontsize',12);
bleeck@4: end
bleeck@4: 
bleeck@4: return
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: % maxpeak=maxstruct(peaks,'y');
bleeck@4: 
bleeck@4: if length(peaks_env_peaks_lr)>1
bleeck@4:     for i=1:length(peaks)
bleeck@4:         % If second highes peak==peak with smallest time intervall -> take
bleeck@4:         % third highest !!
bleeck@4:         if peak_env_peaks{2}.x == peaks_env_peaks_lr{1}.x
bleeck@4:           poix = peak_envX(peak_env_peaks{3}.x);
bleeck@4:         else
bleeck@4:           poix = peak_envX(peak_env_peaks{2}.x);
bleeck@4:         end
bleeck@4:         if peaks{i}.x==poix
bleeck@4:             peak_oi = peaks{i};
bleeck@4:         end
bleeck@4:     end
bleeck@4: else
bleeck@4:     % pure sinusoid ???
bleeck@4:     dominant_time = 0
bleeck@4:     pitchstrength = 0.001;
bleeck@4:     return
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: % Stefan's method on HCL
bleeck@4: % Take second peak of envelope from short time intervals
bleeck@4: if length(peaks_env_peaks_lr)>1
bleeck@4:     for i=1:length(peaks)
bleeck@4:         % If second highes peak==peak with smallest time intervall -> take
bleeck@4:         % third highest !!
bleeck@4:         if peak_env_peaks{2}.x == peaks_env_peaks_lr{1}.x
bleeck@4:           poix = peak_envX(peak_env_peaks{3}.x);
bleeck@4:         else
bleeck@4:           poix = peak_envX(peak_env_peaks{2}.x);
bleeck@4:         end
bleeck@4:         if peaks{i}.x==poix
bleeck@4:             peak_oi = peaks{i};
bleeck@4:         end
bleeck@4:     end
bleeck@4: else
bleeck@4:     % pure sinusoid ???
bleeck@4:     dominant_time = 0
bleeck@4:     pitchstrength = 0.001;
bleeck@4:     return
bleeck@4: end
bleeck@4: 
bleeck@4: if plot_switch
bleeck@4:     plot(tip_x(peak_oi.x), peak_oi.y, 'ro');
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: % peak_oi contains the peak for quantification
bleeck@4: 
bleeck@4: 
bleeck@4: % ++++++++++++ Part two: Quantification ++++++++++
bleeck@4: 
bleeck@4: % **** First method
bleeck@4: % pitchstrength is mean of sum / mean of peaks
bleeck@4: % psum = 0;
bleeck@4: % for i = 1:length(peaks)
bleeck@4: %     psum = psum+peaks{i}.y;
bleeck@4: % end
bleeck@4: % psum = psum / length(peaks);
bleeck@4: % pitchstrength = psum/peaks{1}.y;
bleeck@4: % dominant_time = peaks{1}.x /getsr(ti_profile);
bleeck@4: % if plot_switch
bleeck@4: %     hold off
bleeck@4: %     plot(ti_profile);
bleeck@4: %     hold on;
bleeck@4: %     plot(peaks{1}.x, peaks{1}.y, 'ko');
bleeck@4: % end
bleeck@4: 
bleeck@4: % % ***** Second method
bleeck@4: % % Heigh of neighbour peak / highest peak
bleeck@4: % % First Peaks is the highest as pitchstrength of Peak Picker
bleeck@4: % % Find lower neighbour (with shorter time interval)
bleeck@4: % ioi=1;  % index of interest
bleeck@4: % dist = inf;
bleeck@4: % for i=1:length(peaks)
bleeck@4: %     d = peak_oi.x - peaks{i}.x;
bleeck@4: %     if ((d<dist)&(d>0))
bleeck@4: %         ioi = i;
bleeck@4: %         dist = d;
bleeck@4: %     end
bleeck@4: % end
bleeck@4: % pitchstrength = peaks{ioi}.y/ peak_oi.y;
bleeck@4: % dominant_time = peak_oi.x /getsr(ti_profile);
bleeck@4: % if plot_switch
bleeck@4: %     hold on;
bleeck@4: %     plot(peak_oi.x , peak_oi.y , 'ko');
bleeck@4: %     plot(peaks{ioi}.x, peaks{ioi}.y, 'go');
bleeck@4: % end
bleeck@4: 
bleeck@4: % % Third Method:
bleeck@4: % Peak to mean valley ration - good with HCL
bleeck@4: pitchstrength= mean([peak_oi.left.y peak_oi.right.y])/peak_oi.y;
bleeck@4: dominant_time = peak_oi.x /getsr(ti_profile);
bleeck@4: 
bleeck@4: % Adaptation
bleeck@4: %pitchstrength = 1-pitchstrength;
bleeck@4: 
bleeck@4: % +++++++++++++ Histogram of distances between two peaks ++++++++++
bleeck@4: 
bleeck@4: %
bleeck@4: max_thresh = 1; %0.98;
bleeck@4: % The linear TIP 
bleeck@4: 
bleeck@4: % Test if all other Peaks ar multiple of the highest peak
bleeck@4: 
bleeck@4: % Start with the first peak
bleeck@4: % Test how many peaks are mutiple of the first peak
bleeck@4: %---x---x---x---x---x---x
bleeck@4: % Take first peak of envelope.
bleeck@4: poix = peak_envX(peak_env_peaks{1}.x);
bleeck@4: for i=1:length(peaks_lr)
bleeck@4:     if peaks_lr {i}.x==poix
bleeck@4:         peak_first = peaks_lr{i};
bleeck@4:         first_i=i;
bleeck@4:     end
bleeck@4: end
bleeck@4: nomult = 0; % Number of multible
bleeck@4: max_delta = 0.1;   % 
bleeck@4: for i=first_i:length(peaks_lr)
bleeck@4:     f = peaks_lr{i}.t / peak_first.t;
bleeck@4:     if abs(round(f)-f)<max_delta
bleeck@4:         nomult=nomult+1;
bleeck@4:     end
bleeck@4: end
bleeck@4: nomult;
bleeck@4: 
bleeck@4: if plot_switch
bleeck@4:     figure(savecf)
bleeck@4: end
bleeck@4: 
bleeck@4: %---x---x---x---x---x---x
bleeck@4: 
bleeck@4: % is there is a priori information? 
bleeck@4: % if (nargin>2)
bleeck@4: %     dominant_time = mb(h_index) / getsr(ti_profile);
bleeck@4: %     pitchstrength = 0; 
bleeck@4: %     return
bleeck@4: % end
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: % ------------ subfunctions ---------------------
bleeck@4: 
bleeck@4: % turns a vector (row) upside down
bleeck@4: function y=upsidedown(x)
bleeck@4: y=[];
bleeck@4: for i=length(x):-1:1
bleeck@4:     y=[y x(i)];
bleeck@4: end