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

various bug fixes and changed copyright message

+% generating function for 'aim-mat'
+% 
+%   INPUT VALUES:
+%  
+%   RETURN VALUE:
+%
+% 
+% (c) 2003, University of Cambridge, Medical Research Council 
+% Christoph Lindner
+% (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 ps_result=genpitchstrength(sai,options)
+
+target_frequency=options.target_frequency;
+
+if target_frequency > 0
+	look_for_target_frequency=1;
+	minimum_allowed_frequency=target_frequency/options.allowed_frequency_deviation;
+	maximum_allowed_frequency=target_frequency*options.allowed_frequency_deviation;
+else
+	look_for_target_frequency=0;
+end
+
+if isfield(options,'resolved_harmonic_minimum')
+	resolved_harmonic_minimum=options.resolved_harmonic_minimum;
+end
+
+
+% find out about scaling:
+maxval=-inf;
+maxfreval=-inf;
+maxsumval=-inf;
+
+nr_frames=length(sai);
+for ii=1:nr_frames
+	maxval=max([maxval getmaximumvalue(sai{ii})]);
+	maxsumval=max([maxsumval getscalesumme(sai{ii})]);
+	maxfreval=max([maxfreval getscalefrequency(sai{ii})]);
+end
+
+
+for frame_number=1:nr_frames
+	current_frame = sai{frame_number};
+	% 	current_frame = setallmaxvalue(current_frame, maxval);
+	frame_start_time(frame_number)=getcurrentframestarttime(current_frame);
+	tip(frame_number)=gettimeintervalprofile(current_frame,options);
+end
+frame_duration=frame_start_time(2)-frame_start_time(1);
+
+% calculate the averaged tips
+% the first frames are not average, since not long enough time has
+% passed. After that, the relevant last frames are averaged
+tip_time=options.tip_average_time;
+% time for pitch to drop to 1/2 (6dB) --in seconds
+% this means so much per frame
+ps_decay_constant=1-log(2)/(tip_time/frame_duration);
+
+
+% % and again to sum them up
+% for frame_number=1:nr_frames
+% 	thisav=mute(tip(frame_number)); % take a copy of the original one
+% 	figure(1243234)
+% 	clf
+% 	hold on
+% 	count=0;
+% 	for recent_frames_nr=frame_number-nr_frames_to_average:frame_number
+% 		if recent_frames_nr > 0
+% 			thisav=thisav+tip(recent_frames_nr);
+% 			count=count+1;
+% 			plot(tip(recent_frames_nr),'g');
+% 		end
+% 	end
+% 	thisav=thisav/count;	% the average
+% 	plot(thisav,'r');
+% 	average_tip(frame_number)=thisav;
+% end
+
+% is it necessary to calculate all, or can we start right at the end?
+% if options.return_only_last==1
+% 	start_frame=length(sai);
+% else	% no we 
+start_frame=1;
+% end
+
+emptysig=mute(tip(1));
+% mache ein Histogramm für jede Frequenz mit allem Maxima:
+nr_points=getnrpoints(sai{1});
+histo=zeros(1,nr_points);
+
+waithand = waitbar(0,'Generating pitch strength');
+for frame_number=start_frame:nr_frames
+% for frame_number=1:10
+	
+	waitbar(frame_number/nr_frames, waithand);     
+	
+	current_frame = sai{frame_number};
+	current_frame = setallmaxvalue(current_frame, maxval);
+	current_frame = setscalesumme(current_frame, maxsumval);
+	current_frame = setscalefrequency(current_frame, maxfreval);
+	
+	interval_sig =  tip(frame_number);
+	% 	interval_sig =  average_tip(frame_number);      
+	ps_result{frame_number}.frequency_sum = getfrequencysum(current_frame);      
+	% Normalisation
+	interval_sig =  interval_sig/getnrpoints( ps_result{frame_number}.frequency_sum);
+	ps_result{frame_number}.frequency_sum = ps_result{frame_number}.frequency_sum/getnrpoints(interval_sig)*options.scalefactor;
+	
+	opts.ps_options=options;
+	% do all relevant pitch calculations for this frame
+	result=find_pitches(interval_sig,opts);
+	
+	ps_result{frame_number}.interval_sum =emptysig;
+	ps_result{frame_number}.pitchstrength=0;
+
+
+	% sum all pitches for a histogram
+	peaks=result.peaks;
+	nr_peaks=length(peaks);
+	if nr_peaks>0
+		for current_peak=1:nr_peaks
+			peak_time=peaks{current_peak}.t;
+			peak_x=peaks{current_peak}.x;
+			peak_height=peaks{current_peak}.y;
+			int_x=round(peak_x);
+			histo(int_x)=histo(int_x)+peak_height;
+		end
+		
+		shist=signal(histo,getsr(interval_sig));
+		shistl=lowpass(shist,400);
+		peaks=IPeakPicker(shistl);
+
+		% ignore the first peak on the left:
+% 		speaks=sortstruct(peaks,'x');
+% 		peaks{1}.y=0;
+		
+		speaks=sortstruct(peaks,'y');
+		
+		pitchstrength=speaks{1}.y;
+		final_dominant_frequency=1/speaks{1}.t;
+% 		
+% 		if mod(frame_number,10)==0
+% 			figure(43242)
+% 			if frame_number<30
+% 				clf;
+% 				hold on
+% 			end
+% 			plot(shist/50);
+% % 			ax=axis;
+% 			plot(shistl,'r'); hold on
+% % 			axis(ax);
+% 			fre=1/speaks{1}.t;
+% 			plot(speaks{1}.x,speaks{1}.y,'o','MarkerFaceColor','g','MarkerEdgeColor','w','MarkerSize',10)
+% 			text(speaks{1}.x,speaks{1}.y,sprintf('%2.0fHz %2.2f',final_dominant_frequency,pitchstrength))
+% 			for i=2:length(speaks)/2
+% 				fre=1/speaks{i}.t;
+% 				plot(speaks{i}.x,speaks{i}.y,'o','MarkerFaceColor','m','MarkerEdgeColor','w','MarkerSize',5)
+% 				text(speaks{i}.x,speaks{i}.y,sprintf('%2.0fHz %2.2f',fre,speaks{i}.y))
+% 			end
+% 			time=frame_number*5;
+% 			text(speaks{1}.x+400,speaks{1}.y,sprintf('Time: %3.0fms',time));
+% 			o=0;
+% 		end
+		
+		ps_result{frame_number}.peaks = speaks;
+ 		ps_result{frame_number}.interval_sum =shistl;
+		
+		if look_for_target_frequency
+			count=0;
+			for ii=1:length(speaks) % construct all maxima
+				fre=1/speaks{ii}.t;
+				dist(ii)=abs(target_frequency-fre);
+				if fre> minimum_allowed_frequency && fre < maximum_allowed_frequency
+					count=count+1;
+					allowednumber(count)=ii;
+					allowedheights(count)=speaks{ii}.y;
+				end
+			end
+% 			[minfredif,nearest_peak_number]=min(dist);
+% 			fre=1/speaks{nearest_peak_number}.t;
+			% is the frequency allowed?
+			if count>0
+% 			if fre< minimum_allowed_frequency || fre > maximum_allowed_frequency
+				[highest_ps,best_number]=max(allowedheights);
+				best_peak_number=allowednumber(best_number);
+				
+				% gib das Verhältnis der peakstrength zu allen peaks
+				% zurück (ohne ihm selbst)
+				do_only_relationship=1;
+				if do_only_relationship
+					count=1;
+					for ii=1:length(speaks)
+						if ii~=best_peak_number
+							allps(count)=speaks{ii}.y;
+							count=count+1;
+						end
+					end
+					ps=speaks{best_peak_number}.y/mean(allps);
+					ps_result{frame_number}.pitchstrength=ps;
+					ps_result{frame_number}.dominant_frequency=1/speaks{best_peak_number}.t;
+				else
+ 					ps_result{frame_number}.pitchstrength=speaks{best_peak_number}.y;
+					ps_result{frame_number}.dominant_frequency=1/speaks{best_peak_number}.t;
+				end
+			else % no, nothing there
+% 				ps_result{frame_number}.interval_sum = result.smoothed_signal;
+				ps_result{frame_number}.peaks = [];
+				ps_result{frame_number}.dominant_frequency=0;
+				ps_result{frame_number}.pitchstrength=0;
+			end % frequency allowance
+		else % free search was already done!
+			ps_result{frame_number}.dominant_frequency=final_dominant_frequency;
+			ps_result{frame_number}.pitchstrength=pitchstrength;
+		end % look for target frequency
+	else % no peaks found
+		ps_result{frame_number}.peaks = [];
+		ps_result{frame_number}.dominant_frequency=0;
+		ps_result{frame_number}.pitchstrength=0;
+	end
+	
+	% dynamic decrease of all pitch strengthes
+	histo=histo.*ps_decay_constant; 
+	
+	
+	% Peak Picker for FrequencyProfile
+	p.dyn_thresh = options.dynamic_threshold_FP;
+	p.smooth_sigma = options.smooth_sigma_FP;
+	ps_result{frame_number}.peaks_frequency_sum = FPeakPicker(ps_result{frame_number}.frequency_sum, p);
+	% 	ps_result{frame_number}.channel_center_fq = getcf(sai{frame_number});
+	ps_result{frame_number}.channel_center_fq = getcf(sai{frame_number});
+	
+end
+close(waithand);
+
+