tomwalters@0: %
tomwalters@0: % 	Calculation of NAP 
tomwalters@0: %	(gammatone, halfwave rectification, lowpass, log compression)
tomwalters@0: % 	7 Jan 2002
tomwalters@0: %	Toshio Irino
tomwalters@0: %
tomwalters@0: function  [NAP, NAPparam, BMM] = CalNAPghll(Snd, NAPparam)
tomwalters@0: 
tomwalters@0: fs = NAPparam.fs;
tomwalters@0: if isfield(NAPparam,'NumCh') == 0,  NAPparam.NumCh = [];   end;
tomwalters@0: if length(NAPparam.NumCh) == 0, 
tomwalters@0: 	NAPparam.NumCh  = 75;
tomwalters@0: end;
tomwalters@0: if isfield(NAPparam,'cf_afb') == 0,  NAPparam.cf_afb = [];   end;
tomwalters@0: if length(NAPparam.cf_afb) == 0, 
tomwalters@0: 	NAPparam.cf_afb  = [100 6000];
tomwalters@0: end;
tomwalters@0: if isfield(NAPparam,'SubBase') == 0,  NAPparam.SubBase = [];   end;
tomwalters@0: if length(NAPparam.SubBase) == 0, 
tomwalters@0: 	NAPparam.SubBase  = 0;
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: fs = NAPparam.fs;
tomwalters@0: NumCh = NAPparam.NumCh;
tomwalters@0: cf_afb = NAPparam.cf_afb;
tomwalters@0: 
tomwalters@0: %%%%% Outer-Mid Ear Compensation %%%%
tomwalters@0: CmpnOutMid = OutMidCrctFilt('ELC',fs,0);
tomwalters@0: Snd_om = filter(CmpnOutMid,1,Snd);
tomwalters@0: % Snd_om = Snd; % if unnecessary
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%%% BMM  %%%
tomwalters@0: disp('*** BMM & NAP Calculation ***');
tomwalters@0: disp(NAPparam)
tomwalters@0: tic;
tomwalters@0: %% dERB = (Freq2ERB(max(cf_afb))-Freq2ERB(min(cf_afb)))/(NumCh-1);
tomwalters@0: %% Frs =  ERB2Freq( Freq2ERB(min(cf_afb)):dERB:Freq2ERB(max(cf_afb)) );
tomwalters@0: Frs = FcNch2EqERB(min(cf_afb),max(cf_afb),NumCh);
tomwalters@0: NAPparam.Frs = Frs;
tomwalters@0: NAPparam.b = 1.019; % default gammatone
tomwalters@0: 
tomwalters@0: % IIR implementation
tomwalters@0: % disp('BMM : Start calculation, Wait a minute');
tomwalters@0: % fcoefs = MakeERBFilters98B(fs,Frs,[],b);   % new version
tomwalters@0: % BMM = ERBFilterBank(Snd_om, fcoefs);
tomwalters@0: % disp(['BMM : elapsed_time = ' num2str(toc,3) ' (sec)']);
tomwalters@0: 
tomwalters@0: %%%%% Lowpass filter for representing Phase-lock property %%%
tomwalters@0: flpcut = 1200;
tomwalters@0: [bzLP apLP] = butter(1,flpcut/(fs/2));
tomwalters@0: bzLP2 = [bzLP(1)^2,  2*bzLP(1)*bzLP(2), bzLP(2)^2]; 
tomwalters@0: apLP2 = [apLP(1)^2,  2*apLP(1)*apLP(2), apLP(2)^2]; 
tomwalters@0: 
tomwalters@0: %%%%% NAP %%%%
tomwalters@0: tic;
tomwalters@0: LenSnd = length(Snd_om);
tomwalters@0: bias = 0.1;
tomwalters@0: for nch = 1:NumCh
tomwalters@0:      if rem(nch, 20) == 0 | nch == 1 |  nch == NumCh 
tomwalters@0: 	 disp(['BMM-NAP #' int2str(nch) '/#' int2str(NumCh) ':  ' ...
tomwalters@0: 		 'elapsed_time = ' num2str(toc,3) ' (sec)']);
tomwalters@0:      end;
tomwalters@0: 
tomwalters@0:      gt = GammaChirp(Frs(nch),fs,4,NAPparam.b,0,0,[],'peak');
tomwalters@0:      %  BMM(nch,:) = filter(gt,1,Snd_om);
tomwalters@0:      BMM(nch,:) = fftfilt(gt,Snd_om);
tomwalters@0:      NAPraw = log10(max(BMM(nch,:),bias)) - log10(bias);
tomwalters@0:      NAP(nch,1:LenSnd) = filter(bzLP2,apLP2,NAPraw); % LP filtered 
tomwalters@0: 
tomwalters@0: 	% [Frs(nch) max(max(NAPraw)) max(max(NAP(nch,:)))]
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if NAPparam.SubBase ~= 0,
tomwalters@0: disp([ '=== Baseline subtraction :  Max NAP = ' num2str( max(max(NAP))) ... 
tomwalters@0: 	', NAPparam.SubBase = ' num2str(NAPparam.SubBase) ' ===']);
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: NAP0 = NAP;
tomwalters@0: NAP = max((NAP0 - NAPparam.SubBase),0);
tomwalters@0: 
tomwalters@0: NAPparam.height = max(max(NAP));
tomwalters@0: NAPparam.tms    = (0:LenSnd-1)/fs*1000;
tomwalters@0: 
tomwalters@0: return
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: if 1
tomwalters@0:      ProfNAP0 = mean(NAP0')';
tomwalters@0:      ProfNAP = mean(NAP')';
tomwalters@0:      plot(NAPparam.Frs,ProfNAP0,NAPparam.Frs,ProfNAP0)
tomwalters@0:      subplot(2,1,1)
tomwalters@0:      image(flipud(NAP0(:,1:100:LenSnd))*15)
tomwalters@0:      subplot(2,1,2)
tomwalters@0:      image(flipud(NAP(:,1:100:LenSnd))*20)
tomwalters@0: end;     
tomwalters@0: 
tomwalters@0: 
tomwalters@0: return
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: % [bzGT, apGT, Frs, ERBw] = MakeERBFiltersB(fs,NumCh,Frs,b); % NG for low freq
tomwalters@0: 
tomwalters@0: %%Simple pre-emphasis
tomwalters@0: %%CoefPreEmph = [1, -0.97];
tomwalters@0: %%Snd = filter(CoefPreEmph,1,Snd);
tomwalters@0: 
tomwalters@0: % [frsp, freq] = freqz(bzLP,apLP,1024,fs);
tomwalters@0: % [frsp2, freq] = freqz(bzLP2,apLP2,1024,fs);
tomwalters@0: % plot(freq,20*log10(abs(frsp)),freq,20*log10(abs(frsp2)))
tomwalters@0: % axis([0 2000 -10 2]);
tomwalters@0: % grid;