Mercurial > hg > aimmat
changeset 4:537f939baef0 tip
various bug fixes and changed copyright message
line wrap: on
line diff
--- /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 +