comparison aim-mat/modules/usermodule/pitchstrength/find_pitches2.m @ 4:537f939baef0 tip

various bug fixes and changed copyright message

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