bleeck@4: % function [pitchstrength, dominant_time] = find_pitches(profile, , 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 = find_pitches(profile,options)
bleeck@4: 
bleeck@4: plot_switch =0;
bleeck@4: 
bleeck@4: if nargin < 2
bleeck@4: 	options=[];
bleeck@4: end
bleeck@4: % 
bleeck@4: 
bleeck@4: % peaks must be higher then this of the maximum to be recognised
bleeck@4: ps_threshold=0.1;	
bleeck@4: heighttowidthminimum=0.3;	% height of peak, where the width is measured
bleeck@4: octave_variability=0.35;	% in percent, when the octave is acepted as such
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: % preliminary return values.
bleeck@4: result.final_pitchstrength=0;
bleeck@4: result.final_dominant_frequency=0;
bleeck@4: result.final_dominant_time=0;
bleeck@4: result.smoothed_signal=[];
bleeck@4: result.peaks = [];
bleeck@4: 
bleeck@4: % change the scaling to logarithm, before doing anything else:
bleeck@4: log_profile=logsigx(profile,0.001,0.035);
bleeck@4: 
bleeck@4: 
bleeck@4: peak_found=0;
bleeck@4: current_lowpass_frequency=500;
bleeck@4: 
bleeck@4: while ~peak_found
bleeck@4: 	
bleeck@4: 	smooth_sig=lowpass(log_profile,current_lowpass_frequency);
bleeck@4: 	envpeaks = IPeakPicker(smooth_sig,0.01);
bleeck@4: 	
bleeck@4: 	if isempty(envpeaks) % only, when there is no signal
bleeck@4: 		return
bleeck@4: 	end
bleeck@4: 	
bleeck@4: 	% finde den peak mit der höchsten Pitschstrength
bleeck@4: 	for i=1:length(envpeaks) % construct all maxima
bleeck@4: 		where=bin2time(smooth_sig,envpeaks{i}.x);
bleeck@4: 		[dummy,height,breit,where_widths]=qvalue(smooth_sig,where,heighttowidthminimum);
bleeck@4: 		width=time2bin(smooth_sig,breit);
bleeck@4: 		if width>0
bleeck@4: 			envpeaks{i}.pitchstrength=height/width;
bleeck@4: 		else
bleeck@4: 			envpeaks{i}.pitchstrength=0;
bleeck@4: 		end
bleeck@4: 		envpeaks{i}.width=width;
bleeck@4: 		envpeaks{i}.where_widths=where_widths;
bleeck@4: 		envpeaks{i}.peak_base_height_y=height*(1-heighttowidthminimum);
bleeck@4: 	end
bleeck@4: 	
bleeck@4: 	% sort all for the highest first!
bleeck@4: 	% 	envpeaks=sortstruct(envpeaks,'pitchstrength');
bleeck@4: 	envpeaks=sortstruct(envpeaks,'y');
bleeck@4: 	
bleeck@4: 	if plot_switch
bleeck@4: 		figure(2134)
bleeck@4: 		clf
bleeck@4: 		hold on
bleeck@4: 		plot(log_profile,'b')
bleeck@4: 		plot(smooth_sig,'g')
bleeck@4: 		for i=1:length(envpeaks)
bleeck@4: 			time=envpeaks{i}.t;
bleeck@4: 			x=envpeaks{i}.x;
bleeck@4: 			y=envpeaks{i}.y;
bleeck@4: 			plot(x,y,'Marker','o','Markerfacecolor','g','MarkeredgeColor','g','MarkerSize',2);
bleeck@4: 			time_left=envpeaks{i}.left.t;
bleeck@4: 			x_left=envpeaks{i}.left.x;
bleeck@4: 			y_left=envpeaks{i}.left.y;
bleeck@4: 			time_right=envpeaks{i}.right.t;
bleeck@4: 			x_right=envpeaks{i}.right.x;
bleeck@4: 			y_right=envpeaks{i}.right.y;
bleeck@4: 			plot(x_left,y_left,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',2);
bleeck@4: 			plot(x_right,y_right,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',2);
bleeck@4: 		end
bleeck@4: 		current_time=options.current_time*1000;
bleeck@4: 		y=get(gca,'ylim');
bleeck@4: 		y=y(2);
bleeck@4: 		text(700,y,sprintf('%3.0f ms',current_time),'verticalalignment','top');
bleeck@4: 		
bleeck@4: 		for i=1:length(envpeaks)
bleeck@4: 			pitchstrength=envpeaks{i}.pitchstrength;
bleeck@4: 			if i <4
bleeck@4: 				line([envpeaks{i}.x envpeaks{i}.x],[envpeaks{i}.y envpeaks{i}.y*(1-heighttowidthminimum)],'Color','r');
bleeck@4: 				line([envpeaks{i}.where_widths(1) envpeaks{i}.where_widths(2)],[envpeaks{i}.y*(1-heighttowidthminimum) envpeaks{i}.y*(1-heighttowidthminimum)],'Color','r');
bleeck@4: 				x=envpeaks{i}.x;
bleeck@4: 				fre=envpeaks{i}.fre;
bleeck@4: 				y=envpeaks{i}.y;
bleeck@4: 				text(x,y*1.05,sprintf('%3.0fHz: %2.2f',fre,pitchstrength),'HorizontalAlignment','center');
bleeck@4: 			end
bleeck@4: 		end
bleeck@4: 	end
bleeck@4: 	
bleeck@4: 	% kucke, ob sich das erste Maxima überlappt. Falls ja, nochmal
bleeck@4: 	% berechnen mit niedriger Lowpass frequenz
bleeck@4: 	nr_peaks=length(envpeaks);
bleeck@4: 	if 	nr_peaks==1
bleeck@4: 		peak_found=1;
bleeck@4: 		continue
bleeck@4: 	end
bleeck@4: 	xleft=envpeaks{1}.where_widths(1);
bleeck@4: 	xright=envpeaks{1}.where_widths(2);
bleeck@4: 	for i=2:nr_peaks %
bleeck@4: 		xnull=envpeaks{i}.x;
bleeck@4: 		if xnull < xleft  && xnull > xright
bleeck@4: 			peak_found=0;
bleeck@4: 			current_lowpass_frequency=current_lowpass_frequency/2;
bleeck@4: 			if current_lowpass_frequency<62.5
bleeck@4: 				return
bleeck@4: 			end
bleeck@4: 			break
bleeck@4: 		end
bleeck@4: 		% wenn noch hier, dann ist alles ok
bleeck@4: 		peak_found=1;
bleeck@4: 	end
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: % reduce the peaks to the relevant ones:
bleeck@4: % through out all with pitchstrength smaller then threshold
bleeck@4: count=1;
bleeck@4: min_ps=ps_threshold*envpeaks{1}.pitchstrength;
bleeck@4: for i=1:length(envpeaks)
bleeck@4: 	if envpeaks{i}.pitchstrength>min_ps
bleeck@4: 		rpeaks{count}=envpeaks{i};
bleeck@4: 		count=count+1;
bleeck@4: 	end
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: % final result for the full set
bleeck@4: result.final_pitchstrength=rpeaks{1}.pitchstrength;
bleeck@4: result.final_dominant_frequency=rpeaks{1}.fre;
bleeck@4: result.final_dominant_time=rpeaks{1}.t;
bleeck@4: result.smoothed_signal=smooth_sig;
bleeck@4: result.peaks=rpeaks;
bleeck@4: % return
bleeck@4: 
bleeck@4: % Neuberechnung der Pitchstrength für peaks, die das Kriterium der 
bleeck@4: % Höhe zu Breite nicht erfüllen
bleeck@4: for i=1:length(rpeaks) % construct all maxima
bleeck@4: 	where=bin2time(smooth_sig,rpeaks{i}.x);
bleeck@4: 	[dummy,height,breit,where_widths]=qvalue(smooth_sig,where,heighttowidthminimum);
bleeck@4: 	width=time2bin(smooth_sig,breit);
bleeck@4: 	
bleeck@4: 	xleft=where_widths(2);
bleeck@4: 	xright=where_widths(1);
bleeck@4: 	left_peak=rpeaks{i}.left;
bleeck@4: 	right_peak=rpeaks{i}.right;
bleeck@4: 	
bleeck@4: 	
bleeck@4: 	xnull=rpeaks{i}.x;
bleeck@4: 	if xleft<left_peak.x || xright>right_peak.x
bleeck@4: 		t_xnull=bin2time(smooth_sig,xnull);
bleeck@4: 		[val,height,breit,widthvals,base_peak_y]=qvalue2(smooth_sig,t_xnull);
bleeck@4: 		width=time2bin(smooth_sig,breit);
bleeck@4: 		if width>0
bleeck@4: 			rpeaks{i}.pitchstrength=height/width;
bleeck@4: 		else
bleeck@4: 			rpeaks{i}.pitchstrength=0;
bleeck@4: 		end
bleeck@4: 		rpeaks{i}.where_widths=time2bin(smooth_sig,widthvals);
bleeck@4: 		rpeaks{i}.width=width;
bleeck@4: 		rpeaks{i}.peak_base_height_y=base_peak_y;
bleeck@4: 	end
bleeck@4: end
bleeck@4: 
bleeck@4: % sort again all for the highest first!
bleeck@4: % rpeaks=sortstruct(rpeaks,'pitchstrength');
bleeck@4: 
bleeck@4: % return
bleeck@4: %%%%%%  look for octave relationships
bleeck@4: % 
bleeck@4: % 
bleeck@4: for i=1:length(rpeaks) % construct all maxima
bleeck@4: 	% look, if the octave of the pitch is also present.
bleeck@4: 	fre=rpeaks{i}.fre;
bleeck@4: 	has_octave=0;
bleeck@4: 	for j=1:length(rpeaks)
bleeck@4: 		oct_fre=rpeaks{j}.fre;
bleeck@4: 		% is the octave there?
bleeck@4: 		if oct_fre > (2-octave_variability)*fre && oct_fre < (2+octave_variability)*fre
bleeck@4: 			rpeaks{i}.has_octave=oct_fre;
bleeck@4: 			rpeaks{i}.octave_peak_nr=j;
bleeck@4: 				
bleeck@4: 			% add the pss
bleeck@4: % 			rpeaks{j}.pitchstrength=rpeaks{i}.pitchstrength+rpeaks{j}.pitchstrength;
bleeck@4: 
bleeck@4: 			ps_real=rpeaks{j}.pitchstrength;
bleeck@4: 			ps_oct=rpeaks{i}.pitchstrength;
bleeck@4: 			hight_oct=rpeaks{j}.y;
bleeck@4: 			hight_real=rpeaks{i}.y;
bleeck@4: 			
bleeck@4: % 			if ps_oct>ps_real && hight_oct > hight_real
bleeck@4: % 				rpeaks{i}.pitchstrength=ps_real-1;	% artificially drop down
bleeck@4: % 			end
bleeck@4: 
bleeck@4: 			has_octave=1;
bleeck@4: 			break
bleeck@4: 		end
bleeck@4: 		if ~has_octave
bleeck@4: 			rpeaks{i}.has_octave=0;
bleeck@4: 			rpeaks{i}.octave_peak_nr=0;
bleeck@4: 		end
bleeck@4: 	end
bleeck@4: end
bleeck@4: 
bleeck@4: if plot_switch
bleeck@4: 	figure(2134)
bleeck@4: 	clf
bleeck@4: 	hold on
bleeck@4: 	plot(log_profile,'b')
bleeck@4: 	plot(smooth_sig,'g')
bleeck@4: 	for i=1:length(rpeaks)
bleeck@4: 		time=rpeaks{i}.t;
bleeck@4: 		x=rpeaks{i}.x;
bleeck@4: 		y=rpeaks{i}.y;
bleeck@4: 		plot(x,y,'Marker','o','Markerfacecolor','g','MarkeredgeColor','g','MarkerSize',2);
bleeck@4: 		time_left=rpeaks{i}.left.t;
bleeck@4: 		x_left=rpeaks{i}.left.x;
bleeck@4: 		y_left=rpeaks{i}.left.y;
bleeck@4: 		time_right=rpeaks{i}.right.t;
bleeck@4: 		x_right=rpeaks{i}.right.x;
bleeck@4: 		y_right=rpeaks{i}.right.y;
bleeck@4: 		plot(x_left,y_left,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',2);
bleeck@4: 		plot(x_right,y_right,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',2);
bleeck@4: 	end
bleeck@4: 	current_time=options.current_time*1000;
bleeck@4: 	y=get(gca,'ylim');
bleeck@4: 	y=y(2);
bleeck@4: 	text(700,y,sprintf('%3.0f ms',current_time),'verticalalignment','top');
bleeck@4: 	
bleeck@4: 	for i=1:length(rpeaks)
bleeck@4: 		pitchstrength=rpeaks{i}.pitchstrength;
bleeck@4: 		if i <10
bleeck@4: 			line([rpeaks{i}.x rpeaks{i}.x],[rpeaks{i}.y rpeaks{i}.peak_base_height_y],'Color','m');
bleeck@4: 			line([rpeaks{i}.where_widths(1) rpeaks{i}.where_widths(2)],[rpeaks{i}.peak_base_height_y rpeaks{i}.peak_base_height_y],'Color','m');
bleeck@4: 			x=rpeaks{i}.x;
bleeck@4: 			fre=rpeaks{i}.fre;
bleeck@4: 			y=rpeaks{i}.y;
bleeck@4: 			text(x,y*1.05,sprintf('%3.0fHz: %2.2f',fre,pitchstrength),'HorizontalAlignment','center');
bleeck@4: 		end
bleeck@4: 	end
bleeck@4: end
bleeck@4: if plot_switch==1
bleeck@4: 	for i=1:length(rpeaks)
bleeck@4: 		x=rpeaks{i}.x;
bleeck@4: 		y=rpeaks{i}.y;
bleeck@4: 		plot(x,y,'Marker','o','Markerfacecolor','g','MarkeredgeColor','g','MarkerSize',6);
bleeck@4: 		
bleeck@4: 		% plot the octaves as green lines
bleeck@4: 		octave=rpeaks{i}.has_octave;
bleeck@4: 		fre=rpeaks{i}.fre;
bleeck@4: 		if octave>0
bleeck@4: 			oct_nr=rpeaks{i}.octave_peak_nr;
bleeck@4: 			oct_fre=rpeaks{oct_nr}.fre;
bleeck@4: 			
bleeck@4: 			x=rpeaks{i}.x;
bleeck@4: 			oct_x=rpeaks{oct_nr}.x;
bleeck@4: 			y=rpeaks{i}.y;
bleeck@4: 			oct_y=rpeaks{oct_nr}.y;
bleeck@4: 			line([x oct_x],[y oct_y],'Color','g','LineStyle','--');
bleeck@4: 		end
bleeck@4: 	end
bleeck@4: end
bleeck@4: 
bleeck@4: % rpeaks=sortstruct(rpeaks,'pitchstrength');
bleeck@4: rpeaks=sortstruct(rpeaks,'y');
bleeck@4: 
bleeck@4: % final result for the full set
bleeck@4: result.final_pitchstrength=rpeaks{1}.pitchstrength;
bleeck@4: result.final_dominant_frequency=rpeaks{1}.fre;
bleeck@4: result.final_dominant_time=rpeaks{1}.t;
bleeck@4: result.smoothed_signal=smooth_sig;
bleeck@4: result.peaks=rpeaks;
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: function fre=x2fre(sig,x)
bleeck@4: t_log = bin2time(sig,x);
bleeck@4: t=f2f(t_log,0,0.035,0.001,0.035,'linlog');
bleeck@4: fre=1/t;