changeset 4:537f939baef0 tip

various bug fixes and changed copyright message
author Stefan Bleeck <bleeck@gmail.com>
date Tue, 16 Aug 2011 14:37:17 +0100
parents 20ada0af3d7d
children
files aim-mat/gui/length_questionbox.fig aim-mat/gui/sr_questionbox.fig aim-mat/gui/sr_questionbox.m aim-mat/modules/sai/none/gennosai.m aim-mat/modules/sai/none/gennostrobes.asv aim-mat/modules/sai/none/gennostrobes.m aim-mat/modules/sai/none/parameters.m aim-mat/modules/usermodule/pitchstrength/FPeakPicker.m aim-mat/modules/usermodule/pitchstrength/IPeakPicker.m aim-mat/modules/usermodule/pitchstrength/PeakPicker.m aim-mat/modules/usermodule/pitchstrength/analyse_aim_profiles.m aim-mat/modules/usermodule/pitchstrength/analyse_frequency_profile.m aim-mat/modules/usermodule/pitchstrength/analyse_timeinterval_profile.m aim-mat/modules/usermodule/pitchstrength/displaypitchstrength.asv aim-mat/modules/usermodule/pitchstrength/displaypitchstrength.m aim-mat/modules/usermodule/pitchstrength/find_pitches.asv aim-mat/modules/usermodule/pitchstrength/find_pitches.m aim-mat/modules/usermodule/pitchstrength/find_pitches2.m aim-mat/modules/usermodule/pitchstrength/genpitchstrength.m aim-mat/modules/usermodule/pitchstrength/genpitchstrength2.m aim-mat/modules/usermodule/pitchstrength/genpitchstrength3.m aim-mat/modules/usermodule/pitchstrength/parameters.m aim-mat/modules/usermodule/pitchstrength/reduce_tilt.m aim-mat/tools/xticklabel_rotate.m
diffstat 24 files changed, 3721 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
Binary file aim-mat/gui/length_questionbox.fig has changed
Binary file aim-mat/gui/sr_questionbox.fig has changed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/gui/sr_questionbox.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,145 @@
+% 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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 varargout = sr_questionbox(varargin)
+% SR_QUESTIONBOX M-file for sr_questionbox.fig
+%      SR_QUESTIONBOX, by itself, creates a new SR_QUESTIONBOX or raises the existing
+%      singleton*.
+%
+%      H = SR_QUESTIONBOX returns the handle to a new SR_QUESTIONBOX or the handle to
+%      the existing singleton*.
+%
+%      SR_QUESTIONBOX('CALLBACK',hObject,eventData,handles,...) calls the local
+%      function named CALLBACK in SR_QUESTIONBOX.M with the given input arguments.
+%
+%      SR_QUESTIONBOX('Property','Value',...) creates a new SR_QUESTIONBOX or raises the
+%      existing singleton*.  Starting from the left, property value pairs are
+%      applied to the GUI before sr_questionbox_OpeningFunction gets called.  An
+%      unrecognized property name or invalid value makes property application
+%      stop.  All inputs are passed to sr_questionbox_OpeningFcn via varargin.
+%
+%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
+%      instance to run (singleton)".
+%
+% See also: GUIDE, GUIDATA, GUIHANDLES
+
+% Edit the above text to modify the response to help sr_questionbox
+
+% Last Modified by GUIDE v2.5 14-Mar-2003 17:30:37
+
+% Begin initialization code - DO NOT EDIT
+gui_Singleton = 1;
+gui_State = struct('gui_Name',       mfilename, ...
+                   'gui_Singleton',  gui_Singleton, ...
+                   'gui_OpeningFcn', @sr_questionbox_OpeningFcn, ...
+                   'gui_OutputFcn',  @sr_questionbox_OutputFcn, ...
+                   'gui_LayoutFcn',  [] , ...
+                   'gui_Callback',   []);
+if nargin & isstr(varargin{1})
+    gui_State.gui_Callback = str2func(varargin{1});
+end
+
+if nargout
+    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
+else
+    gui_mainfcn(gui_State, varargin{:});
+end
+% End initialization code - DO NOT EDIT
+
+
+% --- Executes just before sr_questionbox is made visible.
+function sr_questionbox_OpeningFcn(hObject, eventdata, handles, varargin)
+% This function has no output args, see OutputFcn.
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+% varargin   command line arguments to sr_questionbox (see VARARGIN)
+
+% Choose default command line output for sr_questionbox
+handles.output = hObject;
+
+old_sr=varargin{1};
+handles.old_sr=old_sr;
+handles.output=old_sr;
+
+set(handles.edit1,'String',old_sr);
+
+
+% Update handles structure
+guidata(hObject, handles);
+
+% do modal
+% UIWAIT makes sr_questionbox wait for user response (see UIRESUME)
+uiwait(handles.figure1);
+
+
+% --- Outputs from this function are returned to the command line.
+function varargout = sr_questionbox_OutputFcn(hObject, eventdata, handles)
+% varargout  cell array for returning output args (see VARARGOUT);
+% hObject    handle to figure
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Get default command line output from handles structure
+varargout{1} = handles.output;
+close(handles.figure1);
+
+
+% --- Executes during object creation, after setting all properties.
+function edit1_CreateFcn(hObject, eventdata, handles)
+% hObject    handle to edit1 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    empty - handles not created until after all CreateFcns called
+
+% Hint: edit controls usually have a white background on Windows.
+%       See ISPC and COMPUTER.
+if ispc
+    set(hObject,'BackgroundColor','white');
+else
+    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
+end
+
+
+
+function edit1_Callback(hObject, eventdata, handles)
+% hObject    handle to edit1 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+
+% Hints: get(hObject,'String') returns contents of edit1 as text
+%        str2double(get(hObject,'String')) returns contents of edit1 as a double
+
+num=str2num(get(handles.edit1,'String'));
+if isempty(num)
+	num=handles.old_sr;
+end
+if num<1
+	num=handles.old_sr;
+end
+set(handles.edit1,'String',num2str(num));
+
+handles.output=num;
+guidata(hObject, handles);
+
+% --- Executes on button press in pushbutton1.
+function pushbutton1_Callback(hObject, eventdata, handles)
+% hObject    handle to pushbutton1 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+handles.output=str2num(get(handles.edit1,'String'));
+guidata(hObject, handles);
+uiresume;
+
+% --- Executes on button press in pushbutton2.
+function pushbutton2_Callback(hObject, eventdata, handles)
+% hObject    handle to pushbutton2 (see GCBO)
+% eventdata  reserved - to be defined in a future version of MATLAB
+% handles    structure with handles and user data (see GUIDATA)
+handles.output=handles.old_sr;
+guidata(hObject, handles);
+uiresume;
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/sai/none/gennosai.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,16 @@
+% generating function for 'aim-mat'
+% 
+%   INPUT VALUES:
+%  
+%   RETURN VALUE:
+%
+% 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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 [returnframes]=gensai(nap,strobes,options)
+returnframes=[];
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/sai/none/gennostrobes.asv	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,16 @@
+% generating function for 'aim-mat'
+% 
+%   INPUT VALUES:
+%  
+%   RETURN VALUE:
+%
+% 
+% (c) 2003-2008, University of Cambridge, Medical Research Council 
+% Maintained by Tom Walters (tcw24@cam.ac.uk), written by Stefan Bleeck (stefan@bleeck.de)
+% http://www.pdn.cam.ac.uk/cnbh/aim2006
+% $Date: 2008-06-10 18:00:16 +0100 (Tue, 10 Jun 2008) $
+% $Revision: 585 $
+
+function [strobes,thres]=gennosai(nap,options)
+strobes=[];
+thres=[];
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/sai/none/gennostrobes.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,17 @@
+% generating function for 'aim-mat'
+% 
+%   INPUT VALUES:
+%  
+%   RETURN VALUE:
+%
+% 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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 [strobes,thres]=gennosai(nap,options)
+strobes=[];
+thres=[];
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/sai/none/parameters.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,17 @@
+% parameter file for 'aim-mat'
+% 
+% 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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
+
+%%%%%%%%%%%%%
+% no sai
+% hidden parameters
+none.generatingfunction='gennosai';
+none.displayname='no sai';
+none.revision='$Revision: 585 $';
+
+% parameters relevant for the calculation of this module
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/FPeakPicker.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,143 @@
+%  
+% function out = PeakPicker(sig_in, params)
+%
+%   Find the peaks of a signal
+%
+%   INPUT VALUES:
+%       sig_in                  Input signal
+%       params.dyn_thresh       dynamic threshold. Off if not used
+%       params.smooth_sigma     sigma for smoothing
+%                               
+%       
+% 		
+%
+%   RETURN VALUE:
+%       out is an array of a struct
+%       out.x                    x position of the Peak
+%       out.t                    according time value
+%		out.y                    y value of the peak
+%       out.left.[x,t,y]         left Minumum
+%       out.right.[x,t,y]        right Minumum
+% 
+% (c) 2003, University of Cambridge, Medical Research Council 
+% Christoph Lindner 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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 out = PeakPicker(sig_in, params)
+
+
+% ----- Other Parameters -----
+% Lowpass param for the higpassfilter
+LP_sigma_for_HP_filter = getnrpoints(sig_in)/7;
+LP_sigma_for_smooth = 3;
+% min width of a peak: upper_threshold+lower_thresh
+upper_thresh =  0.02*getnrpoints(sig_in);
+lower_thresh =  0.03*getnrpoints(sig_in);
+
+% x is index or position in vector
+% y is value of the
+% t is the time domain wich is assigned to the x dimension
+
+% the original values befor filtering
+orig_values = getdata(sig_in)';
+
+if isfield(params,'smooth_sigma') 
+    if (params.smooth_sigma~=0)
+        % smooth the curve to kill small side peaks 
+        sig_in = smooth(sig_in, params.smooth_sigma);
+    end
+end
+
+values=getdata(sig_in)';
+
+% ------------------- Find the local maxima ------------------------------
+% find x positions of ALL local maxima, incl. zero!!
+max_x = find((values >= [0 values(1:end-1)]) & (values > [values(2:end) 0]));
+if isfield(params,'smooth_sigma') 
+    if (params.smooth_sigma~=0)
+        % smoothing might have shifted the positions of maxima.
+        % Therefore the maximum is the highest of the neighbours in
+        % distance +- smooth_sigma
+        for i=1:length(max_x)
+            start = max_x(i)-params.smooth_sigma;
+            if start<1
+                start=1;
+            end
+            stop = max_x(i)+params.smooth_sigma;
+            if stop>length(orig_values)
+                stop=length(orig_values);
+            end
+            m = find(orig_values(start:stop) == max(orig_values(start:stop)));
+            m=m(1);
+            max_x(i)=start-1+m;
+        end
+    end
+end
+max_y = orig_values(max_x);
+orig_max_y = orig_values(max_x);
+
+% ------------------- Find the local minima -----------------------------
+min_x = find((values < [inf values(1:end-1)]) & (values <= [values(2:end) inf]));
+min_y = values(min_x);
+
+
+% peakpos_x=zeros(1,length(max_x));
+peakpos_x=[];
+for i=1:length(max_x),
+    % only take the highest peak 
+    my = [max_y==max(max_y)];  % find the highest peak
+%     peakpos_x(i) = max_x(my);   % x pos of highest peak
+    peakpos_x = [peakpos_x max_x(my)]; 
+    max_y = max_y([max_y<max(max_y)]);   % del max value in y domain
+    max_x = max_x([max_x ~= peakpos_x(end)]); % and in x domain
+end
+peakpos_y = orig_values(peakpos_x);  % extract the y vector 
+
+% --------------- Dynamic Threshold ---------------------
+% works relativ to the mean
+if isfield(params,'dyn_thresh') 
+    if (params.dyn_thresh~=0)
+        % dynamic thresholding 
+        m = mean(orig_values);  
+        dthr = params.dyn_thresh.*m;
+        peakpos_x = peakpos_x([peakpos_y>=dthr]);
+        peakpos_y = peakpos_y([peakpos_y>=dthr]);
+    end
+end
+
+maxima = cell(1, length(peakpos_x));
+% find the left end right minima that belong to a maximum
+for i=1:length(peakpos_x)
+    maxima{i}.x = peakpos_x(i);
+    maxima{i}.t = bin2time(sig_in, maxima{i}.x);
+    maxima{i}.y = orig_values(peakpos_x(i));
+    
+    % find left and right minimum for this maximum
+    maxima{i}.left.x = max(min_x([min_x < maxima{i}.x]));
+    if isempty(maxima{i}.left.x)
+        maxima{i}.left.x = 1;
+        maxima{i}.left.t = 0;
+        maxima{i}.left.y = orig_values(maxima{i}.left.x);
+    else
+        maxima{i}.left.y = orig_values(maxima{i}.left.x);
+        maxima{i}.left.t = bin2time(sig_in, maxima{i}.left.x);
+    end
+    maxima{i}.right.x = min(min_x([min_x > maxima{i}.x]));
+    if isempty(maxima{i}.right.x)
+        maxima{i}.right.x = length(orig_values);
+        maxima{i}.right.t = 0;
+        maxima{i}.right.y = orig_values(maxima{i}.right.x);
+    else
+        maxima{i}.right.y = orig_values(maxima{i}.right.x);
+        maxima{i}.right.t = bin2time(sig_in, maxima{i}.right.x);
+    end
+end
+
+
+out = maxima;
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/IPeakPicker.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,145 @@
+% function out = PeakPicker(sig_in, params)
+%
+%   Find the peaks of a signal and their neighbours
+%
+%   INPUT VALUES:
+%       sig_in                  Input signal
+%         threshold     dynamic threshold. Off if not used
+%                               
+%
+%   RETURN VALUE:
+%       out is an array of a struct
+%       out.x                    x position of the Peak
+%       out.t                    according time value
+%		out.y                    y value of the peak
+%       out.left.[x,t,y]         left Minumum
+%       out.right.[x,t,y]        right Minumum
+% 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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 out = PeakPicker(sig_in,threshold)
+% threshold is the percentage of the maximum value of the signal,
+% under which a peak is not counted
+
+if nargin <2
+	threshold=0;
+end
+
+plot_switch = 0;
+
+% the original values befor filtering
+orig_values = getdata(sig_in)';
+values=getdata(sig_in)';
+
+% ------------------- Find the local maxima ------------------------------
+% find x positions of ALL local maxima, incl. zero!!
+max_x = find((values >= [0 values(1:end-1)]) & (values > [values(2:end) 0]));
+max_y = orig_values(max_x);
+orig_max_y = orig_values(max_x);
+
+% ------------------- Find the local minima -----------------------------
+min_x = find((values < [inf values(1:end-1)]) & (values <= [values(2:end) inf]));
+min_y = values(min_x);
+
+
+peakpos_x=[];
+for i=1:length(max_x),
+    % only take the highest peak 
+    my = [max_y==max(max_y)];  % find the highest peak
+%     peakpos_x(i) = max_x(my);   % x pos of highest peak
+    peakpos_x = [peakpos_x max_x(my)]; 
+    max_y = max_y([max_y<max(max_y)]);   % del max value in y domain
+    max_x = max_x([max_x ~= peakpos_x(end)]); % and in x domain
+end
+peakpos_y = orig_values(peakpos_x);  % extract the y vector 
+
+% maxima = cell(1, length(peakpos_x));
+maxima = [];
+
+if plot_switch
+	figure(123);
+	clf
+	plot(sig_in,'k');
+	hold on;
+end
+
+threshold_val=threshold*max(sig_in);
+counter=1;
+
+% find the left end right minima that belong to a maximum
+for i=1:length(peakpos_x)
+	y_val= orig_values(peakpos_x(i));
+	
+	if y_val< threshold_val
+		continue
+	end
+
+	maxima{counter}.y =y_val;
+	maxima{counter}.x = peakpos_x(i);
+	maxima{counter}.t = bin2time(sig_in, maxima{i}.x);
+	maxima{counter}.fre = 1/maxima{counter}.t;
+	
+	if plot_switch
+		plot(maxima{counter}.x,maxima{counter}.y,'go');
+	end    
+	% find left and right minimum for this maximum
+	maxima{counter}.left.x = max(min_x([min_x < maxima{counter}.x]));
+	if isempty(maxima{counter}.left.x)
+		maxima{counter}.left.x = 1;
+		maxima{counter}.left.t = 0;
+		maxima{counter}.left.y = orig_values(maxima{counter}.left.x);
+	else
+		maxima{counter}.left.y = orig_values(maxima{counter}.left.x);
+		maxima{counter}.left.t = bin2time(sig_in, maxima{counter}.left.x);
+	end
+	maxima{counter}.right.x = min(min_x([min_x > maxima{counter}.x]));
+	if isempty(maxima{counter}.right.x)
+		maxima{counter}.right.x = length(orig_values);
+		maxima{counter}.right.t = 0;
+		maxima{counter}.right.y = orig_values(maxima{counter}.right.x);
+	else
+		maxima{counter}.right.y = orig_values(maxima{counter}.right.x);
+		maxima{counter}.right.t = bin2time(sig_in, maxima{counter}.right.x);
+	end
+	
+	if plot_switch
+		plot(maxima{counter}.right.x,maxima{counter}.right.y,'ro');
+		plot(maxima{counter}.left.x,maxima{counter}.left.y,'ro');
+	end
+	counter=counter+1;
+end
+
+
+% umrechnung in die Darstellung, die wir brauchen:
+for i=1:counter-1
+	logtime=maxima{i}.t;
+	time=logtime2time(logtime);
+	maxima{i}.t=time;
+	maxima{i}.fre=1/time;
+	maxima{i}.y=gettimevalue(sig_in,logtime);
+
+	left=maxima{i}.left;
+	lefttime=logtime2time(left.t);
+	maxima{i}.left.t=lefttime;
+	maxima{i}.left.fre=1/lefttime;
+	maxima{i}.left.y=gettimevalue(sig_in,left.t);
+
+	right=maxima{i}.right;
+	righttime=logtime2time(right.t);
+	maxima{i}.right.t=righttime;
+	maxima{i}.right.fre=1/righttime;
+	maxima{i}.right.y=gettimevalue(sig_in,right.t);
+end
+
+
+out = maxima;
+
+
+
+function time=logtime2time(logtime)
+	time=f2f(logtime,0,0.035,0.001,0.035,'linlog');
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/PeakPicker.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,146 @@
+%  
+% function out = PeakPicker(sig_in, params)
+%
+%   Find the peaks of a signal
+%
+%   INPUT VALUES:
+%       sig_in                  Input signal
+%       params.dyn_thresh       dynamic threshold. Off if not used
+%       params.smooth_sigma     sigma for smoothing
+%                               
+%       
+% 		
+%
+%   RETURN VALUE:
+%       out is an array of a struct
+%       out.x                    x position of the Peak
+%       out.t                    according time value
+%		out.y                    y value of the peak
+%       out.left.[x,t,y]         left Minumum
+%       out.right.[x,t,y]        right Minumum
+% 
+% (c) 2003, University of Cambridge, Medical Research Council 
+% Christoph Lindner 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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 out = PeakPicker(sig_in, params)
+
+if nargin<2
+    params=[];
+end
+
+% % % % ----- Other Parameters -----
+% % % % Lowpass param for the higpassfilter
+% % % LP_sigma_for_HP_filter = getnrpoints(sig_in)/7;
+% % % LP_sigma_for_smooth = 3;
+% % % % min width of a peak: upper_threshold+lower_thresh
+% % % upper_thresh =  0.02*getnrpoints(sig_in);
+% % % lower_thresh =  0.03*getnrpoints(sig_in);
+
+% x is index or position in vector
+% y is value of the
+% t is the time domain wich is assigned to the x dimension
+
+% the original values befor filtering
+orig_values = getdata(sig_in)';
+
+if isfield(params,'smooth_sigma') 
+    if (params.smooth_sigma~=0)
+        % smooth the curve to kill small side peaks 
+        sig_in = smooth(sig_in, params.smooth_sigma);
+    end
+end
+
+values=getdata(sig_in)';
+
+% ------------------- Find the local maxima ------------------------------
+% find x positions of ALL local maxima, incl. zero!!
+max_x = find((values >= [0 values(1:end-1)]) & (values > [values(2:end) 0]));
+if isfield(params,'smooth_sigma') 
+    if (params.smooth_sigma~=0)
+        % smoothing might have shifted the positions of maxima.
+        % Therefore the maximum is the highest of the neighbours in
+        % distance +- smooth_sigma
+        for i=1:length(max_x)
+            start = max_x(i)-params.smooth_sigma;
+            if start<1
+                start=1;
+            end
+            stop = max_x(i)+params.smooth_sigma;
+            if stop>length(orig_values)
+                stop=length(orig_values);
+            end
+            m = find(orig_values(start:stop) == max(orig_values(start:stop)));
+            m=m(1);
+            max_x(i)=start-1+m;
+        end
+    end
+end
+max_y = orig_values(max_x);
+orig_max_y = orig_values(max_x);
+
+% ------------------- Find the local minima -----------------------------
+min_x = find((values < [inf values(1:end-1)]) & (values <= [values(2:end) inf]));
+min_y = values(min_x);
+
+
+% peakpos_x=zeros(1,length(max_x));
+peakpos_x=[];
+for i=1:length(max_x),
+    % only take the highest peak 
+    my = [max_y==max(max_y)];  % find the highest peak
+%     peakpos_x(i) = max_x(my);   % x pos of highest peak
+    peakpos_x = [peakpos_x max_x(my)]; 
+    max_y = max_y([max_y<max(max_y)]);   % del max value in y domain
+    max_x = max_x([max_x ~= peakpos_x(end)]); % and in x domain
+end
+peakpos_y = orig_values(peakpos_x);  % extract the y vector 
+
+% --------------- Dynamic Threshold ---------------------
+% works relativ to the mean
+if isfield(params,'dyn_thresh') 
+    if (params.dyn_thresh~=0)
+        % dynamic thresholding 
+        m = mean(orig_values);  
+        dthr = params.dyn_thresh.*m;
+        peakpos_x = peakpos_x([peakpos_y>=dthr]);
+        peakpos_y = peakpos_y([peakpos_y>=dthr]);
+    end
+end
+
+maxima = cell(1, length(peakpos_x));
+% find the left end right minima that belong to a maximum
+for i=1:length(peakpos_x)
+    maxima{i}.x = peakpos_x(i);
+    maxima{i}.t = bin2time(sig_in, maxima{i}.x);
+    maxima{i}.y = orig_values(peakpos_x(i));
+    
+    % find left and right minimum for this maximum
+    maxima{i}.left.x = max(min_x([min_x < maxima{i}.x]));
+    if isempty(maxima{i}.left.x)
+        maxima{i}.left.x = 1;
+        maxima{i}.left.t = 0;
+        maxima{i}.left.y = orig_values(maxima{i}.left.x);
+    else
+        maxima{i}.left.y = orig_values(maxima{i}.left.x);
+        maxima{i}.left.t = bin2time(sig_in, maxima{i}.left.x);
+    end
+    maxima{i}.right.x = min(min_x([min_x > maxima{i}.x]));
+    if isempty(maxima{i}.right.x)
+        maxima{i}.right.x = length(orig_values);
+        maxima{i}.right.t = 0;
+        maxima{i}.right.y = orig_values(maxima{i}.right.x);
+    else
+        maxima{i}.right.y = orig_values(maxima{i}.right.x);
+        maxima{i}.right.t = bin2time(sig_in, maxima{i}.right.x);
+    end
+end
+
+
+out = maxima;
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/analyse_aim_profiles.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,44 @@
+% function  result = analyse_aim_profiles(ti_profile, fq_profile, peaks_tip, peaks_fqp, channel_center_fq)
+%
+%   To analyse the time interval and frequency profile of the auditory image
+%
+%   INPUT VALUES:
+%  
+%   RETURN VALUE:
+%
+% 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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_aim_profiles(ti_profile, fq_profile, peaks_tip, peaks_fqp, channel_center_fq )
+
+
+[fqp_result, fqp_fq] = analyse_frequency_profile(fq_profile, peaks_fqp);
+if fqp_result==1
+    [tip_result, dominant_ti] = analyse_timeinterval_profile(ti_profile, peaks_tip, fqp_result, fqp_fq); 
+else
+    [tip_result, dominant_ti] = analyse_timeinterval_profile(ti_profile, peaks_tip);
+end
+
+% return results
+if fqp_fq==0
+    cfq = 0;
+else
+    cfq =  channel_center_fq(fqp_fq);
+end
+if dominant_ti==0
+    dti = 0;
+else
+    dti = 1/dominant_ti;
+end
+result.tip = tip_result;
+result.dti = dti;
+result.fqp = fqp_result;
+result.dfq = cfq;
+
+
+
+%--------------- subfunction -----------------------
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/analyse_frequency_profile.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,170 @@
+% function [result, center_c] = analyse_frequency_profiletunedDP(fq_profile, peaks, apFQP)
+% 
+%   To analyse the time frequency profile of the auditory image
+%   Returns value between 0 and 1 discribing the strength of the spectral
+%   pitch. Used for quantitative analysis of ramped and damped sinusoids
+%   and for sinusoidally amplitude modulated sinusoids
+%
+%
+%   INPUT VALUES:
+%       fq_profile      frequency profile (signal class)
+%       peaks           output of peakpicker
+%       apFQP           a priori information where to find the peak
+%  
+%   RETURN VALUE:
+%       result      all relevant information 
+%
+% 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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_frequency_profile(fq_profile, options)
+
+% for debug reasons -
+plot_switch = 0;
+
+
+% these finally are the results, that we want back: 
+% 1:the hight at the highest point
+result.free.highest_peak_hight=0;
+result.free.highest_peak_frequency=0;
+% hight to the region below value
+result.free.height_width_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;
+
+% other things useful for plotting
+result.smoothed_signal=fq_profile; % no smoothing used here
+result.peaks = [];
+
+
+% first find the peaks of the frequency profile
+peaks = PeakPicker(fq_profile);
+
+% translate the time in the peaks back to a frequency
+nap=options.nap;    % ugly, but we need the frequency informatin somewhere
+for i=1:length(peaks)
+    peaks{i}.t=1/chan2fre(nap,peaks{i}.x);
+end
+
+result.peaks = peaks;
+
+fq_profile_vals = getdata(fq_profile);
+if plot_switch
+    cf_save = gcf;
+    figure(13)
+    cla;
+    plot(fq_profile_vals,'r-');
+    hold on
+    axis auto
+end
+% stop if there are no peaks, e.g. in the first frames
+if length(peaks)<1
+    return
+end
+
+% Peak of interest is highest peak of the profile
+peaks_oi = peaks{1};
+
+apFQP=options.target_frequency;
+if apFQP>0
+    % no peak finding - peak is given a priori
+    % take nearest peak to a apriori frequency
+    dist=+inf;
+    indexOI = 1;
+    for p=1:length(peaks)
+        d=abs(apFQP-peaks{p}.t);
+        if d<dist
+            indexOI=p;     
+            dist=d;
+        end
+    end %p
+    peak_oi = peaks{indexOI};
+end
+
+
+% %%  Method --- works with SAM sounds
+% %%  Highest peak / number of neighbourgh channels which are bigger than -xdB
+% max_attenuation_dB = -6;
+% atten_fac = 10^(max_attenuation_dB/20);  % attenuation as factor
+% % Take highest Peak
+% %poi = peaks2{1};
+% %poi = peaks{1};
+% poi = peaks_oi;
+% 
+% fq_profile_vals = fq_profile_vals./poi.y;
+% maxy =1;
+% % maxy = poi.y;
+% pwidth = 0;
+% x = poi.x;
+% while (x<=length(fq_profile_vals))&&(fq_profile_vals(x)>maxy*atten_fac)
+%     pwidth=pwidth+1; % one more channel
+%     x=x+1;
+% end
+% if plot_switch
+%   line([x x],[0 peaks_oi.y]);
+% end
+% 
+% x = poi.x-1;
+% while (x>=1)&&(fq_profile_vals(x)>maxy*atten_fac)
+%     pwidth=pwidth+1; % one more channel
+%     x=x-1;
+% end
+% if plot_switch
+%   line([x x],[0 peaks_oi.y]);
+% end
+% result_width =  1-(pwidth/length(fq_profile_vals));
+% 
+
+% ----------------------------
+% Method developed with Roy 13/06
+% Calculate the mean of the Channels in a 20 to 80 percent
+% range left of the main peak
+
+start_frequency_integration=options.start_frequency_integration;
+stop_frequency_integration=options.stop_frequency_integration;
+
+startx=floor(start_frequency_integration*peaks_oi.x);
+if startx<=0
+    startx=1;
+end
+stopx=floor(stop_frequency_integration*peaks_oi.x);
+ps_result = 1 - mean(fq_profile_vals(startx:stopx))/fq_profile_vals(peaks_oi.x);
+
+center_c = peaks_oi.x;
+
+if plot_switch
+  plot(peaks_oi.x,peaks_oi.y,'r.');
+  plot(startx,fq_profile_vals(startx),'bx');
+  line([startx startx],[0 fq_profile_vals(startx)])
+  plot(stopx,fq_profile_vals(stopx),'bx');
+  line([stopx stopx],[0 fq_profile_vals(stopx)])
+  figure(cf_save);
+end
+
+% % % For debug plot minima
+% for i=1:nops
+%     plot(peaks{i}.left.x, peaks{i}.left.y, 'ob');
+%     plot(peaks{i}.right.x, peaks{i}.right.y, 'xb');
+% end
+
+
+% finally define the return values 
+result.free.highest_peak_hight=peaks_oi.y; % height of the first peak
+result.free.highest_peak_frequency=chan2fre(nap,peaks_oi.x);
+
+% hight to the region below value
+result.free.height_width_ratio=ps_result;
+
+% 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;
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/analyse_timeinterval_profile.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,338 @@
+% 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 (bleeck@gmail.com)
+% 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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/displaypitchstrength.asv	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,188 @@
+% generating function for 'aim-mat'
+% 
+%   INPUT VALUES:
+%  
+%   RETURN VALUE:
+%
+% 
+% (c) 2003, University of Cambridge, Medical Research Council 
+% Stefan Bleeck (stefan@bleeck.de)
+% http://www.mrc-cbu.cam.ac.uk/cnbh/aimmanual
+% Christoph Lindner
+% $Date: 2003/07/17 10:56:16 $
+% $Revision: 1.5 $
+
+function displaypitchstrength(sai,options,frame_number,ax)
+if nargin<4
+    ax=gca;
+end
+
+% ?????? per Definition ????
+minimum_time_interval=0.1;  % in ms
+maximum_time_interval=35;
+
+ti_result=options.handles.data.usermodule{frame_number}.ti_result;
+f_result=options.handles.data.usermodule{frame_number}.f_result;
+
+fq_sum = f_result.smoothed_signal;
+int_sum = ti_result.smoothed_signal;
+cla;
+hold on
+
+% TIP
+tip=getvalues(int_sum);
+% the tip is in the form of a logarithmic signal, it must therfore be
+% transformed back to linear:
+
+tip_x = bin2time(int_sum, 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(int_sum,tip_x(1)):time2bin(int_sum,tip_x(end)));
+% tip_x is in ms. Change to Hz
+tip_x = 1./tip_x;
+% plot(tip_x, tip, 'b');
+
+% plot the smoothed interval profile, from which the pitch was derived
+smoothed_signal=ti_result.smoothed_signal;
+% if ~isempty(smoothed_signal)
+% 	smoothed_signal=linsigx(smoothed_signal); % change the logarithmic signal back to linear
+% 	% smoothed_signal=reverse(smoothed_signal);
+% 	smoothval=getvalues(smoothed_signal);
+% 	s_x = bin2time(smoothed_signal, 1:length(smoothval));  % Get the times
+% 	s_x = s_x((s_x>=(minimum_time_interval/1000)) & s_x<=(maximum_time_interval/1000));  
+% 	
+% 	smoothval=smoothval(time2bin(smoothed_signal,s_x(1)):time2bin(smoothed_signal,s_x(end)));
+% %     s_x=f2f(s_x,0,0.035,0.0001,0.035);
+% %     s_x=(s_x+0.0001); %??????????
+% %  	s_x=s_x*1.01; %??????????
+% 	s_x=1./s_x;
+% 	plot(ax,s_x,smoothval,'b');
+% end
+
+% set(gca,'XScale','log');    
+plot(smoothed_signal);
+hold on
+
+% xlabel('Frequency [Hz]');
+set(gca, 'YAxisLocation','right');
+pks = [];
+
+nr_labels = 8;
+ti=50*power(2,[0:nr_labels]);
+for i=1:length(pks)
+    ti = ti((ti>(pks(i)+pks(i)*0.1))|(ti<pks(i)-pks(i)*0.1));
+end
+ti = [ti round(pks)];
+ti = sort(ti);
+% set(gca,'XTick', ti);
+% set(gca, 'XLim',[1000/maximum_time_interval sai{frame_number}.channel_center_fq(end)])
+set(options.handles.checkbox6, 'Value',0);
+set(options.handles.checkbox7, 'Value',0);
+
+
+% if there is a target frequency, then plot this one also!
+target_frequency=options.target_frequency;
+% at least three possibilities: take the hight:
+% pitchstrength=ti_result.free.highest_peak_hight;
+% dominant_frequency_found= ti_result.free.highest_peak_frequency;
+
+
+% three different pitch strengths
+% select the one with the highest peak rate
+
+peaks=ti_result.peaks;
+% select the one with the highest neigbouring_ratio
+peaks2=sortstruct(peaks,'v2012_height_base_width_ratio');
+% peaks2=sortstruct(peaks,'peak_base_height_y');
+% peaks2=sortstruct(peaks,'y');
+
+
+if length(peaks2)>1
+% 	for ii=1:length(peaks)
+	for i=1:3
+		t=peaks2{i}.t;
+		ypeak=peaks2{i}.y;
+% 		ps=peaks{ii}.pitchstrength;
+		ps=peaks2{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(peaks2{i}.left.t,peaks2{i}.left.y,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',5);
+        plot(peaks2{i}.right.t,peaks2{i}.right.y,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',5);
+        
+        
+        ybase=peaks2{i}.v2012_base_where_widths;
+		line([t t],[ybase ypeak],'color','m');
+		line([peaks2{i}.left.t peaks2{i}.right.t],[ybase ybase],'color','m');
+	end
+end
+
+if target_frequency>0 % search only at one special frequency
+	found_fre=ti_result.fixed.highest_peak_frequency;
+	if found_fre>0
+		yval=gettimevalue(smoothed_signal,1/found_fre);
+	else 
+		yval=0;
+	end
+	plot(found_fre,yval,'Marker','o','Markerfacecolor','y','MarkeredgeColor','w','MarkerSize',5);
+	
+	found_ps=ti_result.fixed.highest_peak_hight;
+	plot(found_fre,yval,'Marker','o','Markerfacecolor','g','MarkeredgeColor','w','MarkerSize',10);
+	text(found_fre/1.1,yval*1.01, ['Pitchstrength at fixed ' num2str(round(found_fre)) 'Hz: ' num2str(fround(found_ps,2)) ],'VerticalAlignment','bottom','HorizontalAlignment','right','color','g','Fontsize',12);
+	bordery=get(gca,'Ylim');
+	line([target_frequency target_frequency],[bordery(1) bordery(2)],'color','g')
+    min_fre=target_frequency/options.allowed_frequency_deviation;
+    max_fre=target_frequency*options.allowed_frequency_deviation;
+	line([min_fre min_fre],[bordery(1) bordery(2)],'color','g')
+	line([max_fre max_fre],[bordery(1) bordery(2)],'color','g')
+end
+
+% maxy=sai{end}.ti_resultlt.peaks{1}.y;
+maxy=max(int_sum); 
+if maxy==0
+    maxy=1;
+end
+set(ax,'Ylim',[0,maxy*1.1]);
+
+
+% plot the frequency profile
+plot_frequency_profile=0;
+if plot_frequency_profile
+%Plot both profiles into one figure
+% FQP
+fqp = getvalues(fq_sum)';
+% fqp = fqp /1000;
+plot(sai{frame_number}.channel_center_fq, fqp,'r');
+
+peaks=f_result.peaks;
+if length(peaks)>1
+% 	for ii=1:length(peaks)
+	for ii=1:1
+		fre=1/peaks{ii}.t;
+		yval=peaks{ii}.y;
+% 		ps=peaks{ii}.pitchstrength;
+		ps=peaks{ii}.y;
+		if ii==1
+			plot(fre,yval,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',10);
+			text(fre,yval*1.03,sprintf('%3.0f Hz: %4.2f ',fre,ps),'VerticalAlignment','bottom','HorizontalAlignment','center','color','r','Fontsize',12);
+		else
+			plot(fre,yval,'Marker','o','Markerfacecolor','g','MarkeredgeColor','w','MarkerSize',5);
+			text(fre,yval*1.03,sprintf('%3.0f Hz: %4.2f ',fre,ps),'VerticalAlignment','bottom','HorizontalAlignment','center','color','g','Fontsize',12);
+		end
+	end
+end
+hold off
+end
+
+
+return
+
+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;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/displaypitchstrength.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,203 @@
+% generating function for 'aim-mat'
+% 
+%   INPUT VALUES:
+%  
+%   RETURN VALUE:
+%
+% 
+% (c) 2011, University of Southampton
+% Maintained by Stefan Bleeck (bleeck@gmail.com)
+% 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 displaypitchstrength(sai,options,frame_number,ax)
+if nargin<4
+    ax=gca;
+end
+
+
+% ?????? per Definition ????
+% minimum_time_interval=0.1;  % in ms
+% maximum_time_interval=35;
+
+
+ti_result=options.handles.data.usermodule{frame_number}.ti_result;
+% f_result=options.handles.data.usermodule{frame_number}.f_result;
+if ~isfield(ti_result,'peaks')
+    return;
+end
+
+% fq_sum = f_result.smoothed_signal;
+int_sum = ti_result.smoothed_signal;
+cla;
+hold on
+
+% TIP
+% tip=getvalues(int_sum);
+% the tip is in the form of a logarithmic signal, it must therfore be
+% transformed back to linear:
+
+% tip_x = bin2time(int_sum, 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(int_sum,tip_x(1)):time2bin(int_sum,tip_x(end)));
+% % tip_x is in ms. Change to Hz
+% tip_x = 1./tip_x;
+% plot(tip_x, tip, 'b');
+
+% plot the smoothed interval profile, from which the pitch was derived
+smoothed_signal=ti_result.smoothed_signal;
+% if ~isempty(smoothed_signal)
+% 	smoothed_signal=linsigx(smoothed_signal); % change the logarithmic signal back to linear
+% 	% smoothed_signal=reverse(smoothed_signal);
+% 	smoothval=getvalues(smoothed_signal);
+% 	s_x = bin2time(smoothed_signal, 1:length(smoothval));  % Get the times
+% 	s_x = s_x((s_x>=(minimum_time_interval/1000)) & s_x<=(maximum_time_interval/1000));  
+% 	
+% 	smoothval=smoothval(time2bin(smoothed_signal,s_x(1)):time2bin(smoothed_signal,s_x(end)));
+% %     s_x=f2f(s_x,0,0.035,0.0001,0.035);
+% %     s_x=(s_x+0.0001); %??????????
+% %  	s_x=s_x*1.01; %??????????
+% 	s_x=1./s_x;
+% 	plot(ax,s_x,smoothval,'b');
+% end
+
+% set(gca,'XScale','log');    
+plot(smoothed_signal);
+hold on
+
+% xlabel('Frequency [Hz]');
+set(gca, 'YAxisLocation','right');
+% pks = [];
+
+% nr_labels = 8;
+% ti=50*power(2,[0:nr_labels]);
+% for i=1:length(pks)
+%     ti = ti((ti>(pks(i)+pks(i)*0.1))|(ti<pks(i)-pks(i)*0.1));
+% end
+% ti = [ti round(pks)];
+% ti = sort(ti);
+% set(gca,'XTick', ti);
+% set(gca, 'XLim',[1000/maximum_time_interval sai{frame_number}.channel_center_fq(end)])
+set(options.handles.checkbox6, 'Value',0);
+set(options.handles.checkbox7, 'Value',0);
+
+
+% if there is a target frequency, then plot this one also!
+target_frequency=options.target_frequency;
+% at least three possibilities: take the hight:
+% pitchstrength=ti_result.free.highest_peak_hight;
+% dominant_frequency_found= ti_result.free.highest_peak_frequency;
+
+
+% three different pitch strengths
+% select the one with the highest peak rate
+
+peaks=ti_result.peaks;
+% select the one with the highest neigbouring_ratio
+peaks2=sortstruct(peaks,'v2012_height_base_width_ratio');
+% peaks2=sortstruct(peaks,'peak_base_height_y');
+% peaks2=sortstruct(peaks,'y');
+
+
+if length(peaks2)>1
+% 	for ii=1:length(peaks)
+	for i=1:min(3,length(peaks))
+		t=peaks2{i}.t;
+		ypeak=peaks2{i}.y;
+% 		ps=peaks{ii}.pitchstrength;
+		ps=peaks2{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(peaks2{i}.left.t,peaks2{i}.left.y,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',5);
+        plot(peaks2{i}.right.t,peaks2{i}.right.y,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',5);
+        
+        
+        ybase=peaks2{i}.v2012_base_where_widths;
+		line([t t],[ybase ypeak],'color','m');
+		line([peaks2{i}.left.t peaks2{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')
+end
+
+if target_frequency>0 % search only at one special frequency
+	found_fre=ti_result.fixed.highest_peak_frequency;
+	if found_fre>0
+		yval=gettimevalue(smoothed_signal,1/found_fre);
+	else 
+		yval=0;
+	end
+	plot(found_fre,yval,'Marker','o','Markerfacecolor','y','MarkeredgeColor','w','MarkerSize',5);
+	
+	found_ps=ti_result.fixed.highest_peak_hight;
+	plot(found_fre,yval,'Marker','o','Markerfacecolor','g','MarkeredgeColor','w','MarkerSize',10);
+	text(found_fre/1.1,yval*1.01, ['Pitchstrength at fixed ' num2str(round(found_fre)) 'Hz: ' num2str(fround(found_ps,2)) ],'VerticalAlignment','bottom','HorizontalAlignment','right','color','g','Fontsize',12);
+	bordery=get(gca,'Ylim');
+	line([target_frequency target_frequency],[bordery(1) bordery(2)],'color','g')
+    min_fre=target_frequency/options.allowed_frequency_deviation;
+    max_fre=target_frequency*options.allowed_frequency_deviation;
+	line([min_fre min_fre],[bordery(1) bordery(2)],'color','g')
+	line([max_fre max_fre],[bordery(1) bordery(2)],'color','g')
+end
+
+% maxy=sai{end}.ti_resultlt.peaks{1}.y;
+maxy=max(int_sum); 
+if maxy==0
+    maxy=1;
+end
+set(ax,'Ylim',[0,maxy*1.1]);
+
+
+% plot the frequency profile
+plot_frequency_profile=0;
+if plot_frequency_profile
+%Plot both profiles into one figure
+% FQP
+fqp = getvalues(fq_sum)';
+% fqp = fqp /1000;
+plot(sai{frame_number}.channel_center_fq, fqp,'r');
+
+peaks=f_result.peaks;
+if length(peaks)>1
+% 	for ii=1:length(peaks)
+	for ii=1:1
+		fre=1/peaks{ii}.t;
+		yval=peaks{ii}.y;
+% 		ps=peaks{ii}.pitchstrength;
+		ps=peaks{ii}.y;
+		if ii==1
+			plot(fre,yval,'Marker','o','Markerfacecolor','r','MarkeredgeColor','r','MarkerSize',10);
+			text(fre,yval*1.03,sprintf('%3.0f Hz: %4.2f ',fre,ps),'VerticalAlignment','bottom','HorizontalAlignment','center','color','r','Fontsize',12);
+		else
+			plot(fre,yval,'Marker','o','Markerfacecolor','g','MarkeredgeColor','w','MarkerSize',5);
+			text(fre,yval*1.03,sprintf('%3.0f Hz: %4.2f ',fre,ps),'VerticalAlignment','bottom','HorizontalAlignment','center','color','g','Fontsize',12);
+		end
+	end
+end
+hold off
+end
+
+
+return
+
+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;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/find_pitches.asv	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,472 @@
+% 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 (stefan@bleeck.de)
+% 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;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/find_pitches.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,473 @@
+% 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 (bleeck@gmail.com)
+% 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)
+% 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 =0;
+
+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.02;
+        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;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/find_pitches2.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,301 @@
+% 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 (bleeck@gmail.com)
+% 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;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/genpitchstrength.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,116 @@
+% 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 (bleeck@gmail.com)
+% 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)
+
+nr_frames=length(sai);
+
+target_frequency=options.target_frequency;
+
+% if only one frequency range is important, then find out, which
+% channels are important and which not
+nr_channels=getnrchannels(sai{1});
+
+% find out about scaling:
+maxval=-inf;
+maxfreval=-inf;
+maxsumval=-inf;
+
+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};
+	frame_start_time(frame_number)=getcurrentframestarttime(current_frame);
+	tip(frame_number)=gettimeintervalprofile(current_frame,options);
+	freqp(frame_number)=getfrequencyprofile(current_frame);
+    
+%     figure(234324)
+%     clf
+%     plot(tip(frame_number));
+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);
+
+start_frame=1;
+histo=mute(tip(1));
+% mache ein Histogramm für jede Frequenz mit allem Maxima:
+nr_points=getnrpoints(sai{1});
+
+waithand = waitbar(0,'Generating pitch strength');
+
+for frame_number=start_frame:nr_frames
+% for frame_number=nr_frames:nr_frames
+	
+	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);
+    % Normalisation
+ 	interval_sig =  interval_sig/nr_channels;
+	
+	% sum the interval sum to the floating average:
+	histo=histo+interval_sig;
+	% dynamic decrease of all pitch strengthes
+	histo=histo*ps_decay_constant; 
+% 
+%         figure(234324)
+%     clf
+%     plot(histo);
+    
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	% Peak Picker for IntervalProfile
+	opts.ps_options=options;
+	opts.current_time=getcurrentframestarttime(current_frame);
+	ti_result=find_pitches(histo,opts);
+	ps_result{frame_number}.ti_result=ti_result;
+    ps_result{frame_number}.ti_result.ti_profile=interval_sig;
+    ps_result{frame_number}.ti_result.averaged_signal=histo;
+	
+	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	% Peak Picker for FrequencyProfile
+    frequency_sig=freqp(frame_number);
+    % Normalisation
+    frequency_sig=frequency_sig/getnrpoints(interval_sig);
+    
+ 	ps_result{frame_number}.channel_center_fq = getcf(sai{frame_number});
+
+    % Analyse the dualprofile
+    options.nap=sai{end}; % just for the frequency information
+    f_result= analyse_frequency_profile(frequency_sig,options);
+	ps_result{frame_number}.f_result=f_result;
+
+end
+
+close(waithand);
+
+return
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/genpitchstrength2.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,239 @@
+% 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 (bleeck@gmail.com)
+% 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);
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/genpitchstrength3.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,213 @@
+% 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 (bleeck@gmail.com)
+% 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)
+
+plot_switch=0;
+
+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
+
+
+% 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)=getsum(current_frame);
+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
+
+
+tip_average=mute(tip(1)); % take a copy of the original one
+
+% 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
+	
+	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
+
+
+	% add the profile to the average
+	tip_average=tip_average+tip(frame_number);
+	result=find_pitches(tip_average,opts);
+	tip_average=tip_average*ps_decay_constant; 
+
+	
+	% 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);
+ 		speaks=sortstruct(peaks,'y');
+		
+		pitchstrength=speaks{1}.y/1000000;
+		final_dominant_frequency=1/speaks{1}.t;
+		
+% 						if mod(frame_number,20)==0
+% 							% 		plot(shist);
+% 							% 		ax=axis;
+% 							figure(43242)
+% 							plot(shistl,'r'); hold on
+% 							% 		axis(ax);
+% 							for i=1:1%length(speaks)
+% 								fre=1/speaks{i}.t;
+% 								plot(speaks{i}.x,speaks{i}.y,'o','MarkerFaceColor','g','MarkerEdgeColor','g')
+% 								text(speaks{i}.x,speaks{i}.y,sprintf('%2.2fHz %2.2f',fre,speaks{i}.y))
+% 							end
+% 							time=frame_number*5;
+% 							text(speaks{1}.x+200,speaks{1}.y,sprintf('Time: %3.0fms',time));
+% 							o=0;
+% 						end
+		
+		
+		ps_result{frame_number}.peaks = speaks;
+		ps_result{frame_number}.interval_sum = result.smoothed_signal;
+		
+		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);
+				
+				ps_result{frame_number}.pitchstrength=speaks{best_peak_number}.y/1000;
+				ps_result{frame_number}.dominant_frequency=1/speaks{best_peak_number}.t;
+			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}.interval_sum = result.smoothed_signal;
+		ps_result{frame_number}.peaks = [];
+		ps_result{frame_number}.dominant_frequency=0;
+		ps_result{frame_number}.pitchstrength=0;
+	end
+	
+	
+	% 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);
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/parameters.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,27 @@
+%%%%%%%%%%%%%
+% usermodules
+pitchstrength.generatingfunction='genpitchstrength';
+pitchstrength.displayname='Pitch strength';
+pitchstrength.displayfunction='displaypitchstrength';
+pitchstrength.revision='$Revision: 1.4 $';
+
+
+% The frequency in the temporal profile, that we are interested in
+pitchstrength.target_frequency=0;	% if zero, then look for the highest pitch strength on your own!
+% if we have a target frequency, then the found frequency can deviate by so much:
+pitchstrength.allowed_frequency_deviation=1.12; %two halftones
+
+% how long the temporal profiles are averaged to get the peak
+pitchstrength.tip_average_time=0.05;	% in sec (50 ms)
+pitchstrength.low_pass_frequency=1000;	% Hz is used in the logarithmic(!) tip for averaging
+
+
+% Parameters for the different ways to define the pitch strength
+% Parameters for the qualitly measurement (hight/width)
+pitchstrength.height_width_ratio=0.3;   % at witch point of the height, the width is measured
+
+
+% parameters for the frequency profile
+pitchstrength.start_frequency_integration=0.2;  % start the ratio calculation at a region at 20% of the maximum
+pitchstrength.stop_frequency_integration=0.8;  % stop the ratio calculation at 80% below the maximum
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/usermodule/pitchstrength/reduce_tilt.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,43 @@
+function sig=reduce_tilt(sig,options)
+
+
+plot_switch=1;
+
+lowpassfactor=100;	% relationship between Samplefrequency and lowpassfrequency
+
+% shift all values that are zero at the beginning to the first
+% positive value. This makes the lowpass continuous for high
+% frequencies
+sr=getsr(sig);
+stip=getvalues(sig);
+
+indx=find(stip>0);
+if ~isempty(indx)
+	first_non_zero=stip(indx(1));
+	stip(1:indx(1))=first_non_zero;
+	stip=[ones(200,1)*first_non_zero;stip];
+end
+s=signal(stip,sr);
+s=reverse(s);
+s2=lowpass(s,sr/lowpassfactor);
+diff=s-s2;
+
+% s2=getpart(s2,0,getlength(sig));
+
+
+if plot_switch
+	figure(4)
+	clf
+	plot(s);
+	hold on
+	plot(s2,'r');
+	plot(diff,'g')
+	axis([0 getnrpoints(sig) min(diff)-10 max(sig)+10]);
+end
+
+diff=getpart(diff,0,getlength(sig));
+diff=reverse(diff);
+
+
+% return value
+sig=diff;
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/tools/xticklabel_rotate.m	Tue Aug 16 14:37:17 2011 +0100
@@ -0,0 +1,249 @@
+function hText = xticklabel_rotate(XTick,rot,varargin)
+%hText = xticklabel_rotate(XTick,rot,XTickLabel,varargin)     Rotate XTickLabel
+%
+% Syntax: xticklabel_rotate
+%
+% Input:    
+% {opt}     XTick       - vector array of XTick positions & values (numeric) 
+%                           uses current XTick values or XTickLabel cell array by
+%                           default (if empty) 
+% {opt}     rot         - angle of rotation in degrees, 90° by default
+% {opt}     XTickLabel  - cell array of label strings
+% {opt}     [var]       - "Property-value" pairs passed to text generator
+%                           ex: 'interpreter','none'
+%                               'Color','m','Fontweight','bold'
+%
+% Output:   hText       - handle vector to text labels
+%
+% Example 1:  Rotate existing XTickLabels at their current position by 90°
+%    xticklabel_rotate
+%
+% Example 2:  Rotate existing XTickLabels at their current position by 45° and change
+% font size
+%    xticklabel_rotate([],45,[],'Fontsize',14)
+%
+% Example 3:  Set the positions of the XTicks and rotate them 90°
+%    figure;  plot([1960:2004],randn(45,1)); xlim([1960 2004]);
+%    xticklabel_rotate([1960:2:2004]);
+%
+% Example 4:  Use text labels at XTick positions rotated 45° without tex interpreter
+%    xticklabel_rotate(XTick,45,NameFields,'interpreter','none');
+%
+% Example 5:  Use text labels rotated 90° at current positions
+%    xticklabel_rotate([],90,NameFields);
+%
+% Note : you can not RE-RUN xticklabel_rotate on the same graph. 
+%
+
+
+
+% This is a modified version of xticklabel_rotate90 by Denis Gilbert
+% Modifications include Text labels (in the form of cell array)
+%                       Arbitrary angle rotation
+%                       Output of text handles
+%                       Resizing of axes and title/xlabel/ylabel positions to maintain same overall size 
+%                           and keep text on plot
+%                           (handles small window resizing after, but not well due to proportional placement with 
+%                           fixed font size. To fix this would require a serious resize function)
+%                       Uses current XTick by default
+%                       Uses current XTickLabel is different from XTick values (meaning has been already defined)
+
+% Brian FG Katz
+% bfgkatz@hotmail.com
+% 23-05-03
+% Modified 03-11-06 after user comment
+%	Allow for exisiting XTickLabel cell array
+% Modified 03-03-2006 
+%   Allow for labels top located (after user comment)
+%   Allow case for single XTickLabelName (after user comment)
+%   Reduced the degree of resizing
+% Modified 11-jun-2010
+%   Response to numerous suggestions on MatlabCentral to improve certain
+%   errors.
+
+% Other m-files required: cell2mat
+% Subfunctions: none
+% MAT-files required: none
+%
+% See also: xticklabel_rotate90, TEXT,  SET
+
+% Based on xticklabel_rotate90
+%   Author: Denis Gilbert, Ph.D., physical oceanography
+%   Maurice Lamontagne Institute, Dept. of Fisheries and Oceans Canada
+%   email: gilbertd@dfo-mpo.gc.ca  Web: http://www.qc.dfo-mpo.gc.ca/iml/
+%   February 1998; Last revision: 24-Mar-2003
+
+% check to see if xticklabel_rotate has already been here (no other reason for this to happen)
+if isempty(get(gca,'XTickLabel')),
+    error('xticklabel_rotate : can not process, either xticklabel_rotate has already been run or XTickLabel field has been erased')  ;
+end
+
+% if no XTickLabel AND no XTick are defined use the current XTickLabel
+%if nargin < 3 & (~exist('XTick') | isempty(XTick)),
+% Modified with forum comment by "Nathan Pust" allow the current text labels to be used and property value pairs to be changed for those labels
+if (nargin < 3 || isempty(varargin{1})) & (~exist('XTick') | isempty(XTick)),
+	xTickLabels = get(gca,'XTickLabel')  ; % use current XTickLabel
+	if ~iscell(xTickLabels)
+		% remove trailing spaces if exist (typical with auto generated XTickLabel)
+		temp1 = num2cell(xTickLabels,2)         ;
+		for loop = 1:length(temp1),
+			temp1{loop} = deblank(temp1{loop})  ;
+		end
+		xTickLabels = temp1                     ;
+	end
+varargin = varargin(2:length(varargin));	
+end
+
+% if no XTick is defined use the current XTick
+if (~exist('XTick') | isempty(XTick)),
+    XTick = get(gca,'XTick')        ; % use current XTick 
+end
+
+%Make XTick a column vector
+XTick = XTick(:);
+
+if ~exist('xTickLabels'),
+	% Define the xtickLabels 
+	% If XtickLabel is passed as a cell array then use the text
+	if (length(varargin)>0) & (iscell(varargin{1})),
+        xTickLabels = varargin{1};
+        varargin = varargin(2:length(varargin));
+	else
+        xTickLabels = num2str(XTick);
+	end
+end    
+
+if length(XTick) ~= length(xTickLabels),
+    error('xticklabel_rotate : must have same number of elements in "XTick" and "XTickLabel"')  ;
+end
+
+%Set the Xtick locations and set XTicklabel to an empty string
+set(gca,'XTick',XTick,'XTickLabel','')
+
+if nargin < 2,
+    rot = 90 ;
+end
+
+% Determine the location of the labels based on the position
+% of the xlabel
+hxLabel = get(gca,'XLabel');  % Handle to xlabel
+xLabelString = get(hxLabel,'String');
+
+% if ~isempty(xLabelString)
+%    warning('You may need to manually reset the XLABEL vertical position')
+% end
+
+set(hxLabel,'Units','data');
+xLabelPosition = get(hxLabel,'Position');
+y = xLabelPosition(2);
+
+%CODE below was modified following suggestions from Urs Schwarz
+y=repmat(y,size(XTick,1),1);
+% retrieve current axis' fontsize
+fs = get(gca,'fontsize');
+
+% Place the new xTickLabels by creating TEXT objects
+hText = text(XTick, y, xTickLabels,'fontsize',fs);
+
+% Rotate the text objects by ROT degrees
+%set(hText,'Rotation',rot,'HorizontalAlignment','right',varargin{:})
+% Modified with modified forum comment by "Korey Y" to deal with labels at top
+% Further edits added for axis position
+xAxisLocation = get(gca, 'XAxisLocation');  
+if strcmp(xAxisLocation,'bottom')  
+    set(hText,'Rotation',rot,'HorizontalAlignment','right',varargin{:})  
+else  
+    set(hText,'Rotation',rot,'HorizontalAlignment','left',varargin{:})  
+end
+
+% Adjust the size of the axis to accomodate for longest label (like if they are text ones)
+% This approach keeps the top of the graph at the same place and tries to keep xlabel at the same place
+% This approach keeps the right side of the graph at the same place 
+
+set(get(gca,'xlabel'),'units','data')           ;
+    labxorigpos_data = get(get(gca,'xlabel'),'position')  ;
+set(get(gca,'ylabel'),'units','data')           ;
+    labyorigpos_data = get(get(gca,'ylabel'),'position')  ;
+set(get(gca,'title'),'units','data')           ;
+    labtorigpos_data = get(get(gca,'title'),'position')  ;
+
+set(gca,'units','pixel')                        ;
+set(hText,'units','pixel')                      ;
+set(get(gca,'xlabel'),'units','pixel')          ;
+set(get(gca,'ylabel'),'units','pixel')          ;
+
+origpos = get(gca,'position')                   ;
+
+% textsizes = cell2mat(get(hText,'extent'))       ;
+% Modified with forum comment from "Peter Pan" to deal with case when only one XTickLabelName is given. 
+x = get( hText, 'extent' );  
+if iscell( x ) == true  
+    textsizes = cell2mat( x ) ;  
+else  
+    textsizes = x;  
+end  
+
+largest =  max(textsizes(:,3))                  ;
+longest =  max(textsizes(:,4))                  ;
+
+laborigext = get(get(gca,'xlabel'),'extent')    ;
+laborigpos = get(get(gca,'xlabel'),'position')  ;
+
+labyorigext = get(get(gca,'ylabel'),'extent')   ;
+labyorigpos = get(get(gca,'ylabel'),'position') ;
+leftlabdist = labyorigpos(1) + labyorigext(1)   ;
+
+% assume first entry is the farthest left
+leftpos = get(hText(1),'position')              ;
+leftext = get(hText(1),'extent')                ;
+leftdist = leftpos(1) + leftext(1)              ;
+if leftdist > 0,    leftdist = 0 ; end          % only correct for off screen problems
+
+% botdist = origpos(2) + laborigpos(2)            ;
+% newpos = [origpos(1)-leftdist longest+botdist origpos(3)+leftdist origpos(4)-longest+origpos(2)-botdist]  
+%
+% Modified to allow for top axis labels and to minimize axis resizing
+if strcmp(xAxisLocation,'bottom')  
+    newpos = [origpos(1)-(min(leftdist,labyorigpos(1)))+labyorigpos(1) ...
+            origpos(2)+((longest+laborigpos(2))-get(gca,'FontSize')) ...
+            origpos(3)-(min(leftdist,labyorigpos(1)))+labyorigpos(1)-largest ...
+            origpos(4)-((longest+laborigpos(2))-get(gca,'FontSize'))]  
+else
+    newpos = [origpos(1)-(min(leftdist,labyorigpos(1)))+labyorigpos(1) ...
+            origpos(2) ...
+            origpos(3)-(min(leftdist,labyorigpos(1)))+labyorigpos(1)-largest ...
+            origpos(4)-(longest)+get(gca,'FontSize')]  
+end
+set(gca,'position',newpos)                      ;
+
+% readjust position of text labels after resize of plot
+set(hText,'units','data')                       ;
+for loop= 1:length(hText),
+    set(hText(loop),'position',[XTick(loop), y(loop)])  ;
+end
+
+% adjust position of xlabel and ylabel
+laborigpos = get(get(gca,'xlabel'),'position')  ;
+set(get(gca,'xlabel'),'position',[laborigpos(1) laborigpos(2)-longest 0])   ;
+
+% switch to data coord and fix it all
+set(get(gca,'ylabel'),'units','data')                   ;
+set(get(gca,'ylabel'),'position',labyorigpos_data)      ;
+set(get(gca,'title'),'position',labtorigpos_data)       ;
+
+set(get(gca,'xlabel'),'units','data')                   ;
+    labxorigpos_data_new = get(get(gca,'xlabel'),'position')  ;
+set(get(gca,'xlabel'),'position',[labxorigpos_data(1) labxorigpos_data_new(2)])   ;
+
+
+% Reset all units to normalized to allow future resizing
+set(get(gca,'xlabel'),'units','normalized')          ;
+set(get(gca,'ylabel'),'units','normalized')          ;
+set(get(gca,'title'),'units','normalized')          ;
+set(hText,'units','normalized')                      ;
+set(gca,'units','normalized')                        ;
+
+if nargout < 1,
+    clear hText
+end
+