Skip to content

GitLab

Commit bba7b439 authored by dualberger's avatar dualberger
Browse files

add all required files from JHF original folders

parent ad2c0a66
Hide whitespace changes
Inline Side-by-side
Showing with 2227 additions and 0 deletions
BAMQ/BAMQ.m 0 → 100644
View file @ bba7b439
function binQ = BAMQ(mRef, mTest, fs)
%BAMQ calculate perceptual binaural difference between mRef and mTest
% -------------------------------------------------------------------------
%
% Usage: binQ = BAMQ(mRef, mTest, fs)
%
% Input: mRef ... reference signal matrix (Nx2)
% mTest ... aligned test signal matrix (Nx2)
% fs ... sampling frequency (files are resampled to 44.1 kHz)
%
% 0 dB FS should correspond to 115 dB SPL
%
% Output: binQ ... binaural perceptual quality value
%
% Values: 100 ... no difference
% 0 ... obvious difference
% -X ... even bigger difference
%
%
% Copyright (C) University of Oldenburg and HrTech 2012-2015
% Author : flja
% Date : 12-02-2019
funcFile = mfilename('fullpath');
pathstr = fileparts(funcFile);
addpath([pathstr filesep 'functions'])
if nargin < 3,
error('Not enough input arguments.');
end
if fs ~= 44100
mRef = resample(mRef,44100,fs);
mTest = resample(mTest,44100,fs);
fs = 44100;
end
if size(mRef,1) < fs*0.4,
error('signals need to be at least 0.4 s long')
end
% avoid NaNs in IVS calculation:
mRef = [0.00001,0.00001;mRef(1:end-1,:)];
mTest = [0.00001,0.00001;mTest(1:end-1,:)];
%% FRONT END: modified Dietz model
[mITDref, mIVSref, mILDref, mPowRef1, mPowRef2] = BAMQFrontEnd(mRef,fs);
[mITDtest, mIVStest, mILDtest, mPowTest1, mPowTest2, fc, fsNew] = BAMQFrontEnd(mTest,fs);
%% BACK END
binQ = BAMQBackEnd(mITDref, mIVSref, mILDref, mPowRef1, mPowRef2, ...
mITDtest, mIVStest, mILDtest, mPowTest1, mPowTest2, fc, fsNew);
binQ = ceil(binQ);
BAMQ/functions/BAMQBackEnd.m 0 → 100644
View file @ bba7b439
function binQ = BAMQBackEnd(mITDref, mIVSref, mILDref, mPowRef1, mPowRef2, ...
mITDtest, mIVStest, mILDtest, mPowTest1, mPowTest2, fc, fsNew)
%% Initialization
ivsThres = 0.98;
% parameters for equation 14 & 15
k.r = 2;
k.s = 0.0012;
k.t = 666;
k.a = 0.023;
n.r = 3.17;
n.s = -0.0015;
n.t = 560;
n.a = -2.75;
% initialize block computation
blockLenFac = 0.4;
blockLen = round(blockLenFac.*fsNew);
nShifts = floor(size(mITDref,1)/blockLen);
nResiSamples = mod(size(mITDref,1),blockLen);
if nResiSamples ~= 0,
nShifts = nShifts + 1;
end
vILDdiffs = [];
vITDdiffs = [];
vIVSdiffs = [];
%% Computation
for iShift = 1:nShifts,
% find frequency bands with most power
if iShift == nShifts,
vPowShiftRef1 = mean(mPowRef1( (iShift-1)*blockLen+1:end,:));
vPowShiftTest1 = mean(mPowTest1((iShift-1)*blockLen+1:end,:));
vPowShiftRef2 = mean(mPowRef2( (iShift-1)*blockLen+1:end,:));
vPowShiftTest2 = mean(mPowTest2((iShift-1)*blockLen+1:end,:));
else
vPowShiftRef1 = mean(mPowRef1( (iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,:));
vPowShiftTest1 = mean(mPowTest1((iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,:));
vPowShiftRef2 = mean(mPowRef2( (iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,:));
vPowShiftTest2 = mean(mPowTest2((iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,:));
end
vPowShiftRef1 = real(10.*log10(vPowShiftRef1));
vPowShiftTest1 = real(10.*log10(vPowShiftTest1));
vPowShiftRef2 = real(10.*log10(vPowShiftRef2));
vPowShiftTest2 = real(10.*log10(vPowShiftTest2));
maxPowRef = max(max(vPowShiftRef1,vPowShiftRef2));
maxPowTest = max(max(vPowShiftTest1,vPowShiftTest2));
idxRef1 = find(vPowShiftRef1 > maxPowRef-10);
idxRef2 = find(vPowShiftRef2 > maxPowRef-10);
idxTest1 = find(vPowShiftTest1 > maxPowTest-10);
idxTest2 = find(vPowShiftTest2 > maxPowTest-10);
idxUnion = union([idxRef1 idxRef2],[idxTest1 idxTest2]);
for jBand = idxUnion,
% load IVS
idxBandDiff = abs(idxUnion-jBand);
idxNeighborBands = find(idxBandDiff == 1 | idxBandDiff == 0); %
if iShift == nShifts,
vIVSbandRef = mean(mIVSref( (iShift-1)*blockLen+1:end,idxUnion(idxNeighborBands)),2);
vIVSbandTest = mean(mIVStest((iShift-1)*blockLen+1:end,idxUnion(idxNeighborBands)),2);
else
vIVSbandRef = mean(mIVSref( (iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,idxUnion(idxNeighborBands)),2);
vIVSbandTest = mean(mIVStest((iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,idxUnion(idxNeighborBands)),2);
end
% ILD sub-measure, equation 6
% -------------------------------------------------------------------------
if iShift == nShifts,
mILDbandRef = mILDref( (iShift-1)*blockLen+1:end,jBand);
mILDbandTest = mILDtest((iShift-1)*blockLen+1:end,jBand);
else
mILDbandRef = mILDref( (iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,jBand);
mILDbandTest = mILDtest((iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,jBand);
end
mILDbandRef(mILDbandRef > 10) = 10;
mILDbandRef(mILDbandRef < -10) = -10;
mILDbandTest(mILDbandTest > 10) = 10;
mILDbandTest(mILDbandTest < -10) = -10;
ILDsumPlusRef = sum(abs(mILDbandRef(mILDbandRef > 0)));
ILDsumMinusRef = sum(abs(mILDbandRef(mILDbandRef < 0)));
ILDsumPlusTest = sum(abs(mILDbandTest(mILDbandTest > 0)));
ILDsumMinusTest = sum(abs(mILDbandTest(mILDbandTest < 0)));
ILDdiff = abs(ILDsumPlusRef - ILDsumPlusTest) + abs(ILDsumMinusRef - ILDsumMinusTest);
% ITD sub-measure, equation 7-9
% -------------------------------------------------------------------------
if iShift == nShifts,
vITDbandRef = mITDref( (iShift-1)*blockLen+1:end,jBand);
vITDbandTest = mITDtest((iShift-1)*blockLen+1:end,jBand);
else
vITDbandRef = mITDref( (iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,jBand);
vITDbandTest = mITDtest((iShift-1)*blockLen+1:(iShift-1)*blockLen+blockLen,jBand);
end
vITDbandRef(abs(vITDbandRef) > 700e-6) = NaN;
idxHighIVSref = vIVSbandRef > ivsThres & [diff(vIVSbandRef)>=0 ; zeros(1,1)];
vITDbandRef(~idxHighIVSref) = NaN;
vITDbandTest(abs(vITDbandTest) > 700e-6) = NaN;
idxHighIVStest = vIVSbandTest > ivsThres & [diff(vIVSbandTest)>=0 ; zeros(1,1)];
vITDbandTest(~idxHighIVStest) = NaN;
ITDsumPlusRef = nansum(abs(vITDbandRef(vITDbandRef >= 0)))./(blockLen-sum(isnan(vITDbandRef))+eps);
ITDsumMinusRef = nansum(abs(vITDbandRef(vITDbandRef < 0)))./(blockLen-sum(isnan(vITDbandRef))+eps);
ITDsumPlusTest = nansum(abs(vITDbandTest(vITDbandTest >= 0)))./(blockLen-sum(isnan(vITDbandTest))+eps);
ITDsumMinusTest = nansum(abs(vITDbandTest(vITDbandTest < 0)))./(blockLen-sum(isnan(vITDbandTest))+eps);
ITDdiff = abs(ITDsumPlusRef - ITDsumPlusTest) + abs(ITDsumMinusRef - ITDsumMinusTest);
% IVS sub-measure, equation 12-15
% -------------------------------------------------------------------------
K = k.r./(1+exp(k.s.*(fc(jBand)-k.t))) + k.a;
N = n.r./(1+exp(n.s.*(fc(jBand)-n.t))) + n.a;
dIVSref = exp(K+N) - exp(K.*vIVSbandRef+N);
dIVStest = exp(K+N) - exp(K.*vIVSbandTest+N);
IVSdiff = abs(mean(dIVSref)-mean(dIVStest));
% -------------------------------------------------------------------------
if iShift == nShifts && nResiSamples > 0,
ILDdiff = ILDdiff.*(blockLen/nResiSamples);
ITDdiff = ITDdiff.*(blockLen/nResiSamples);
end
vILDdiffs = [vILDdiffs ILDdiff];
vITDdiffs = [vITDdiffs ITDdiff];
vIVSdiffs = [vIVSdiffs IVSdiff];
end
end
ILDdiff = mean(vILDdiffs);
ITDdiff = mean(vITDdiffs);
IVSdiff = mean(vIVSdiffs);
load ARESmodel
binQ = arespredict(stBinModel, [ILDdiff ITDdiff IVSdiff]);
BAMQ/functions/BAMQFrontEnd.m 0 → 100644
View file @ bba7b439
function [mITD, mIVS, mILD, mPow1, mPow2, fc, fsNew] = BAMQFrontEnd(insig,fs)
%DIETZ2011 Dietz 2011 binaural model changed for BAMQ
% Usage: [...] = BinQFrontEnd(insig,fs);
%
% Input parameters:
% insig : binaural signal for which values should be calculated
% fs : sampling rate (Hz)
%
% Reference:
% M. Dietz, S. D. Ewert, and V. Hohmann. Auditory model based direction
% estimation of concurrent speakers from binaural signals. Speech
% Communication, 53(5):592 - 605, 2011.
%% Initialization
fsNew = 6000;
% Preprocessing parameters
vMiddleEarThr = [500 2000]; % Bandpass freqencies for middle ear transfer
middleEarOrder = 1;
compressionPower = 0.4;
haircellLowOrder = 5;
haircellLowCutoffHz = 770;
% Gammatone filterbank parameters
fLow = 50;
fHigh = 14000;
fBase = 700;
filtersPerERB = 1;
gammaOrder = 4;
bandwidthFactor = 1;
% modulation filter parameters
fModLow1 = 500;
fModLow2 = 600;
orderModLow1 = 2;
orderModLow2 = 10;
modGammaCenterFreqHz = 300;
modGammaBandwidthHz = 300;
modGammaAttenuationDB = 3;
modGammaOrder = 1;
% temporal resolution
tauCycles = 5; % in cycles
%% middle ear band pass filtering
[z,p,k] = butter(middleEarOrder,vMiddleEarThr/(fs/2));
[sos,g] = zp2sos(z,p,k);
h=dfilt.df2sos(sos,g);
mSignalME = filter(h,insig);
%% Filterbank analysis
% create filterbank
analyzer = gfb_analyzer_new(fs, fLow, fBase, fHigh,...
filtersPerERB,gammaOrder,bandwidthFactor);
fc = analyzer.center_frequencies_hz;
analyzer2 = analyzer;
idxHighestFine = find(fc < 1400,1,'last');
% apply filterbank
[mSignalIE, analyzer] = gfb_analyzer_process(analyzer, mSignalME(:,1));
[mSignalIE2, analyzer2] = gfb_analyzer_process(analyzer2, mSignalME(:,2));
mSignalIE = mSignalIE.';
mSignalIE2 = mSignalIE2.';
clear vSignalME
%% haircell processing
% half-wave rectification
mSignalRect = max(real(mSignalIE),0);
mSignalRect2 = max(real(mSignalIE2),0);
clear mSignalIE mSignalIE2
% compression
mSignalHC = mSignalRect.^compressionPower;
mSignalHC2 = mSignalRect2.^compressionPower;
clear mSignalRect mSignalRect2
% audio signal power calculation
mSignalHC1re = resample(mSignalHC,fsNew,fs);
mSignalHC2re = resample(mSignalHC2,fsNew,fs);
mPow1 = abs(mSignalHC1re.^(0.4^-1)).^2;
mPow2 = abs(mSignalHC2re.^(0.4^-1)).^2;
% lowpass
[z,p,k] = butter(haircellLowOrder,haircellLowCutoffHz/(fs/2));
[sos,g] = zp2sos(z,p,k);
h=dfilt.df2sos(sos,g);
mSignalHC = filter(h,mSignalHC);
mSignalHC2 = filter(h,mSignalHC2);
%% finestructure filter
% processing the haircell output with a fine structure filter
[mSignalFine, mSignalFine2] =...
gfb_envelope_filter(mSignalHC(:,1:idxHighestFine), mSignalHC2(:,1:idxHighestFine), fs, fc(1:idxHighestFine),...
bandwidthFactor,0,gammaOrder);
% modulation filter part 1, derivation to remove DC
mSignalModNoDC = filter([1 -1],1,mSignalHC(:,idxHighestFine+1:end));
mSignalModNoDC2 = filter([1 -1],1,mSignalHC2(:,idxHighestFine+1:end));
clear mSignalHC mSignalHC2
% calculation of interaural functions from haircell fine structure
stInterauralParamFine = dietzInterauralFunctions4binQ(...
mSignalFine, mSignalFine2, fc(1:idxHighestFine), tauCycles, fs, fsNew);
%% modulation filter part2
% mod filter lowpass 1 and 2
[z,p,k] = butter(orderModLow1,fModLow1/(fs/2));
[sos,g] = zp2sos(z,p,k);
h=dfilt.df2sos(sos,g);
mSignalModLow1 = filter(h,mSignalModNoDC);
mSignalModLow12 = filter(h,mSignalModNoDC2);
clear mSignalModNoDC mSignalModNoDC2
[z,p,k] = butter(orderModLow2,fModLow2/(fs/2));
[sos,g] = zp2sos(z,p,k);
h=dfilt.df2sos(sos,g);
mSignalModLow2 = filter(h,mSignalModLow1);
mSignalModLow22 = filter(h,mSignalModLow12);
clear mSignalModLow1 mSignalModLow12
[mSignalMod, mSignalMod2] =...
gfb_envelope_filter(mSignalModLow2, mSignalModLow22, fs,...
modGammaCenterFreqHz, modGammaBandwidthHz, ...
modGammaAttenuationDB, modGammaOrder);
clear mSignalModLow2 mSignalModLow22
% calculation of interaural functions from haircell modulation
stInterauralParamMod = dietzInterauralFunctions4binQ(...
mSignalMod, mSignalMod2, modGammaCenterFreqHz + zeros(1,length(mSignalMod(1,:))),...
tauCycles, fs, fsNew);
clear mSignalMod mSignalMod2
%% combine fine and mod
mIPD = [stInterauralParamFine.ipd_lp stInterauralParamMod.ipd_lp];
mIVS = [stInterauralParamFine.ic stInterauralParamMod.ic];
mFinst = [stInterauralParamFine.f_inst stInterauralParamMod.f_inst];
%% IPD2ITD
mITD = 1/(2*pi)*mIPD./mFinst;
%% ILD
[z,p,k] = butter(2,30/(fsNew/2));
[sos,g] = zp2sos(z,p,k);
h=dfilt.df2sos(sos,g);
mSignalLP = filter(h,mSignalHC1re);
mSignalLP2 = filter(h,mSignalHC2re);
mILD = 20*log10(max(mSignalLP2,1e-4)./max(mSignalLP,1e-4))/0.4;
BAMQ/functions/arespredict.m 0 → 100644
View file @ bba7b439
function Yq = arespredict(model, Xq)
% arespredict
% Predicts output values for the given query points Xq using an ARES model.
%
% Call:
% Yq = arespredict(model, Xq)
%
% Input:
% model : ARES model
% Xq : Inputs of query data points (Xq(i,:)), i = 1,...,nq
%
% Output:
% Yq : Predicted outputs of the query data points (Yq(i)),
% i = 1,...,nq
% =========================================================================
% ARESLab: Adaptive Regression Splines toolbox for Matlab/Octave
% Author: Gints Jekabsons (gints.jekabsons@rtu.lv)
% URL: http://www.cs.rtu.lv/jekabsons/
%
% Copyright (C) 2009-2011 Gints Jekabsons
%
% 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/>.
% =========================================================================
% Last update: November 9, 2009
if nargin < 2
error('Too few input arguments.');
end
X = ones(size(Xq,1),length(model.knotdims)+1);
if model.trainParams.cubic
for i = 1 : length(model.knotdims)
X(:,i+1) = createbasisfunction(Xq, X, model.knotdims{i}, model.knotsites{i}, ...
model.knotdirs{i}, model.parents(i), model.minX, model.maxX, model.t1(i,:), model.t2(i,:));
end
else
for i = 1 : length(model.knotdims)
X(:,i+1) = createbasisfunction(Xq, X, model.knotdims{i}, model.knotsites{i}, ...
model.knotdirs{i}, model.parents(i), model.minX, model.maxX);
end
end
Yq = X * model.coefs;
return
BAMQ/functions/audspacebw.m 0 → 100644
View file @ bba7b439
function [y,n] = audspacebw(flow,fhigh,varargin)
%AUDSPACEBW Auditory scale points specified by bandwidth
% Usage: y=audspacebw(flow,fhigh,bw,hitme);
% y=audspacebw(flow,fhigh,bw);
% y=audspacebw(flow,fhigh);
% [y,n]=audspacebw(...);
%
% AUDSPACEBW(flow,fhigh,bw,scale) computes a vector containing values
% equistantly scaled between frequencies flow and fhigh on the
% selected auditory scale. All frequencies are specified in Hz.The
% distance between two consecutive values is bw on the selected scale,
% and the points will be centered on the scale between flow and fhigh.
%
% See the help on FREQTOAUD to get a list of the supported values of the
% scale parameter.
%
% AUDSPACEBW(flow,fhigh,bw,hitme,scale) will do as above, but one of
% the points is quaranteed to be the frequency hitme.
%
% [y,n]=AUDSPACEBW(...) additionally returns the number of points n in
% the output vector y.
%
% See also: freqtoaud, audspace, audfiltbw
%
% Url: http://ltfat.sourceforge.net/doc/auditory/audspacebw.php
% Copyright (C) 2005-2012 Peter L. Soendergaard.
% This file is part of LTFAT version 1.1.1
% 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/>.
% AUTHOR : Peter L. Soendergaard
% changed: Feb 2013 (jhf)
% ------ Checking of input parameters ---------
if nargin<2
error('%s: Too few input parameters.',upper(mfilename));
end;
if ~isnumeric(flow) || ~isscalar(flow) || flow<0
error('%s: flow must be a non-negative scalar.',upper(mfilename));
end;
if ~isnumeric(fhigh) || ~isscalar(fhigh) || fhigh<0
error('%s: fhigh must be a non-negative scalar.',upper(mfilename));
end;
if flow>fhigh
error('%s: flow must be less than or equal to fhigh.',upper(mfilename));
end;
% import = {'freqtoaud'};
hitme = varargin{2};
bw = varargin{1};
if ~isnumeric(bw) || ~isscalar(bw) || bw<=0
error('%s: bw must be a positive scalar.',upper(mfilename));
end;
%% ------ Computation --------------------------
if isempty(hitme)
% Convert the frequency limits to auds.
% audlimits = freqtoaud([flow,fhigh],flags.audscale);
audlimits = 9.2645*sign([flow,fhigh]).*log(1+abs([flow,fhigh])*0.00437);
audrange = audlimits(2)-audlimits(1);
% Calculate number of points, excluding final point
n = floor(audrange/bw);
% The remainder is calculated in order to center the points
% correctly between flow and fhigh.
remainder = audrange-n*bw;
audpoints = audlimits(1)+(0:n)*bw+remainder/2;
% Add the final point
n=n+1;
else
% Convert the frequency limits to auds.
% audlimits = freqtoaud([flow,fhigh,hitme],flags.audscale);
audlimits = 9.2645*sign([flow,fhigh,hitme]).*log(1+abs([flow,fhigh,hitme])*0.00437);
audrangelow = audlimits(3)-audlimits(1);
audrangehigh = audlimits(2)-audlimits(3);
% Calculate number of points, exluding final point.
nlow = floor(audrangelow/bw);
nhigh = floor(audrangehigh/bw);
audpoints=(-nlow:nhigh)*bw+audlimits(3);
n=nlow+nhigh+1;
end;
% y = audtofreq(audpoints,flags.audscale);
y = (1/0.00437)*sign(audpoints).*(exp(abs(audpoints)/9.2645)-1);
BAMQ/functions/createbasisfunction.m 0 → 100644
View file @ bba7b439
function Xn = createbasisfunction(X, Xtmp, knotdims, knotsites, knotdirs, parent, minX, maxX, t1, t2)
% Creates a list of response values of a defined basis function for either
% a piecewise-linear or piecewise-cubic model.
% =========================================================================
% ARESLab: Adaptive Regression Splines toolbox for Matlab/Octave
% Author: Gints Jekabsons (gints.jekabsons@rtu.lv)
% URL: http://www.cs.rtu.lv/jekabsons/
%
% Copyright (C) 2009-2011 Gints Jekabsons
%
% 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/>.
% =========================================================================
% Last update: November 9, 2009
n = size(X,1);
% if the basis function will always be zero, return NaN
if ((knotdirs(end) > 0) && (knotsites(end) >= maxX(knotdims(end)))) || ...
((knotdirs(end) < 0) && (knotsites(end) <= minX(knotdims(end))))
Xn = NaN(n,1);
return
end
% if the basis function has a direct parent in the model, use it to speed-up the calculations
if parent > 0
Xn = Xtmp(:,parent+1);
start = length(knotdims);
else
Xn = ones(n,1);
start = 1;
end
if nargin < 10 % piecewise-linear
for i = start : length(knotdims)
if knotdirs(i) > 0
z = X(:,knotdims(i)) - knotsites(i);
else
z = knotsites(i) - X(:,knotdims(i));
end
Xn = Xn .* max(0,z);