tomwalters@0: % method of class @frame
tomwalters@0: % 
tomwalters@0: %   INPUT VALUES:
tomwalters@0: %  
tomwalters@0: %   RETURN VALUE:
tomwalters@0: %
tomwalters@0: % 
tomwalters@0: % (c) 2003, University of Cambridge, Medical Research Council 
tomwalters@0: % Stefan Bleeck (stefan@bleeck.de)
tomwalters@0: % http://www.mrc-cbu.cam.ac.uk/cnbh/aimmanual
tomwalters@0: % $Date: 2003/01/17 16:57:46 $
tomwalters@0: % $Revision: 1.3 $
tomwalters@0: 
tomwalters@0: function ret=frequencysharpening(current_frame,options)
tomwalters@0: % frequency sharpening does local lateral inhibition according to its
tomwalters@0: % parameters
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if nargin < 2
tomwalters@0:     options=[];
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if isfield(options,'excitatory_area') % how many erbs the excitation lasts
tomwalters@0:     excitatory_area=options.excitatory_area;
tomwalters@0: else
tomwalters@0:     excitatory_area=0.5;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if isfield(options,'grafix') % how many erbs the excitation lasts
tomwalters@0:     grafix=options.grafix;
tomwalters@0: else
tomwalters@0:     grafix=0;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if isfield(options,'inhibitory_area') % how many erbs the inhibition lasts
tomwalters@0:     inhibitory_area=options.inhibitory_area;
tomwalters@0: else
tomwalters@0:     inhibitory_area=1.5;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if isfield(options,'inhibitory_influence') % how strong the inhibition is
tomwalters@0:     inhibitory_influence=options.inhibitory_influence;
tomwalters@0: else
tomwalters@0:     inhibitory_influence=0.2;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: cfs=getcf(current_frame);
tomwalters@0: EarQ = 9.26449;				%  Glasberg and Moore Parameters
tomwalters@0: minBW = 24.7;
tomwalters@0: order = 1;
tomwalters@0: ERB = ((cfs/EarQ).^order + minBW^order).^(1/order);
tomwalters@0: B=1.019*2*pi.*ERB;
tomwalters@0: 
tomwalters@0: vals=getvalues(current_frame);
tomwalters@0: nr_chan=size(vals,1);
tomwalters@0: nr_dots=size(vals,2);
tomwalters@0: new_vals=zeros(size(vals));
tomwalters@0: 
tomwalters@0: % calculate the interchannel crossinfluence for each channel
tomwalters@0: erb1=21.4*log10(4.73e-3*cfs(1)+1);  % from Glasberg,Moore 1990
tomwalters@0: erb2=21.4*log10(4.73e-3*cfs(nr_chan)+1);
tomwalters@0: nr_erbs=erb2-erb1;
tomwalters@0: erb_density=nr_chan/nr_erbs;
tomwalters@0: neighbourinfluence=zeros(nr_chan);
tomwalters@0: if grafix
tomwalters@0:     figure(2)
tomwalters@0:     clf
tomwalters@0: end
tomwalters@0: for i=1:nr_chan  % through all channels: prepare working variable
tomwalters@0:     count=1;
tomwalters@0:     for j=i-nr_chan:i+nr_chan  % calculate the influence of all channels
tomwalters@0:         if j>0 & j< nr_chan 
tomwalters@0:             dist=abs(i-j);  % we are symmetric
tomwalters@0:             disterbs=dist/erb_density;  % so many erbs between the channel and me
tomwalters@0:             
tomwalters@0:             % the function(disterbs) to describe the influence of
tomwalters@0:             % neighbouring channels:
tomwalters@0:             %             final=1/disterbs;
tomwalters@0:             
tomwalters@0:             % Mexican hat function:
tomwalters@0:             %             p=disterbs*disterbs/(threshold_decay_halflife*threshold_decay_halflife);
tomwalters@0:             % version with power
tomwalters@0:             %             final=(2-p)*exp(-0.5*p);
tomwalters@0:             % version with difference from gaussians:
tomwalters@0:             p1=disterbs*disterbs/(excitatory_area*excitatory_area);
tomwalters@0:             p2=disterbs*disterbs/(inhibitory_area*inhibitory_area);
tomwalters@0:             final=exp(-p1)-inhibitory_influence*exp(-p2);
tomwalters@0:             
tomwalters@0:             neighbourinfluence(i,count)=final;
tomwalters@0:             count=count+1;
tomwalters@0:         end
tomwalters@0:     end
tomwalters@0:     % normalize it
tomwalters@0:     %     scale=1/sum(neighbourinfluence(i,:));
tomwalters@0:     %     neighbourinfluence(i,:)=neighbourinfluence(i,:)*scale;
tomwalters@0:     if grafix
tomwalters@0:         figure(2)
tomwalters@0:         plot(neighbourinfluence(i,:)); hold on
tomwalters@0:         if i==30
tomwalters@0:             plot(neighbourinfluence(i,:),'.-r','linewidth',3); hold on
tomwalters@0:         end
tomwalters@0:     end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: for i=1:nr_dots % through the time
tomwalters@0:     for j=1:nr_chan  % through all channels: prepare working variable
tomwalters@0:         new_vals(:,i)=new_vals(:,i)+vals(:,i).*neighbourinfluence(:,j);
tomwalters@0:     end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: s=getfrequencysum(current_frame);
tomwalters@0: current_frame=setvalues(current_frame,new_vals);
tomwalters@0: 
tomwalters@0: ret=current_frame;
tomwalters@0: 
tomwalters@0: if grafix
tomwalters@0:     figure(3) ;clf;
tomwalters@0:     s2=getfrequencysum(current_frame);
tomwalters@0:     plot(s); hold on
tomwalters@0:     ax=axis;
tomwalters@0:     plot(s2,'.-r');
tomwalters@0:     axis([0 58 0 40000]);
tomwalters@0: end