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