The Auditory Modeling Toolbox
Go to function
KELVASA2015 - Localization model in cochlear-implant listeners (Kelvasa and Dietz, 2015)
Program code:
function [results] = kelvasa2015(insig,fs,varargin)
%Kelvasa and Dietz 2015 Bilateral CI localization model
% Usage: [results] = kelvasa2015(insig,fs);
%
% KELVASA2015(insig,fs) implements the ACE signal processing strategy
% upon the two channel input signal to produce bilateral electrodograms.
% This is further processed through an electrode nerve interface to generate
% spike times of a population of AN neurons. A chosen localization model from Kelvasa
% and Dietz 2015 is then used to map the two channel (right and left) outputs
% to a predicted azimuthal position.
%
% Input parameters:
% insig : Can be either: [N x 2] two channel audio signal
% or results structure of preProcessed data
%
% fs : sampling rate (Hz)
%
% varargin : structure with all parameters required for model. If
% this is not included, default paramters are loaded.
%
% Output parameters:
% results : A structure containing the processed electrodograms,
% AN spike times, and model predicted azimuthal locations
%
% The output structure "results" has the following fields:
%
% electrodogramCHAN1 : [NxM] matrix of CI electrode current output
% in mA(???) with N = number of CI
% electrodes and M = time
%
% APvecCHAN1 : [Nx2] matrix of [Nx1] indices of spiking AN
% fibers [Nx2] spike times in seconds
%
% electrodogramCHAN2 : same but for second channel
%
% APvecCHAN2 : same but for second channel
%
% SpkSumPerBin : [NxM] matrix of Right and Left Spike Rate
% differences in spikes per second with
% N= number of AN frequency bands and
% M = time bins
%
% SpkDiffPerBin : same but for Right and Left spike rate differences
%
% ANbinPredictions : [NxM] matrix of azimuthal angle bin predictions in
% degrees with N= number of AN frequency
% bands and M = time bins
%
% weightedPredictions : [1xM] matrix of bin weighted azimuthal angle
% bin predictions in degrees with
% M = time bins
%
% mappingData : Structure containing data used to calibrate and
% implement the chosen localization
% model as detailed in
% kelvasa2015calibratemodels.m
%
% The steps of the binaural model to calculate the result are the
% following :
%
% 1) Process two channel input signal through a CI strategy as detailed
% in (Hamacher, 2003) and (Fredelake and Hohmann, 2012)
% to produce bilateral electrodograms.
%
% 2) Process electrodogram through an electrode nerve interface and
% auditory nerve model as detailed in (Fredelake and
% Hohmann, 2012)
%
% 3) Compute bilateral spike rate differences over chosen AN frequency
% bands and time windows as detailed in (Kelvasa and
% Dietz, 2015)
%
% 4) Calibrate the chosen localization model with a chosen calibration
% signal. This step can take several hours so
% preProcessed calibration is loaded for "Speech
% Shaped Noise" at 55dB as detailed in (Kelvasa and
% Dietz, 2015)
%
% 5) Map the spike rate differences for each AN frequency band to a
% predicted azimuthal angle using the chosen
% localization model as detailed in (Kelvasa and
% Dietz, 2015)
%
% Parameters implemented in the model processing stages are set through
% AMT using the following files: arg_kelvasa2015ciparams.m,
% arg_kelvasa2015anparams.m, and arg_kelvasa2015locationmodelparams.m in
% which thorough descriptions of the input parameters are given.
%
% Kelvasa, D., & Dietz, M. (2015). Auditory model-based sound direction
% estimation with bilateral cochlear implants. Trends in hearing, 19.
%
% Fredelake, S., & Hohmann, V. (2012). Factors affecting predicted speech
% intelligibility with cochlear implants in an auditory model for
% electrical stimulation. Hearing research, 287(1), 76-90.
%
% Hamacher, V. (2003). PhD Thesis. "Signalverarbeitungsmodelle des
% elektrisch stimulierten Geh�rs."
%
% References
% 1. http://www.sciencedirect.com/science/article/pii/S016763931000097X
%
% Url: http://amtoolbox.sourceforge.net/doc/experiments/exp_kelvasa2015.php
%
% Copyright (C) 2009-2016 Piotr Majdak and Peter L. Søndergaard.
% This file is part of AMToolbox version 0.9.7
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Authors:
% Daryl Kelvasa (daryl.kelvasa@uni-oldenburg.de) 2016
% Mathias Dietz (mdietz@uwo.ca) 2016
%
% See also: arg_kelvasa2015ciparams.m, arg_kelvasa2015anparams.m,
% arg_kelvasa2015locationmodelparams.m, kelvasa2015anprocessing.m,
% kelvasa2015ciprocessing.m, kelvasa2015anbinning.m
% kelvasa2015calibratemodels.m, kelvasa2015localize.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Retrieve and compute model paramters
%
% Url: http://amtoolbox.sourceforge.net/data/amt-test/htdocs/amt-0.9.8/doc/models/kelvasa2015.php
% Copyright (C) 2009-2015 Piotr Majdak and Peter L. Søndergaard.
% This file is part of AMToolbox version 0.9.8
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% Set flags
definput.flags.debugging = {'debug','no_debug'};
definput.flags.plot = {'no_plot','plot'};
definput.flags.plot_stage_fig = {'plot_stage_fig','no_plot_stage_fig',};
% Import default arguments from other functions
definput.import={'kelvasa2015','amtcache'};
[flags,kv] = ltfatarghelper({},definput,varargin);
% Check Inputs
if nargin<2
error('%s: Too few input parameters.',upper(mfilename));
end;
if isstruct(insig);
results = insig; doPreProcessing = 0;
else doPreProcessing = 1;
end
if isnumeric(insig) ; if min(size(insig))~=2
error('%s: insig has to be a numeric two channel signal!',upper(mfilename));
end; end
if ~isnumeric(fs) || ~isscalar(fs) || fs<=0
error('%s: fs has to be a positive scalar!',upper(mfilename));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model processing starts here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if doPreProcessing
%Signal resampling to CI processing
if kv.FS_ACE ~= fs
insig = resample(insig,kv.FS_ACE,fs);
debug('Resampling stimulus to ACE sampling frequency', flags); end
results.sigLengthSec = (size(insig,1)/kv.FS_ACE);
results.signal = insig;
results.fs = kv.FS_ACE;
results.modelParameters = kv;
for chan = 1 : 2
signal = insig(:,chan);
%CI Model adapted from (Hamacher, 2003) and
%(Fredelake and Hohmann, 2012)
debug('Running ACE signal processing strategy',...
flags);
[electrodogram, vTime] = ...
kelvasa2015ciprocessing(signal,...
kv.FS_ACE,'argimport',flags,kv);
%CI/Auditory Nerve interface and AN model adapted from
%(Fredelake and Hohmann, 2012)
debug('Running Electrode/Nerve interface and AN model',...
flags);
[APvec] = ...
kelvasa2015anprocessing(electrodogram,...
vTime,'argimport',flags,kv);
if chan == 1
results.vTime = vTime;
results.electrodogramCHAN1 = electrodogram;
results.APvecCHAN1 = APvec; else
results.electrodogramCHAN2 = electrodogram;
results.APvecCHAN2 = APvec; end
end
else debug('Using preProcessed AN and CI data from struct input',...
flags);
end
%Processing starts here if preProcessing has already been performed
%AN frequency band binning of AN spike rate outputs as
%described in Kelvasa and Dietz 2015
debug('Binning AN outputs into AN frequency bands along cochlea',...
flags);
[~,spikeRatePerBin1] = ...
kelvasa2015anbinning(results.APvecCHAN1,...
results.sigLengthSec,'argimport',flags,kv);
[~,spikeRatePerBin2] = ...
kelvasa2015anbinning(results.APvecCHAN2,...
results.sigLengthSec,'argimport',flags,kv);
results.SpkSumPerBin = spikeRatePerBin2+spikeRatePerBin1;
results.SpkDiffPerBin = spikeRatePerBin2-spikeRatePerBin1;
%Calibration of the linear and statistical localization models as
%described in Kelvasa and Dietz 2015.
calibDataFilename = [kv.localizationModelCalibWav(1:end-4),'_',...
num2str(kv.localizationModelCalibStimulusLevelDB),...
'dB_calibration_', kv.identifier];
%Check for preprocessed calibration data
[mappingData] = amtcache('get', calibDataFilename, flags.cachemode);
if isempty(mappingData)
debug('Recomputing localization model calibration', flags);
[mappingData] = kelvasa2015calibratemapping('argimport',...
flags,kv);
amtcache('set',calibDataFilename,mappingData);
end
%Mapping of Right-Left spike rate for each bin to a predicted
%azimuthal angle as described in Kelvasa and Dietz 2015
debug(['Implementing ',kv.localizationModel,...
' localization model to predict azimuthal source locations.'],...
flags);
[ANbinPredictions,...
weightedPredictions,...
mappingData] = ...
kelvasa2015localize(mappingData, ...
results.SpkDiffPerBin,...
results.SpkSumPerBin,...
'argimport',...
flags,kv);
results.ANbinPredictions = ANbinPredictions;
results.weightedPredictions = weightedPredictions;
results.mappingData = mappingData;
if flags.do_plot_stage_fig
plot_kelvasa2015mainFigure(results);
end
end
%%
function debug(text_str,flags)
if flags.do_debug
amtdisp(text_str); end
end
%%
function plot_kelvasa2015mainFigure(results)
%% Plot main results
%Generate Main Figure
mainFig = figure('units','centimeters',...
'position',[0.5 0.5 21 17.5]);
rows = 4;
columns = 2;
xIndent = 3;
yIndent = 1.5;
ySpace = 1;
axisWidth = 8;
axisHeight = 3;
xTix = linspace(0,results.sigLengthSec,4);
xTixLab = [num2str(round(100.*(xTix'))./100),repmat(' s',numel(xTix),1)];
for axisInd = 2 : rows*columns
%Generate axis
mapInd = [0 1 2 3 0 1 2 3];
cornerX = xIndent + axisWidth*(ceil(axisInd/rows)-1) + ...
ceil(axisInd/rows)-1;
cornerY = yIndent + mapInd(axisInd)*(axisHeight+ySpace);
if axisInd == 5; cornerX = cornerX - axisWidth- 1; end
ax(axisInd) = axes('Parent',mainFig,...
'Units','centimeters',...
'Fontsize', 14,...
'XTickLabel',[],...
'YTickLabel',[],...
'Position',...
[cornerX cornerY axisWidth axisHeight]);
hold(ax(axisInd),'on')
axis(ax(axisInd),'manual') ;
%Plot data
switch axisInd
case {2,6}
if axisInd == 2;
data = results.APvecCHAN1;
set(ax(axisInd),'YTick',[1,500],...
'YTickLabel',round([1,500].*(35/500)));
ylabel({'Distance in','cochlea from Apex,','(mm)'});
else data = results.APvecCHAN2;
end
a = scatter(data(:,2),data(:,1));
set(a, 'MarkerFaceColor','black',...
'MarkerEdgeColor','black',...
'SizeData',1);
set(ax(axisInd),'ylim',[1 500],...
'xlim',[0 max(results.sigLengthSec)],...
'XTick',xTix,...
'XTickLabel', xTixLab);
case {3,7}
elecCHAN1 = results.electrodogramCHAN1;
elecCHAN2 = results.electrodogramCHAN2;
elecRange = [min(min(elecCHAN1(:)), min(elecCHAN2(:))),...
max(max(elecCHAN1(:)), max(elecCHAN2(:)))];
if axisInd == 3;
data = elecCHAN1;
ylabel({'CI Electrode','Current,','(\mu A)'});
set(ax(axisInd),...
'YTick',[1,6,11,17,22],...
'YTickLabel',[1,6,11,17,22]);
else data = elecCHAN2;
end
imagesc(results.vTime,[],data,elecRange);
set(ax(axisInd),'ylim',[1 22],...
'xlim',[0 max(results.vTime)],...
'XTickLabel',[])
case {4,8}
WINDOW = round(0.012 * results.fs);
NOVERLAP = round(0.5*WINDOW); NFFT = 10000;
[S1,~,~] = spectrogram(results.signal(:,1),WINDOW,...
NOVERLAP,NFFT,results.fs);
[S2 F T] = spectrogram(results.signal(:,2),WINDOW,...
NOVERLAP,NFFT,results.fs); S = [S1;S2];
colorRange = [min(20.*log10(abs(S(:)))),...
max(20.*log10(abs(S(:))))];
if axisInd == 4;
data = S1;
title('Chan 1 (Left)')
ylabel({'Frequency,','(kHz)',''});
set(ax(axisInd),...
'YTick',0:2000:8000,...
'YTickLabel',(0:2000:8000)./1000);
else data = S2;
title('Chan 2 (Right)')
end
imagesc(T,F,20.*log10(abs(data)),colorRange)
set(ax(axisInd),'ylim',[0 results.fs/2],...
'xlim',[0 max(T)],...
'XTickLabel',[])
case 5
%Set colors for each bin plot
col = hsv(numel(results.modelParameters.binPos));
col(2,:) = col(2,:) + [0.2 -0.5 0];
timeWin = results.modelParameters.timeWindowSec;
timeBins = (1:numel(results.weightedPredictions)).*timeWin;
numBins = numel(results.modelParameters.binPos);
for bin = 1 : numBins + 1
if ~isempty(results.ANbinPredictions) && bin <= numBins
a = scatter(timeBins,...
results.ANbinPredictions(bin,:),100);
set(a,'MarkerFaceColor',col(bin,:),...
'Marker','o');
else
a = scatter(timeBins,...
results.weightedPredictions,100);
set(a,'MarkerFaceColor','black',...
'Marker','o');
end
end
set(ax(axisInd),'ylim',[-10 100],...
'xlim',[0 results.sigLengthSec+timeWin],...
'YTickLabel',[],...
'YTick',[0 30 60 90],...
'XTick',timeBins,...
'YTickLabel',[0 30 60 90],...
'YGrid','on',...
'XGrid','on')
ylabel({'Predicted','Azimuthal','Angle, (�)'})
xlabel(['Time Bins, ',num2str(timeWin),'s'])
leg = cellstr([num2str(results.modelParameters.binPos'),...
repmat(' mm',numBins,1)]);
leg{end+1} = 'weighted';
a = legend(leg); b = get(a,'position'); b(1) = b(1)+4.5;
b(2) = b(2)+0.25; set(a,'Position',b)
end
end
end
Build with Bootstrap