tomwalters@0: % tomwalters@0: % Dynamic Compressive Gammachirp Filterbank tomwalters@0: % Version 2.05 tomwalters@0: % Toshio IRINO tomwalters@0: % Created: 6 Sep 2003 tomwalters@0: % Modified: 7 Jun 2004 tomwalters@0: % Modified: 12 Jul 2004 (PpgcEstShiftERB) tomwalters@0: % Modified: 14 Jul 2004 (LinPpgc) tomwalters@0: % Modified: 4 Aug 2004 (introducing GCresp) tomwalters@0: % Modified: 16 Aug 2004 (ExpDecayVal) tomwalters@0: % Modified: 31 Aug 2004 (introducing GCFBv2_SetParam) tomwalters@0: % Modified: 8 Sep 2004 (TTS. tidy up the names. 2.00 -> 2.01) tomwalters@0: % Modified: 10 Sep 2004 (Normalization at Level estimation path) tomwalters@0: % Modified: 7 Oct 2004 (c2val is level dependent 2.02) tomwalters@0: % Modified: 22 Oct 2004 (level estimation 2.03) tomwalters@0: % Modified: 8 Nov 2004 (error detection of SndIn) tomwalters@0: % Modified: 30 Nov 2004 (c2val control) tomwalters@0: % Modified: 23 May 2005 (v205. Pc == average of two input, RMS2dBSPL, tomwalters@0: % Fast filtering when 'fix' : under construction) tomwalters@0: % Modified: 24 May 2005 (v205 Mod in LinLvl1 =..., LvldB= ...) tomwalters@0: % Modified: 3 Jun 2005 (v205) tomwalters@0: % Modified: 1 Jun 2005 (v205, GCparam.GainCmpnstdB) tomwalters@0: % Modified: 14 Jul 2005 (v205, GCparam.LvlEst.RefdB, Pwr, Weight) tomwalters@0: % Modified: 15 Sep 2005 (v205, rename GCparam.LvlRefdB --> GainRefdB) tomwalters@0: % tomwalters@0: % tomwalters@0: % function [cGCout, pGCout, Ppgc, GCparam, GCresp] = GCFB2(Snd,GCparam) tomwalters@0: % INPUT: Snd: Input Sound tomwalters@0: % GCparam: Gammachirp parameters tomwalters@0: % GCparam.fs: Sampling rate (48000) tomwalters@0: % GCparam.NumCh: Number of Channels (75) tomwalters@0: % GCparam.FRange: Frequency Range of GCFB [100 6000] tomwalters@0: % specifying asymptotic freq. of passive GC (Fr1) tomwalters@0: % tomwalters@0: % OUTPUT: cGCout: Compressive GammaChirp Filter Output tomwalters@0: % pGCout: Passive GammaChirp Filter Output tomwalters@0: % Ppgc: power at the output of passive GC tomwalters@0: % GCparam: GCparam values tomwalters@0: % GCresp : GC response result tomwalters@0: % tomwalters@0: % Note tomwalters@0: % 1) This version is completely different from GCFB v.1.04 (obsolete). tomwalters@0: % We introduced the "compressive gammachirp" to accomodate both the tomwalters@0: % psychoacoustical simultaneous masking and the compressive tomwalters@0: % characteristics (Irino and Patterson, 2001). The parameters were tomwalters@0: % determined from large dataset (See Patterson, Unoki, and Irino, 2003.) tomwalters@0: % tomwalters@0: % tomwalters@0: % References: tomwalters@0: % Irino, T. and Unoki, M.: IEEE ICASSP'98 paper, AE4.4 (12-15, May, 1998) tomwalters@0: % Irino, T. and Patterson, R.D. : JASA, Vol.101, pp.412-419, 1997. tomwalters@0: % Irino, T. and Patterson, R.D. : JASA, Vol.109, pp.2008-2022, 2001. tomwalters@0: % Patterson, R.D., Unoki, M. and Irino, T. : JASA, Vol.114,pp.1529-1542,2003. tomwalters@0: % tomwalters@0: % tomwalters@0: function [cGCout, pGCout, GCparam, GCresp] = GCFBv2(SndIn,GCparam) tomwalters@0: tomwalters@0: %%%% Handling Input Parameters %%%%% tomwalters@0: if nargin < 2, help GCFBv205; end; tomwalters@0: tomwalters@0: [nc, LenSnd] = size(SndIn); tomwalters@0: if nc ~= 1, error('Check SndIn. It should be 1 ch (Monaural).' ); end; tomwalters@0: tomwalters@0: GCparam = GCFBv205_SetParam(GCparam); tomwalters@0: tomwalters@0: %%%%%%%%%%%%% tomwalters@0: fs = GCparam.fs; tomwalters@0: NumCh = GCparam.NumCh; tomwalters@0: if length(GCparam.b1) == 1 & length(GCparam.c1) == 1 tomwalters@0: b1 = GCparam.b1(1); % Freq. independent tomwalters@0: c1 = GCparam.c1(1); % Freq. independent tomwalters@0: else tomwalters@0: error('Not prepared yet: Freq. dependent b1, c1'); tomwalters@0: end; tomwalters@0: tomwalters@0: [Fr1, ERBrate1] = EqualFreqScale('ERB',NumCh,GCparam.FRange); tomwalters@0: Fr1 = Fr1(:); tomwalters@0: ERBspace1 = mean(diff(ERBrate1)); tomwalters@0: tomwalters@0: GCresp.Fr1 = Fr1; tomwalters@0: Fp1 = Fr2Fpeak(GCparam.n,b1,c1,Fr1); tomwalters@0: GCresp.Fp1 = Fp1; tomwalters@0: tomwalters@0: [ERBrate ERBw] = Freq2ERB(Fr1); tomwalters@0: [ERBrate1kHz ERBw1kHz] = Freq2ERB(1000); tomwalters@0: Ef = ERBrate/ERBrate1kHz - 1; tomwalters@0: GCresp.Ef = Ef; tomwalters@0: %%%% tomwalters@0: % fratVal = frat(1,1) + frat(1,2)*Ef + tomwalters@0: % frat(2,1)*Ppgc + frat(2,2)*Ef*Ppgc; tomwalters@0: %%% tomwalters@0: tomwalters@0: %%%%% Outer-Mid Ear Compensation %%%% tomwalters@0: if length(GCparam.OutMidCrct) > 2 tomwalters@0: % disp(['*** Outer/Middle Ear correction: ' GCparam.OutMidCrct ' ***']); tomwalters@0: CmpnOutMid = OutMidCrctFilt(GCparam.OutMidCrct,fs,0); tomwalters@0: % 1kHz: -4 dB, 2kHz: -1 dB, 4kHz: +4 dB tomwalters@0: Snd = filter(CmpnOutMid,1,SndIn); tomwalters@0: else tomwalters@0: % disp('*** No Outer/Middle Ear correction ***'); tomwalters@0: Snd = SndIn; tomwalters@0: end; tomwalters@0: tomwalters@0: % for compensation filer, use OutMidCrctFilt('ELC',fs,0,1); tomwalters@0: tomwalters@0: %%%%% Gammachirp %%% tomwalters@0: %disp('*** Gammmachirp Calculation ***'); tomwalters@0: waithand=waitbar(0, 'Generating Basilar Membrane Motion - dcGC'); tomwalters@0: if 0, disp(GCparam), end; tomwalters@0: tic; tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: %%%% Start calculation %%%%% tomwalters@0: tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: %%%%% Passive Gammachirp filtering %%%%% tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: tomwalters@0: Tstart = clock; tomwalters@0: cGCout = zeros(NumCh, LenSnd); tomwalters@0: pGCout = zeros(NumCh, LenSnd); tomwalters@0: Ppgc = zeros(NumCh, LenSnd); tomwalters@0: tomwalters@0: for nch=1:NumCh tomwalters@0: tomwalters@0: % passive gammachirp tomwalters@0: pgc = GammaChirp(Fr1(nch),fs,GCparam.n,b1,c1,0,'','peak'); % pGC tomwalters@0: pGCout(nch,1:LenSnd)=fftfilt(pgc,Snd); % fast fft based filtering tomwalters@0: tomwalters@0: end tomwalters@0: tomwalters@0: tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: %%%%% Compressive Gammachirp filtering %%%%% tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: tomwalters@0: %%%%% Initial settings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: tomwalters@0: GCresp.Fr2 = zeros(NumCh,LenSnd); tomwalters@0: GCresp.fratVal = zeros(NumCh,LenSnd); tomwalters@0: tomwalters@0: %%%% Level independent b2 & c2 for Signal path %%% tomwalters@0: b2val = GCparam.b2(1,1)*ones(NumCh,1) + GCparam.b2(1,2)*Ef(:); tomwalters@0: c2val = GCparam.c2(1,1)*ones(NumCh,1) + GCparam.c2(1,2)*Ef(:); tomwalters@0: GCresp.b2val = b2val; tomwalters@0: GCresp.c2val = c2val; tomwalters@0: tomwalters@0: nDisp = 20*fs/1000; % display every 20 ms tomwalters@0: cGCout = zeros(NumCh,LenSnd); tomwalters@0: GCresp.Fr2 = zeros(NumCh,LenSnd); tomwalters@0: GCresp.fratVal = zeros(NumCh,LenSnd); tomwalters@0: LvldB = zeros(NumCh,LenSnd); tomwalters@0: LvlLinPrev = zeros(NumCh,2); tomwalters@0: ExpDecayVal = exp(-1/(GCparam.LvlEst.DecayHL*fs/1000)*log(2)); % decay exp. tomwalters@0: tomwalters@0: NchShift = round(GCparam.LvlEst.LctERB/ERBspace1); tomwalters@0: NchLvlEst = min(max(1, (1:NumCh)'+NchShift),NumCh); % shift in NumCh [1:NumCh] tomwalters@0: Fp1LvlEst = Fp1(NchLvlEst(:)); tomwalters@0: zrs = zeros(NumCh,1); tomwalters@0: LvlLinMinLim = 10^(-GCparam.LvlEst.RMStoSPLdB/20); % minimum should be 0 dBSPL tomwalters@0: LvlLinRef = 10.^(( GCparam.LvlEst.RefdB - GCparam.LvlEst.RMStoSPLdB)/20); tomwalters@0: tomwalters@0: tomwalters@0: %%%%% Sample-by-sample processing %%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: tomwalters@0: Tstart = clock; tomwalters@0: tomwalters@0: %%%%% These lines moved from outside the inner loop tomwalters@0: % TCW August 2006 tomwalters@0: if strcmp(GCparam.Ctrl(1:3),'tim') == 1 tomwalters@0: Fr2LvlEst = GCparam.LvlEst.frat * Fp1LvlEst; tomwalters@0: [ACFcoefLvlEst] = ... tomwalters@0: MakeAsymCmpFiltersV2(fs,Fr2LvlEst,GCparam.LvlEst.b2, GCparam.LvlEst.c2); tomwalters@0: end tomwalters@0: %%%% tomwalters@0: tomwalters@0: tomwalters@0: for nsmpl=1:LenSnd tomwalters@0: tomwalters@0: if strcmp(GCparam.Ctrl,'fix') == 1 tomwalters@0: LvldB(:,nsmpl) = GCparam.GainRefdB*ones(NumCh,1); % fixed value tomwalters@0: tomwalters@0: elseif strcmp(GCparam.Ctrl(1:3),'tim') == 1, % when time-varying tomwalters@0: tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: %%%%%%% Level estimation path %%%%%%%% tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: tomwalters@0: % Moved the following lines outside the main loop as they're not updated tomwalters@0: % each time. Improves speed by about 25% tomwalters@0: % TCW August 2006 tomwalters@0: % Fr2LvlEst = GCparam.LvlEst.frat * Fp1LvlEst; tomwalters@0: % [ACFcoefLvlEst] = ... tomwalters@0: % MakeAsymCmpFiltersV2(fs,Fr2LvlEst,GCparam.LvlEst.b2, GCparam.LvlEst.c2); tomwalters@0: tomwalters@0: tomwalters@0: if nsmpl == 1, %% initialization tomwalters@0: [dummy,ACFstatusLvlEst] = ACFilterBank(ACFcoefLvlEst,[]); tomwalters@0: end; tomwalters@0: [cGCLvlEst,ACFstatusLvlEst] =... tomwalters@0: ACFilterBank(ACFcoefLvlEst,ACFstatusLvlEst,pGCout(NchLvlEst,nsmpl)); tomwalters@0: tomwalters@0: %%%%% Modified: 24 May 05 tomwalters@0: LvlLin(1:NumCh,1) = ... tomwalters@0: max([max(pGCout(NchLvlEst,nsmpl),0), LvlLinPrev(:,1)*ExpDecayVal]')'; tomwalters@0: LvlLin(1:NumCh,2) = ... tomwalters@0: max([max(cGCLvlEst,0), LvlLinPrev(:,2)*ExpDecayVal]')'; tomwalters@0: LvlLinPrev = LvlLin; tomwalters@0: tomwalters@0: %%%%% Modified: 14 July 05 tomwalters@0: LvlLinTtl = GCparam.LvlEst.Weight * ... tomwalters@0: LvlLinRef.*(LvlLin(:,1)/LvlLinRef).^GCparam.LvlEst.Pwr(1) ... tomwalters@0: + ( 1 - GCparam.LvlEst.Weight ) * ... tomwalters@0: LvlLinRef.*(LvlLin(:,2)/LvlLinRef).^GCparam.LvlEst.Pwr(2); tomwalters@0: tomwalters@0: LvldB(:,nsmpl) = 20*log10( max(LvlLinTtl,LvlLinMinLim) ) ... tomwalters@0: + GCparam.LvlEst.RMStoSPLdB; tomwalters@0: else tomwalters@0: error([ 'GCparam.Ctrl should be "fix" or "tim[e-varying]" ']) tomwalters@0: end; tomwalters@0: tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: %%%%%%%%%%% Signal path %%%%%%%%%%%%%% tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tomwalters@0: % Filtering High-Pass Asym. Comp. Filter tomwalters@0: fratVal = GCparam.frat(1,1) + GCparam.frat(1,2)*Ef(:) + ... tomwalters@0: (GCparam.frat(2,1) + GCparam.frat(2,2)*Ef(:)).*LvldB(:,nsmpl); tomwalters@0: Fr2Val = Fp1(:).*fratVal; tomwalters@0: tomwalters@0: [ACFcoef] = MakeAsymCmpFiltersV2(fs,Fr2Val,b2val,c2val); tomwalters@0: if nsmpl == 1, tomwalters@0: [dummy,ACFstatus] = ACFilterBank(ACFcoef,[]); % initiallization tomwalters@0: end; tomwalters@0: tomwalters@0: [SigOut,ACFstatus] = ACFilterBank(ACFcoef,ACFstatus,pGCout(:,nsmpl)); tomwalters@0: cGCout(:,nsmpl) = SigOut; tomwalters@0: tomwalters@0: GCresp.Fr2(:,nsmpl) = Fr2Val; tomwalters@0: GCresp.fratVal(:,nsmpl) = fratVal; tomwalters@0: tomwalters@0: if nsmpl==1 || rem(nsmpl,nDisp)==0, tomwalters@0: waitbar(nsmpl./LenSnd, waithand); tomwalters@0: end % waitbar tomwalters@0: tomwalters@0: end % loop over samples tomwalters@0: tomwalters@0: tomwalters@0: %%%% Signal path Gain Normalization at Reference Level (GainRefdB) %%% tomwalters@0: tomwalters@0: fratRef = GCparam.frat(1,1) + GCparam.frat(1,2)*Ef(:) + ... tomwalters@0: (GCparam.frat(2,1) + GCparam.frat(2,2)*Ef(:)).*GCparam.GainRefdB; tomwalters@0: tomwalters@0: cGCRef = CmprsGCFrsp(Fr1,fs,GCparam.n,b1,c1,fratRef,b2val,c2val); tomwalters@0: GCresp.cGCRef = cGCRef; tomwalters@0: GCresp.LvldB = LvldB; tomwalters@0: tomwalters@0: GainFactor = 10^(GCparam.GainCmpnstdB/20)*(cGCRef.NormFctFp2 * ones(1,LenSnd)); tomwalters@0: cGCout = GainFactor.*cGCout; tomwalters@0: tomwalters@0: %%%%%%%%%% tomwalters@0: close(waithand); tomwalters@0: return