tomwalters@0: % Generating function for AIM-MAT
tomwalters@0: % Based on Dick Lyon's code for the Pole-Zero filter cascade
bleeck@3: % (c) 2011, University of Southampton
bleeck@3: % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@3: % download of current version is on the soundsoftware site: 
bleeck@3: % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@3: % documentation and everything is on http://www.acousticscale.org
tomwalters@0: function fr=gen_pzfc(sig,options)
tomwalters@0: 
tomwalters@0: % Unpack the options and get the raw signal from the signal class
tomwalters@0: CFMin=options.lowest_frequency;
tomwalters@0: CFMax=options.highest_frequency;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % Default parameters (these defaults should be the same as in the parameter
tomwalters@0: % file)
tomwalters@0: if ~isfield(options,'segment_length')
tomwalters@0:     segLength=200;
tomwalters@0: else
tomwalters@0:     segLength=options.segment_length;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options,'stepfactor')
tomwalters@0:     stepfactor=1/3;
tomwalters@0: else
tomwalters@0:     stepfactor=options.stepfactor;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options,'pdamp')
tomwalters@0:     pdamp=0.12;
tomwalters@0: else
tomwalters@0:     pdamp=options.pdamp;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options,'zdamp')
tomwalters@0:     zdamp=0.2;
tomwalters@0: else
tomwalters@0:     zdamp=options.zdamp;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options,'zfactor')
tomwalters@0:     zfactor=1.4;
tomwalters@0: else
tomwalters@0:     zfactor=options.zfactor;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'bw_over_cf')
tomwalters@0:     bwOverCF = 0.11;
tomwalters@0: else
tomwalters@0:     bwOverCF = options.bw_over_cf;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'bw_min_hz');
tomwalters@0:     bwMinHz = 27; 
tomwalters@0: else
tomwalters@0:     bwMinHz=options.bw_min_hz;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'use_fitted_params')
tomwalters@0:     useFit=0;
tomwalters@0: else
tomwalters@0:     useFit=options.use_fitted_params;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'agcfactor')
tomwalters@0:     agcfactor=12;
tomwalters@0: else
tomwalters@0:     agcfactor=options.agcfactor;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'do_agc')
tomwalters@0:     doagcstep=1;
tomwalters@0: else
tomwalters@0:     doagcstep=options.do_agc;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'pre_excite_with_noise')
tomwalters@0:     do_pre_excite=0;
tomwalters@0: else
tomwalters@0:     do_pre_excite=options.pre_excite_with_noise;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'pre_excite_time')
tomwalters@0:     pre_excite_time=0.2; % 200ms
tomwalters@0: else
tomwalters@0:     pre_excite_time=options.pre_excite_time;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'pre_excite_level_dB')
tomwalters@0:     pre_excite_level_dB=-10;
tomwalters@0: else
tomwalters@0:     pre_excite_level_dB=options.pre_excite_level_dB;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if ~isfield(options, 'doplot')
tomwalters@0:     doplot=0;
tomwalters@0: else
tomwalters@0:     doplot=options.doplot;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if useFit~=0 
tomwalters@0:    if isfield(options, 'fitted_params')
tomwalters@0:         eval(options.fitted_params); % sets the varaible ValParam
tomwalters@0:         if ~isvector('ValParam')
tomwalters@0:             error('PZFC param file specified is not in the right format. Should set variable ValParam.');
tomwalters@0:         end
tomwalters@0:    else
tomwalters@0:         error('If use_fitted_params is nonzero (i.e. true) then fitted_params must be a real filename');
tomwalters@0:    end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: fs=getsr(sig);
tomwalters@0: inputColumn=getvalues(sig);
tomwalters@0: Nt = length(inputColumn);
tomwalters@0: 
tomwalters@0: if useFit==0
tomwalters@0:     if options.erb_scale == 0
tomwalters@0:         [Coeffs Nch, freqs] = PZBankCoeffsHiLow(fs, CFMax, CFMin, pdamp, zdamp, zfactor, stepfactor, bwOverCF, bwMinHz);
tomwalters@0:     else
tomwalters@0:         [Coeffs Nch, freqs] = PZBankCoeffsERB(fs, CFMax, CFMin, pdamp, zdamp, zfactor, stepfactor);
tomwalters@0:     end
tomwalters@0: else
tomwalters@0:     % Use the parameters from the fit to set up the filterbank (using ERB scale)
tomwalters@0:     [Coeffs,  Nch, freqs] = PZBankCoeffsERBFitted(fs, CFMax, CFMin, ValParam);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: deskewDelay = 0;
tomwalters@0: doStereo = 0; % no stereo option for AIM-MAT version
tomwalters@0: Ntracks=1; % ditto - don't try and make this any more than 1!
tomwalters@0: 
tomwalters@0: pfreqs = Coeffs(:,1);
tomwalters@0: pdamps = Coeffs(:,2);
tomwalters@0: za0 = Coeffs(:,3);
tomwalters@0: za1 = Coeffs(:,4);
tomwalters@0: za2 = Coeffs(:,5);
tomwalters@0: Nch = length(pfreqs);
tomwalters@0: 
tomwalters@0: mindamp = .18; % (Q = 3 max end of interp)
tomwalters@0: maxdamp = .4; % (Q = 1.25 min end of interp)
tomwalters@0: % min and max refer to the damping, not the x and r coordinates
tomwalters@0: rmin = exp(-mindamp*pfreqs);
tomwalters@0: rmax = exp(-maxdamp*pfreqs);
tomwalters@0: xmin = rmin.*cos(pfreqs.*(1-mindamp^2).^0.5);
tomwalters@0: xmax = rmax.*cos(pfreqs.*(1-maxdamp^2).^0.5);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: agcepsilons = [0.0064, 0.0016, 0.0004, 0.0001]; % roughly 3 ms, 12 ms, 48 ms, 200 ms
tomwalters@0: agcgains = [1, 1.4, 2, 2.8];
tomwalters@0: agcgains = agcgains/mean(agcgains);
tomwalters@0: Nstages = length(agcepsilons);
tomwalters@0: 
tomwalters@0: outarray = []; %zeros(Nch*Ntracks,Nt);
tomwalters@0: for j = 1:Ntracks
tomwalters@0:     inputColumn(:,j) = filter([0.5 0.5], 1, inputColumn(:,j)); % puts a zero in the TF
tomwalters@0: end
tomwalters@0: 
tomwalters@0: state1 = zeros(Nch,Ntracks);
tomwalters@0: state2 = zeros(Nch,Ntracks);
tomwalters@0: prevout = zeros(Nch,Ntracks);
tomwalters@0: agcstate = zeros(Nch, Nstages*Ntracks);
tomwalters@0: % initialize with double dampling:
tomwalters@0: [pdampsmod agcstate] = AGCdampStep([ ], pdamps, agcepsilons, agcgains, agcstate); % initialize pdampsmod, agcstate
tomwalters@0: % make safer place to start than full gain, in case input is abruptly loud:
tomwalters@0: pdampsmod = pdampsmod + 0.05;
tomwalters@0: agcstate = agcstate + 0.05;
tomwalters@0: 
tomwalters@0: % Noise pre-excitation
tomwalters@0: if do_pre_excite == 1
tomwalters@0:     wbh=waitbar(0, 'Pre-calculating noise excitaiton for Pole-Zero filter cascade'); 
tomwalters@0:     Nt_noise=pre_excite_time.*fs;
tomwalters@0:     doagcstep_noise=1; % we want to run the AGC in this section, even if not in the main section.
tomwalters@0:     
tomwalters@0:     nsegs = floor(Nt_noise/segLength);
tomwalters@0:     %wbh=waitbar(0, 'Calculating Pole-Zero filter cascade');
tomwalters@0:     for segno = 1:nsegs
tomwalters@0:         inputSegment = rand(segLength,1);
tomwalters@0:         inputSegment = inputSegment.*10.^(pre_excite_level_dB./20);
tomwalters@0: 
tomwalters@0:         [outsegment, prevout, agcstate, state1, state2, pdampsmod] = ...
tomwalters@0:             RunPZBankSegment(inputSegment, prevout, agcstate, state1, state2, pfreqs, pdamps, pdampsmod, ...
tomwalters@0:             mindamp, maxdamp, xmin, xmax, rmin, rmax, za0, za1, za2, agcepsilons, agcgains, agcfactor, doagcstep_noise, doplot);
tomwalters@0: 
tomwalters@0:         waitbar(segno/nsegs, wbh);
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     end
tomwalters@0:     close(wbh); 
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % Reset the previous output of the filterbank
tomwalters@0: state1 = zeros(Nch,Ntracks);
tomwalters@0: state2 = zeros(Nch,Ntracks);
tomwalters@0: prevout = zeros(Nch,Ntracks);
tomwalters@0: 
tomwalters@0: nsegs = floor(Nt/segLength);
tomwalters@0: wbh=waitbar(0, 'Calculating Pole-Zero filter cascade');
tomwalters@0: for segno = 1:nsegs
tomwalters@0:     inputSegment = inputColumn((segno-1)*segLength + (1:segLength), :);
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     [outsegment, prevout, agcstate, state1, state2, pdampsmod] = ...
tomwalters@0:         RunPZBankSegment(inputSegment, prevout, agcstate, state1, state2, pfreqs, pdamps, pdampsmod, ...
tomwalters@0:         mindamp, maxdamp, xmin, xmax, rmin, rmax, za0, za1, za2, agcepsilons, agcgains, agcfactor, doagcstep, doplot);
tomwalters@0: 
tomwalters@0:     if isempty(outarray)
tomwalters@0:         outarray=outsegment;
tomwalters@0:     else
tomwalters@0:         outarray = [outarray, outsegment]; % gather up the segments
tomwalters@0:     end
tomwalters@0: 
tomwalters@0:     waitbar(segno/nsegs, wbh);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: % Convert all the stuff into the form that AIM-MAT wants it
tomwalters@0: outarray=flipdim(outarray, 1);
tomwalters@0: freqs=flipdim(freqs, 1);
tomwalters@0: 
tomwalters@0: fr=frame(outarray);
tomwalters@0: fr=setsr(fr,fs);
tomwalters@0: fr=setstarttime(fr,getminimumtime(sig));
tomwalters@0: fr=setcf(fr,freqs);
tomwalters@0: close(wbh);
tomwalters@0: 
tomwalters@0: return;
tomwalters@0: