tomwalters@0: %       Inverse Filterbank with Fixed Passive/Compressive Gammachirp
tomwalters@0: %       No dynamic resynthesis
tomwalters@0: %	Version 2.05
tomwalters@0: %       Toshio IRINO
tomwalters@0: %       Created:   9 Nov 2004 
tomwalters@0: %       Modified:  11 Nov 2004  Renamed and based on InvPGCFBv2.m
tomwalters@0: %       Modified:  18 July 2005  v205
tomwalters@0: %
tomwalters@0: %   function [SynSnd, SynGCout] = InvFxGCFB2(GCout,GCparam,GCresp,SwMethod)
tomwalters@0: %      INPUT:  GCout: pGCout or cGCout
tomwalters@0: %              GCparam:  Gammachirp parameters 
tomwalters@0: %              Fp: peak frequency vector for synthesis filter 
tomwalters@0: %              SwMethod: 1: pGCout - pGCsyn
tomwalters@0: %                        2: cGCout - pGCsyn
tomwalters@0: %                        3: cGCout - Fixed(Ref) cGCsyn
tomwalters@0: %        
tomwalters@0: %      OUTPUT: SynSnd:  Synthetic Sound
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 [SynSnd,SynGCout] = InvFxGCFB2(GCout,GCparam,GCresp,SwMethod)
tomwalters@0: 
tomwalters@0: %%%% Handling Input Parameters %%%%%
tomwalters@0: if nargin < 4,         help InvFxGCFBv2;           end;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: GCparam = GCFBv205_SetParam(GCparam);
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%
tomwalters@0: fs = GCparam.fs;
tomwalters@0: NumCh0 = 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: [NumCh, LenSnd] = size(GCout);
tomwalters@0: if NumCh ~= NumCh0, error('Mismatch in NumCh'); end;
tomwalters@0: 
tomwalters@0: if SwMethod == 1 
tomwalters@0:     disp('*** Inverse Passive Gammmachirp Calculation ***');
tomwalters@0:     Fr1 = GCresp.Fr1;
tomwalters@0: elseif SwMethod == 2;
tomwalters@0:     disp('*** Inverse Passive Gammmachirp Calculation ***');
tomwalters@0:     Fr1 = Fpeak2Fr(GCparam.n,b1,c1,GCresp.cGCRef.Fp2); % pGC peak == Fp2
tomwalters@0: elseif  SwMethod == 3
tomwalters@0:     disp('*** Inverse Fixed Compressive Gammmachirp Calculation ***');
tomwalters@0:     Fr1 = GCresp.Fr1;
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: LenInv = 50/1000*fs;  % 50 ms zero filling for processing
tomwalters@0: zz = zeros(NumCh,LenInv);
tomwalters@0: TrGCout = fliplr(GCout);  %%% time-reversal %%
tomwalters@0: TrGCout = [zz TrGCout zz];
tomwalters@0: [NumCh,LenTr] = size(TrGCout);
tomwalters@0: SynGCout = zeros(NumCh,LenTr);
tomwalters@0: TrGCout1 = zeros(NumCh,LenTr);
tomwalters@0: 
tomwalters@0: %%%%% Time-Reversal ACfilter cGC %%%
tomwalters@0: 
tomwalters@0: Tstart = clock;
tomwalters@0: nDisp = 20*fs/1000; % display every 20 ms
tomwalters@0: if SwMethod == 3,  % The same function is applied like wavelet trans.
tomwalters@0:    [ACFcoef] = MakeAsymCmpFiltersV2(fs,...
tomwalters@0:                GCresp.cGCRef.Fr2,GCresp.cGCRef.b2,GCresp.cGCRef.c2);   
tomwalters@0:    [dummy,ACFstatus] =  ACFilterBank(ACFcoef,[]);  % initiallization
tomwalters@0:    for nsmpl = 1:LenTr
tomwalters@0:      [SigOut,ACFstatus] = ACFilterBank(ACFcoef,ACFstatus,TrGCout(:,nsmpl));
tomwalters@0:      TrGCout1(:,nsmpl) = GCresp.cGCRef.NormFctFp2.*SigOut;
tomwalters@0:      if nsmpl==1 | rem(nsmpl,nDisp)==0,
tomwalters@0:        disp(['Time-Reversal ACF-cGC : Time ' num2str(nsmpl/fs*1000,3) ...
tomwalters@0: 	    ' (ms) / ' num2str(LenSnd/fs*1000) ' (ms).    elapsed time = ' ...
tomwalters@0: 	    num2str(fix(etime(clock,Tstart)*10)/10) ' (sec)']);
tomwalters@0:      end;
tomwalters@0:    end;
tomwalters@0:    TrGCout = TrGCout1;
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: %%%%% Time-Reversal Passive Gammachirp  %%%
tomwalters@0: 
tomwalters@0: for nch=1:NumCh
tomwalters@0:     pgc = GammaChirp(Fr1(nch),fs,GCparam.n,b1,c1,0,'','peak');
tomwalters@0:     SynGCout(nch,1:LenTr) =  fftfilt(pgc,TrGCout(nch,:));
tomwalters@0:     if nch == 1 | rem(nch,20)==0  
tomwalters@0:       disp(['Time-Reversal Passive-Gammachirp  ch #' num2str(nch) ...
tomwalters@0: 	    ' / #' num2str(NumCh) '.    elapsed time = ' ...
tomwalters@0: 	    num2str(fix(etime(clock,Tstart)*10)/10) ' (sec)']);
tomwalters@0:     end;
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: SynGCout = fliplr(SynGCout);
tomwalters@0: SynSnd = mean(SynGCout)*sqrt(NumCh); % Approximately all right
tomwalters@0:                                      % Maybe depend on filter shape
tomwalters@0:                                      % 11 Nov 2004
tomwalters@0: 
tomwalters@0: %%%%% Inverse of Outer-Mid Ear Compensation %%%%
tomwalters@0: if length(GCparam.OutMidCrct) > 2
tomwalters@0:  disp(['*** Inverse Outer/Middle Ear correction: ' GCparam.OutMidCrct ' ***']);
tomwalters@0:  InvOutMid = OutMidCrctFilt(GCparam.OutMidCrct,fs,0,1); % Inverse
tomwalters@0:  %% OutMidCrct--> 1kHz: -4 dB, 2kHz: -1 dB, 4kHz: +4 dB
tomwalters@0:  %% [dummy NpIOM] = max(abs(InvOutMid));
tomwalters@0:  SynSnd = filter(InvOutMid,1,SynSnd);  
tomwalters@0:  Nout = LenInv+length(InvOutMid)+(1:LenSnd);
tomwalters@0: else 
tomwalters@0:  disp('*** No Outer/Middle Ear correction ***');
tomwalters@0:  Nout = LenInv+(1:LenSnd);
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: SynSnd   = SynSnd(Nout);
tomwalters@0: SynGCout = SynGCout(:,Nout);
tomwalters@0: 
tomwalters@0: disp(' ');
tomwalters@0: return