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

various bug fixes and changed copyright message

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