function varargout = baumgartner2016_spectralanalysis(sig,spl,varargin)
%BAUMGARTNER2016_SPECTRALANALYSIS Spectral analysis
% Usage: [mp,fc] = baumgartner2016_spectralanalysis(sig,spl)
%
% Input parameters:
% sig : incoming time-domain signal
% spl : sound pressure level (re 20e-6 Pa) in dB. Set to NAN in
% order to bypass internal level adjustment.
% type : flag for target (default) or template
% name : identifying string for caching (e.g., 'NH12_baseline')
%
% Output parameters:
% mp : spectral magintude profile. Dimensions (4-6 optional):
% 1) frequency, 2) position (polar angle), 3) channel, 4)
% fiber type, 5) time frame.
% fc : center frequencies of auditory filters
%
% BAUMGARTNER2016_SPECTRALANALYSIS(...) computes temporally integrated
% spectral magnitude profiles.
%
% BAUMGARTNER2016_SPECTRALANALYSIS accepts the following optional parameters:
%
% 'tiwin' Set temporal integration window in seconds. Default is Inf.
% 'ID' Listener's ID (important for caching).
% 'Condition' Label of experimental condition (also for caching).
%
% BAUMGARTNER2016_SPECTRALANALYSIS accepts these optional flags:
%
% 'target' Processing of a target sound (for caching). This is the
% default.
% 'template' Processing of a template sound (for caching).
% 'gammatone' To apply a Gammatone filterbank instead of the Zilany et
% al. (2014) model of the auditory periphery.
%
% Url: http://amtoolbox.org/amt-1.5.0/doc/modelstages/baumgartner2016_spectralanalysis.php
% #StatusDoc: Perfect
% #StatusCode: Perfect
% #Verification: Verified
% #Requirements: SOFA M-Signal M-Stats O-Statistics
% #Author: Robert Baumgartner (2016), Acoustics Research Institute, Vienna, Austria
% This file is licensed unter the GNU General Public License (GPL) either
% version 3 of the license, or any later version as published by the Free Software
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and
% at <https://www.gnu.org/licenses/gpl-3.0.html>.
% You can redistribute this file and/or modify it under the terms of the GPLv3.
% This file is distributed without any warranty; without even the implied warranty
% of merchantability or fitness for a particular purpose.
definput.import={'baumgartner2016'}; % defaults from arg_baumgartner2016
definput.flags.type={'target','template'};
[flags,kv]=ltfatarghelper({},definput,varargin);
% Set cachename
if flags.do_target
cachenameprefix = 'ireptar_';
if size(sig,1) > kv.tiwin*kv.fs; cachenameprefix = [cachenameprefix 'tiwin' num2str(kv.tiwin*1e3) 'ms_']; end
if not(isempty(kv.ID))
cachenameprefix = [cachenameprefix kv.ID '_'];
end
if not(isempty(kv.Condition))
cachenameprefix = [cachenameprefix kv.Condition '_'];
end
%%% equalized conditions (don't run periphery model again on it)
cachenameprefix = strrep(cachenameprefix,'BB','baseline');
cachenameprefix = strrep(cachenameprefix,'CL','baseline');
%%%
else
cachenameprefix = 'ireptem_';
if not(isempty(kv.ID))
cachenameprefix = [cachenameprefix kv.ID '_'];
end
end
cachenameprefix = [cachenameprefix 'lat' num2str(kv.lat) '_' num2str(spl) 'dB'];
%% Remove pausings (at beginning and end)
Nch = size(sig,3);
if not(isnan(spl))
idnz = mean(sig(:,:).^2,2) ~=0 | [diff(mean(sig(:,:).^2,2)) ~= 0;1];
sig = sig(idnz,:,:);
% and evaluate broadband ILD
if Nch == 2
ILD = dbspl(sig(:,:,2))./dbspl(sig(:,:,1));
if mean(ILD) < 1
chdamp = 2;
else
ILD = 1./ILD;
chdamp = 1;
end
end
end
%% Gammatone
if flags.do_gammatone
cachename = [cachenameprefix '_gammatone_' num2str(1/kv.space,'%u') 'bpERB'];
if flags.do_middleear; cachename = [cachename '_middleear']; end
if flags.do_ihc; cachename = [cachename '_ihc']; end
[mp,fc] = amt_cache('get',cachename,flags.cachemode);
if isempty(mp)
% Set level by single factor, not individually for every direction or
% channel (directions and channels are pooled)
if not(isnan(spl))
sigsize = size(sig);
sig = scaletodbspl(sig(:),spl, 100);
sig = reshape(sig,sigsize);
end
if flags.do_middleear
miearfilt = middleearfilter(kv.fs);
sig = lconv(sig,miearfilt(:));
end
Lpad = round(5/kv.flow*kv.fs); % add 5 cycles of lowest freq
sig = postpad(sig,size(sig,1)+Lpad);
if kv.space == 1
[mp,fc] = auditoryfilterbank(sig(:,:),kv.fs,...
'flow',kv.flow,'fhigh',kv.fhigh);
else
fc = audspacebw(kv.flow,kv.fhigh,kv.space,'erb');
[bgt,agt] = gammatone(fc,kv.fs,'complex');
mp = 2*real(ufilterbankz(bgt,agt,sig(:,:))); % channel (3rd) dimension resolved!
end
Nfc = length(fc); % # bands
% IHC transduction
if flags.do_ihc
mp = ihcenvelope(mp,kv.fs,'ihc_dau1996');
end
% % Set back the channel dimension
% mp = reshape(mp,[size(mp,1),Nfc,size(sig,2),size(sig,3)]);
% Averaging over time (RMS)
if flags.do_target && size(mp,1) > kv.tiwin*kv.fs
Lframe = round(kv.tiwin*kv.fs); % length of each frame
Nframes = ceil(size(mp,1)/Lframe); % # frames
mp = postpad(mp,Lframe*Nframes,0,1);
mp = reshape(mp,[Lframe,Nframes,Nfc,size(sig,2),Nch]);
mp = permute(mp,[1,3,4,5,6,2]); % frames as last dimension
mp = shiftdim(rms(mp));
time = (0:Lframe:Lframe*Nframes-1)/kv.fs;
else % integrate over whole duration
mp = reshape(mp,[size(mp,1),Nfc,size(sig,2),Nch]); % Set back the channel dimension
mp = shiftdim(rms(mp));
time = 0;
end
% Logarithmic transformation (dB)
mp = 100 + 20*log10(mp);
amt_cache('set',cachename,mp,fc);
end
% Limit dynamic range
if not(isnan(spl))
mp = min(mp,kv.GT_maxSPL); % maybe +30 because dynamic range was evaluated with broadband noise (40 auditory bands)
mp = max(mp,kv.GT_minSPL);
end
end
%% Zilany Model
if flags.do_zilany2007 || flags.do_zilany2014
ftd = [0.16,0.23,0.61]; % Liberman (1978)
ftweights = ftd(kv.fiberTypes)/sum(ftd(kv.fiberTypes));
for tt = 1:length(kv.fiberTypes)
ft = kv.fiberTypes(tt);
cachename = [cachenameprefix '_fiberType' num2str(ft)];
if kv.cihc < 1; cachename = [cachename '_cihc' num2str(kv.cihc)]; end
if kv.cohc < 1; cachename = [cachename '_cohc' num2str(kv.cohc)]; end
try
[mp,fc,time] = amt_cache('get',cachename,flags.cachemode);
catch
[mp,fc] = amt_cache('get',cachename,flags.cachemode);
time = 0:kv.tiwin:(size(mp,5)-1)*kv.tiwin;
end
if isempty(mp)
Nmin = .05*kv.fsmod; % decay time at 700-Hz in response to click
if length(sig)/kv.fs*kv.fsmod < Nmin
sig = postpad(sig,ceil(Nmin*kv.fs/kv.fsmod),0,1);
end
amt_disp(['Compute: ' cachename]);
Ntar = size(sig,2); % # target angles
len = ceil(length(sig)/kv.fs*kv.fsmod);
ANresp = zeros(len,kv.numCF,Ntar,2);
for ch = 1:Nch
for ii = 1:Ntar
if Nch > 1 && ch == chdamp;
spl_mod = ILD(ii)*spl;
else
spl_mod = spl;
end
if flags.do_zilany2007
[ANout,fc] = zilany2007(spl_mod,sig(:,ii,ch),kv.fs,...
kv.fsmod,'flow',kv.flow,'fhigh',kv.fhigh,'numCF',kv.numCF);
ANout = ANout';
else % zilany2014
fc = audspace(kv.flow,kv.fhigh,kv.numCF,'erb');
ANout = zilany2014(scaletodbspl(sig(:,ii,ch),spl_mod),kv.fs,...
fc,'fiberType',ft,'cohc',kv.cohc,'cihc',kv.cihc);
ANout=ANout';
end
% Compensate for cochlear delay?!?
% Check stimulus onset
ionset = find(diff(mean(ANout.^2,1)) ~= 0,1,'first');
ANresp(:,:,ii,ch) = ANout(:,(1:len)+ionset-1)';
amt_disp([num2str(ii+(ch-1)*Ntar) ' of ' num2str(Ntar*Nch) ' done'],'volatile');
end
end
amt_disp();
% Averaging over time (RMS)
if flags.do_target && size(ANresp,1) > kv.tiwin*kv.fsmod
Lframe = kv.tiwin*kv.fsmod; % length of each frame
Nframes = ceil(size(ANresp,1)/Lframe); % # frames
ANresp = postpad(ANresp,Lframe*Nframes,0,1);
ANresp = reshape(ANresp,[Lframe,Nframes,length(fc),size(sig,2),size(sig,3)]);
ANresp = permute(ANresp,[1,3,4,5,6,2]); % frames as last dimension
time = (0:Lframe:Lframe*Nframes-1)/kv.fsmod;
else % integrate over whole duration
% ANresp = reshape(ANresp,[size(len,length(fc),size(sig,2),size(sig,3)]); % retreive polar dimension if squeezed out
time = 0;
end
mp = shiftdim(mean(ANresp),1);
amt_cache('set',cachename,mp,fc,time);
end
% if size(mp,2) ~= size(sig,2) % retreive polar dimension if squeezed out
% mp = reshape(mp,[size(mp,1),size(sig,2),size(sig,3)]);
% end
if tt == 1 % init
mp_cum = zeros(kv.numCF,size(sig,2),size(sig,3),1,length(time));
mp_sep = zeros(kv.numCF,size(sig,2),size(sig,3),length(kv.fiberTypes),length(time));
end
mp_cum = mp_cum + ftweights(tt)*mp;
mp_sep(:,:,:,tt,:) = mp;
end
if flags.do_ftcum
mp = mp_cum;
else % flags.do_ftopt
mp = mp_sep;
end
% adjust frequency range for cached data
idf = fc <= kv.fhigh & fc >= kv.flow;
fc = fc(idf);
mp = mp(idf,:,:,:,:);
end
% fiber activity gating
% mp = round(mp*0.3);
varargout{1} = mp;
if nargout > 1
varargout{2} = fc;
if nargout > 2
varargout{2} = time;
end
end
end