tomwalters@0: %   CalMICoef
tomwalters@0: %   SAI ->  Mellin Image Coefficient Direct
tomwalters@0: %   IRINO T.
tomwalters@0: %   18 Jan 01
tomwalters@0: %   27 Jun 01 (modified to MFCC type, Not on LogFrq)
tomwalters@0: %   11 Jan 02 (NAPparam, MIparam)
tomwalters@0: %
tomwalters@0: %   Modified for the size shape image
tomwalters@0: %   Marc A. Al-Hames
tomwalters@0: %   April 2003
tomwalters@0: %
tomwalters@0: %   function [ MICoef ] = CalMellinCoef(SAIPhsCmp,NAPparam,MIparam)
tomwalters@0: %   INPUT:  SAIPhsCmp : SAI val with Phase Compensation
tomwalters@0: %	    NAPparam:
tomwalters@0: %	           fs       : Sampling Frequency
tomwalters@0: %	           Frs      : Channel frequencies
tomwalters@0: %	    MIparam:
tomwalters@0: %	           RangeAudFig:   Range of Auditory Figure 
tomwalters@0: %                         [ZERO, Boundary] in sampling-point
tomwalters@0: %	           TFval    : TFval   == Hval          (--> abscissa of MI)
tomwalters@0: %	           c_2pi    : Kernel spatial frequeny  (--> ordinate of MI)
tomwalters@0: %	           Mu       : Kernel spatial weighting
tomwalters@0: %   OUTPUT: MICoef   : MI value
tomwalters@0: 
tomwalters@0: function [ ssi ] = Calssicoef(SAIPhsCmp,NAPparam,MIparam)
tomwalters@0: 
tomwalters@0: fs    = NAPparam.fs;
tomwalters@0: Frs   = NAPparam.Frs;
tomwalters@0: RangeAudFig =  MIparam.RangeAudFig;
tomwalters@0: TFval = MIparam.TFval;
tomwalters@0: c_2pi = MIparam.c_2pi;
tomwalters@0: Mu    = MIparam.Mu;
tomwalters@0: 
tomwalters@0: [NumCh LenSAI] = size(SAIPhsCmp);
tomwalters@0: LenAF =  diff(RangeAudFig)+1;
tomwalters@0: LenTaper = round(0.5*NAPparam.fs/1000); % 0.5 ms taper
tomwalters@0: WinAF = TaperWindow(LenAF+LenTaper,'han',LenTaper);
tomwalters@0: WinAF = ones(NumCh,1)*WinAF(LenTaper+(1:LenAF)); 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: AFval = WinAF .* SAIPhsCmp(:,RangeAudFig(1):RangeAudFig(2));
tomwalters@0: [NumCh,LenAF] = size(AFval);
tomwalters@0: %MICoef=AFval; %added line
tomwalters@0: %%%%%%%%%%
tomwalters@0: %% LogFrs = log10(Frs(:)/min(Frs))/log10(6000/100); % normalized in [100 6000]
tomwalters@0: %% NormFrq = LogFrs;
tomwalters@0: %% Change to MFCC type, DCT on ERB domain  on 27 Jun 2001
tomwalters@0: NormFreq = ( (0:NumCh-1) + 0.5 )/NumCh;
tomwalters@0: c_pi = 2*c_2pi;
tomwalters@0: 
tomwalters@0: amp = exp((Mu-0.5)*0.5)*sqrt(2/NumCh); 	% mag. norm. when NormFreq == 0.5
tomwalters@0: 					% when using ERB --> Frs~=760 Hz
tomwalters@0: Kernel = amp*exp( ( i*pi*c_pi(:) - (Mu-0.5)) * NormFreq(:)');
tomwalters@0: Kernel(1,1:NumCh) = Kernel(1,1:NumCh)/sqrt(2);
tomwalters@0: %
tomwalters@0: % for confirmation (15 Jan 2002)
tomwalters@0: % Kernel1 = Kernel;
tomwalters@0: % Kernel2 = DCTWarpFreq(0,length(NormFreq),length(c_pi),0);
tomwalters@0: % plot(real(Kernel1(7,:)))
tomwalters@0: % hold on
tomwalters@0: % plot(1:NumCh,real(Kernel(4,:)),'r--', 1:NumCh,abs(Kernel(4,:)))
tomwalters@0: % sum(Kernel1(3,:)-Kernel2(2,:))
tomwalters@0: % hold off
tomwalters@0: % pause
tomwalters@0: %
tomwalters@0: % 	Kernel 	Mag at 100Hz	Mag at 6000Hz
tomwalters@0: % Mu = 2:	2.1170    	0.4724  % lowpass
tomwalters@0: % Mu = 1:	1.2840    	0.7788  % lowpass
tomwalters@0: % Mu = 0.5:	1.0		1.0	% flat
tomwalters@0: % Mu = 0:	0.7788    	1.2840  % high pass
tomwalters@0: %
tomwalters@0: % clf
tomwalters@0: % plot(LogFrs,Kernel); 
tomwalters@0: % amp*exp((-Mu+0.5)*[0 1])
tomwalters@0: 
tomwalters@0: %%%%%%%%%%
tomwalters@0: 
tomwalters@0: TFmargin = 0.1;
tomwalters@0: AFave = zeros(NumCh,length(TFval));
tomwalters@0: MICoef = zeros(length(c_2pi),length(TFval));
tomwalters@0: 
tomwalters@0: ValNorm = 1;
tomwalters@0: 
tomwalters@0: % It's worth notin what's going on here, becuase it's not completely
tomwalters@0: % obvious. AFval contains the 
tomwalters@0: for cntTF = 1:length(TFval);
tomwalters@0:    cntCh = 0;
tomwalters@0:    for nch = 1:NumCh
tomwalters@0: 	nSAImin = max(1,fix((TFval(cntTF)-TFmargin)/Frs(nch)*fs));
tomwalters@0: 	nSAImax = min(ceil((TFval(cntTF)+TFmargin)/Frs(nch)*fs), LenAF);
tomwalters@0: 	% aaa(nch,1:5) = [cntTF, nch, Frs(nch), nSAImin nSAImax];
tomwalters@0: 	if nSAImin <= nSAImax,
tomwalters@0: 	  AFave(nch,cntTF) = mean(AFval(nch,nSAImin:nSAImax))/ValNorm;
tomwalters@0: 	  cntCh = cntCh +1;
tomwalters@0:         end;
tomwalters@0:    end;
tomwalters@0:    % NumValidCh(cntTF) = cntCh;
tomwalters@0: %   ChNorm = NumCh/cntCh;  % 10 Oct 01 --> Too much normalization 
tomwalters@0:    % ChNorm = 1;	 	  % it seems the best at 15 Jan 02   
tomwalters@0:    % MICoef(1:length(c_2pi),cntTF) = ( Kernel*AFave(:,cntTF) )*ChNorm;
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: ssi = AFave;
tomwalters@0: 
tomwalters@0: % Normalization by Number of channels within one AF,
tomwalters@0: % which depends on h values.
tomwalters@0: 
tomwalters@0: % ValRatio = sum(NumValidCh)/(length(TFval)*NumCh);
tomwalters@0: % MICoef = MICoef*ValRatio;
tomwalters@0: 
tomwalters@0: return
tomwalters@0: 
tomwalters@0: % Do not apply this. (12 Mar. 2001)
tomwalters@0: % if length(SAIval) > 1,
tomwalters@0: %	ValNorm = max(mean(SAIval));  % Mag. of strobing point
tomwalters@0: %	if ValNorm < 0.01, ValNorm = 1; end;
tomwalters@0: % end;