view aim-mat/modules/usermodule/pitchstrength/find_pitches.asv @ 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) 2003, University of Cambridge, Medical Research Council 
% Stefan Bleeck ([e-mail])
% http://www.mrc-cbu.cam.ac.uk/cnbh/aimmanual
% $Date: 2003/07/17 10:56:16 $
% $Revision: 1.5 $
function result = find_pitches(profile,options)
% different ways to define the pitch strength:
% 1: the absolute height of the highest peak
% 2: ratio between peak hight and width devided at base
% % % % 2: the ration of the highest peak divided by the width at a certain point (given
% % % %     by the parameter height_to_width_ratio
% 3: the height of the highest peak divided by the hight of the peak just one to 
%     the right. This measurement is very successfull in the ramped/damped stimuli
% 
% 4: Further we want to know all these values at a fixed point called target_frequency
%     and in the range given by allowed_frequency_deviation in %
%     
% 5: We are also interested in the frequency value of the highest peak, because
%     we hope, that this is finally the pitch.
%     
% 6: for the dual profile model, we also give back two more values concerning similar
%     properties for the spectral profile. Here, the pitch strenght is defined as the 
%     height of the highest peak divided by the range that is given in two parameters 
%     (usually 20% to 80% of the maximum)


plot_switch =1;

if nargin < 2
    options=[];
end
% 

% peaks must be higher then this of the maximum to be recognised
ps_threshold=0.1;   
height_width_ratio=options.ps_options.height_width_ratio;	% height of peak, where the width is measured


% preliminary return values.
% height of the highest peak
% result.free.highest_peak_hight=0;
% result.free.highest_peak_frequency=0;
% hight to width ratio
% result.free.height_width_ratio=0;
% highest peak divided by next highest
% result.free.neigbouring_ratio=0;

% now all these at a fixed frequency:
% % result.fixed.highest_peak_hight=0;
% result.fixed.highest_peak_frequency=0;
% result.fixed.height_width_ratio=0;
% result.fixed.neigbouring_ratio=0;
% 
% % other things useful for plotting
% result.smoothed_signal=[];
% result.peaks = [];


% now start the show

% % change the scaling to logarithm, before doing anything else:
% log_profile=logsigx(profile,0.001,0.035);
% result.smoothed_signal=log_profile;

% don't do this change
log_profile=profile;

current_lowpass_frequency=options.ps_options.low_pass_frequency;
smooth_sig=lowpass(log_profile,current_lowpass_frequency);
envpeaks = PeakPicker(smooth_sig);
result.smoothed_signal=smooth_sig;

if isempty(envpeaks) % only, when there is no signal
    return
end

% reject impossible peaks
ep=envpeaks;ep2=[];c=1;
for i=1:length(envpeaks);
    if envpeaks{i}.t>0.002 && envpeaks{i}.t<0.03 && envpeaks{i}.y>0.01
        ep2{c}=ep{i};
        c=c+1;
    end
end
envpeaks=ep2; % replace
if isempty(envpeaks) % only, when there is no signal
    return
end
%
% % translate times back to times:
% for i=1:length(envpeaks) % construct all maxima
%     envpeaks{i}.t=1/x2fre(smooth_sig,envpeaks{i}.x);
% end


% nr1: highest peak value
% find the highest peak and its frequency
% sort all for the highest first!
% envpeaks=sortstruct(envpeaks,'y');
% result.free.highest_peak_frequency=1/envpeaks{1}.t;
% result.free.highest_peak_hight=envpeaks{1}.y;



% SB 08.2012: adjusted the method to make it simpler for Daniels signals of
% IRN
% nr 2: height to width of each peak. Height=height of peak. Width = width
% of the two adjacent minima
for i=1:length(envpeaks) % construct all maxima
    where=envpeaks{i}.t;
    hp=envpeaks{i}.y; % height peak
    hb=(envpeaks{i}.left.y+envpeaks{i}.right.y)/2;% base peak
    diffheight=hp-hb;
    w=log(envpeaks{i}.right.t)-log(envpeaks{i}.left.t); % width at base
    if w>0 && diffheight>0.
        envpeaks{i}.v2012_height_base_width_ratio=diffheight/w;
    else
        envpeaks{i}.v2012_height_base_width_ratio=0;
    end
    envpeaks{i}.v2012_base_width=w;
    envpeaks{i}.v2012_base_where_widths=hb; %base height
end

envpeaks=sortstruct(envpeaks,'v2012_height_base_width_ratio');

if plot_switch
    figure(2134)
    clf
    hold on
    plot(log_profile,'r');
    plot(smooth_sig,'b')
    for i=1:min(5,length(envpeaks))
%         time=envpeaks{i}.t;
%         x=envpeaks{i}.x;
%         y=envpeaks{i}.y;
%         plot(time,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(time_left,y_left,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',2);
%         plot(time_right,y_right,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',2);
        
		t=envpeaks{i}.t;
		ypeak=envpeaks{i}.y;
% 		ps=peaks{ii}.pitchstrength;
		ps=envpeaks{i}.v2012_height_base_width_ratio;
        
		if i==1
			plot(t,ypeak,'Marker','o','Markerfacecolor','b','MarkeredgeColor','b','MarkerSize',10);
 			text(t,ypeak*1.05,sprintf('%3.0f Hz: %4.2f ',1/t,ps),'VerticalAlignment','bottom','HorizontalAlignment','center','color','b','Fontsize',12);
		else
			plot(t,ypeak,'Marker','o','Markerfacecolor','g','MarkeredgeColor','w','MarkerSize',5);
			text(t,ypeak*1.05,sprintf('%3.0f Hz: %4.2f ',1/t,ps),'VerticalAlignment','bottom','HorizontalAlignment','center','color','g','Fontsize',12);
        end
        plot(envpeaks{i}.left.t,envpeaks{i}.left.y,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',5);
        plot(envpeaks{i}.right.t,envpeaks{i}.right.y,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',5);
        
        
        ybase=envpeaks{i}.v2012_base_where_widths;
		line([t t],[ybase ypeak],'color','m');
		line([envpeaks{i}.left.t envpeaks{i}.right.t],[ybase ybase],'color','m');
        
    end
    
%          set(gca,'xscale','log')
    set(gca,'xlim',[0.001 0.03])
    
    fres=[500 300 200 150 100 70 50 20];
    set(gca,'xtick',1./fres);
    set(gca,'xticklabel',fres);
    xlabel('Frequency (Hz)')
    ylabel('arbitrary normalized units')
    
%     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)
%         height_width_ratio=envpeaks{i}.height_width_ratio;
%         if i <4
%             line([envpeaks{i}.t envpeaks{i}.x],[envpeaks{i}.y envpeaks{i}.y*(1-height_width_ratio)],'Color','r');
%             line([envpeaks{i}.where_widths(1) envpeaks{i}.where_widths(2)],[envpeaks{i}.y*(1-height_width_ratio) envpeaks{i}.y*(1-height_width_ratio)],'Color','r');
%             x=envpeaks{i}.t;
%             fre=1/envpeaks{i}.t;
%             y=envpeaks{i}.y;
%             text(x,y*1.05,sprintf('%3.0fHz: %2.2f',fre,height_width_ratio),'HorizontalAlignment','center');
%         end
%     end
end


% and return the final result, only highest peak
peaks2=sortstruct(envpeaks,'v2012_height_base_width_ratio');
result.free.ps=peaks2{1}.v2012_height_base_width_ratio;
result.free.fre=1/peaks2{1}.t;



% 
% 
% % nr 2: height to width of highest peak
% for i=1:length(envpeaks) % construct all maxima
% %     where=bin2time(smooth_sig,envpeaks{i}.x);
%     where=envpeaks{i}.t;
%     [~,height,breit,where_widths]=qvalue(smooth_sig,where,height_width_ratio);
%     width=time2bin(smooth_sig,breit);
%     if width>0
%         envpeaks{i}.height_width_ratio=height/width;
%     else
%         envpeaks{i}.height_width_ratio=0;
%     end
%     envpeaks{i}.width=width;
%     envpeaks{i}.where_widths=where_widths;
%     envpeaks{i}.peak_base_height_y=height*(1-height_width_ratio);
% end
% % and return the results
% % result.free.height_width_ratio=envpeaks{1}.height_width_ratio;

% 
% % nr 3: height of highest / right neigbour (right when time is towards
% % left, sorry)
% for i=1:length(envpeaks) % construct all maxima
%     left=envpeaks{i}.left;
%     if isfield(left,'y') && left.y>0
%         envpeaks{i}.neigbouring_ratio=result.free.highest_peak_hight/left.y;
%     else
%         envpeaks{i}.neigbouring_ratio=0;
%     end
% end
% result.free.neigbouring_ratio=envpeaks{1}.neigbouring_ratio;


target_frequency=options.ps_options.target_frequency;
% now find all values for the fixed pitch strengh in target_frequency
if target_frequency>0 % only, when wanted
    min_fre=target_frequency/options.ps_options.allowed_frequency_deviation;
    max_fre=target_frequency*options.ps_options.allowed_frequency_deviation;
    
    for i=1:length(envpeaks) % look through all peaks, which one we need
        fre_peak=1/envpeaks{i}.t;
        if fre_peak > min_fre && fre_peak < max_fre
            % we assume for the moment, that we only have one allowed here
            % nr 1: height
            result.fixed.highest_peak_frequency=fre_peak;
            result.fixed.highest_peak_hight=envpeaks{i}.y;
            result.fixed.height_width_ratio=envpeaks{i}.height_width_ratio;
            result.fixed.neigbouring_ratio=envpeaks{i}.neigbouring_ratio;
        end
    end
end

result.peaks =envpeaks;



return





% 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;