Skip to content

GitLab

C
CombinedAudioQualityModel
Commit 5efefaec authored by dualberger's avatar dualberger
Browse files
Options

reorganize file structure part 1

parent bba7b439
Hide whitespace changes
Inline Side-by-side
Showing with 930 additions and 216 deletions
Example_combAudioQual.asv deleted 100644 → 0
View file @ bba7b439
% EXAMPLE for the usage of the combined audio quality model that is described
% in Flener et al. (2019). It is a reference-based audio quality model and based on
% monaural (intensity and amplitude modulation) and binaural (interaural level and
% phase differences, interaural vector strength) cues to predict subjective audio
% quality ratings.
cur_dir=pwd;
addpath(genpath(cur_dir)); % include all folders from the current path
%% input signals
if verLessThan('matlab','8.0')
% reference signal (clean)
[RefSig, fsRef] = wavread('Stimuli/REF2_short.wav');
% test signal (processed)
[TestSig, fsTest] = wavread('Stimuli/RAI_P2_2_short.wav');
else
% reference signal (clean)
[RefSig, fsRef] = audioread('Stimuli/REF2_short.wav');
% test signal (processed)
[TestSig, fsTest] = audioread('Stimuli/RAI_P2_2_short.wav');
end
% compare sampling frequencies
if fsTest ~= fsRef,
error('signals have different sampling frequencies')
else
fs = fsTest;
end
minlen = min([length(TestSig), length(RefSig)]);
TestSig = TestSig(1:minlen,1);
RefSig = RefSig(1:minlen,1);
% calculate objective perceptual measures
stOut = GPSMqBin(RefSig, TestSig, fs); % monaural model output
binQ(aa,1) = BAMQ(RefSig, TestSig, fs); % binaural model output
% **********************************************
% **********Combine outputs of BAMq and GPSMq***
% **********************************************
obj_meas_mon=stOut.OPM(:,1); % monaural
obj_meas_bin=binQ(:,1); %
%*********f3************
mon_tmp=0.0528.*obj_meas_mon;
bin_tmp=0.0078*obj_meas_bin;
obj_meas=min(log10(mon_tmp),bin_tmp);
disp(stOut)
% OPM: 5.7203
% OPM_raw: 1.7207
% SNR_dc: 0.3909
% SNR_ac: 1.0952
\ No newline at end of file
Example_combAudioQual.m
View file @ 5efefaec
......@@ -3,7 +3,7 @@
% monaural (intensity and amplitude modulation) and binaural (interaural level and
% phase differences, interaural vector strength) cues to predict subjective audio
% quality ratings.
tic
cur_dir=pwd;
addpath(genpath(cur_dir)); % include all folders from the current path
%% input signals
......@@ -45,21 +45,25 @@ mon_tmp=0.0528.*obj_meas_mon;
bin_tmp=0.0078*obj_meas_bin;
obj_meas=min(log10(mon_tmp),bin_tmp);
toc
disp('***********************')
disp('****monaural measures**')
disp('***********************')
disp(stOut)
% out: 5.7203
% outFix: 1.7207
% SNR_dc: 0.3909
% SNR_dc_fix: 1.0952
% SNR_ac: 0.3909
% SNR_ac_fix: 1.0952
% SNR_dc: 0.0179
% SNR_ac: 0.0225
% SNR_dc_fix: 0.0179
% SNR_ac_fix: 0.0225
% out: 132.6339
% outFix: 132.7494
disp('***********************')
disp('***binaural measures***')
disp('***********************')
disp(binQ)
% binQ: 100
disp('***********************')
disp('****overall measures***')
disp('***********************')
disp(obj_meas)
\ No newline at end of file
disp(obj_meas)
% overall measure: 0.7800
\ No newline at end of file
GPSMq/GPSMqBin.m
View file @ 5efefaec
function stOut = GPSMqBin(RefSig, TestSig, fs)
% input: RefSig ... reference signal, 1-dim
% TestSig ... test signal, 1-dim, time-aligned with RefSig
% GPSMqBin Main function for calculation of audio quality differences between an
% unprocessed reference signal and a processed test signal. At first model
% parameters are defined, followed by the front end and back end
% calculations resulting in the objective perceptual measure. For
% more detailed information see Biberger et al. (2018) and
% Flener et al. (2019)
%
% input: RefSig ... reference signal, 2-dim
% TestSig ... test signal, 2-dim, time-aligned with RefSig
% fs ... sampling frequency
%
% output:
% output: stOut ... struct that contains the objective perceptual
% measure (OPM) as it is described in Flener et al.(2019)
% in Eq. 15 and 16, intensity SNRs (SNR_dc), and envelope
% power SNRs (SNR_ac)
%
%
% Usage: stOut = GPSMqBin(RefSig, TestSig, fs)
% authors: jan-hendrik.flessner@uni-oldenburg.de
% thomas.biberger@uni-oldenburg.de
% updated 2019-06-02
% Back-end related paramter settings
params.auditory_filt_range.start=8; % audio freqs=[63 80 100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500 16000]
......@@ -54,12 +69,12 @@ stOut.SNR_ac = 0.5 * (SNR_ac_l + SNR_ac_r);
stOut.SNR_dc_fix = 0.5 * (SNR_dc_fix_l + SNR_dc_fix_r);
stOut.SNR_ac_fix = 0.5 * (SNR_ac_fix_l + SNR_ac_fix_r);
%% transformation of a single SNR-value into a perceptual measure
%% transformation of a single SNR-value into a perceptual measure (OPM*-> see Eq. 15 in Flener et al.,2019)
f = 10*log10(stOut.SNR_dc + stOut.SNR_ac + eps);
b = [-4.08695250298565;75.6467755339438];
stOut.out = [f 1]*b;
%% transformation of a single SNR-value into a perceptual measure FIX
%% transformation of a single SNR-value into a perceptual measure FIX (OPM*dual-> see Eq. 16 in Flener et al.,2019)
f = 10*log10(stOut.SNR_dc_fix + stOut.SNR_ac_fix + eps);
b = [-4.15742483856597;74.7791007678067];
stOut.outFix = [f 1]*b;
......
GPSMq/accessoires/mfb2_GPSM.m 0 → 100644
View file @ 5efefaec
% mfb2_GPSM.m - applies modulation filterbank, mfb2_GPSM.m is based on
% mfb2.m from Stephan Ewert and Torsten Dau
%
% Usage: [inf, out] = mfb2_GPSM(in,style,fs)
%
% in = input column vector
% 1 = 3-dB overlap, 0.5 = every second filter, 2 = two times as much filters ... .
% style = toggles the overall 150 Hz lowpass on (1) / off (2).
% fs = sampling rate in Hz,
% should be greater than 8000 Hz to avoid errors due to the bilinear transform.
%
% inf = center frequencies of the modulation filters.
% [out1,out2, ...,outn] = each column of matrix 'out' contains the output of
% a single modulation filter.
%
% copyright (c) 1999 Stephan Ewert and Torsten Dau, Universitaet Oldenburg
% modified by thomas.biberger@uni-oldenburg.de, 2018-10-11
function [inf, out] = mfb2_GPSM(in,style,fs)
% if fs < 8000
% warning('mfb2: sample rate lower than 8000 Hz')
% end
Q = 1; % default
lpcut = 150; % default 150 Hz
sw = 1;
mf=[2 4 8 16 32 64 128 256]; % add from thomas 03.04
b2lpcut=1; % lowpass cutoff freq
wlp = 2*pi*lpcut/fs;
wb2lp = 2*pi*b2lpcut/fs; % orig. from Stephan & Torsten
[b1,a1] = folp(wlp);
[b2,a2] = solp(wb2lp,1/sqrt(2)); % orig. from Stephan & Torsten
switch style % switch overall lowpass
case 1
outtmp = filter(b1,a1,in);
case 2
outtmp = in;
end
switch sw
case 0 % only one modulation filter
w0 = 2*pi*mf/fs;
[b3,a3] = sobp(w0,Q);
out = filter(b3,a3,outtmp);
inf = mf;
case 1 % lowpass and modulation filter(s)
out = zeros(length(in),length(mf)+1);
out(:,1) = filter(b2,a2,outtmp); % lowpass for dc-component
for i=1:length(mf)
w0 = 2*pi*mf(i)/fs;
[b3,a3] = sobp(w0,Q);
out(:,i+1) = filter(b3,a3,outtmp);
end
inf = [0 mf];
case 2 % only modulation filters
out=zeros(length(in),length(mf));
for i=1:length(mf)
w0 = 2*pi*mf(i)/fs;
[b3,a3] = sobp(w0,Q);
out(:,i) = filter(b3,a3,outtmp);
end
inf = mf;
case 3 % only lowpass
out = filter(b2,a2,outtmp);
inf = 0;
end
% subfunctions
% second order bandpass filter
function [b,a] = sobp(w0,Q)
W0 = tan(w0/2);
B0 = W0/Q;
b = [B0; 0; -B0];
a = [1 + B0 + W0^2; 2*W0^2 - 2; 1 - B0 + W0^2];
b = b/a(1);
a = a/a(1);
% first order lowpass filter
function [b,a] = folp(w0);
W0 = tan(w0/2);
b = [W0, W0]/(1 + W0);
a = [1,(W0 - 1)/(1 + W0)];
% second order Butterworth lowpass filter
function [b,a] = solp(w0,Q)
W0 = tan(w0/2);
b = [1; 2; 1];
a = [1 + 1/(Q*W0) + 1/W0^2; -2/W0^2 + 2; 1 - 1/(Q*W0) + 1/W0^2];
b = b/a(1);
a = a/a(1);
\ No newline at end of file
GPSMq/accessoires/moving_average.m 0 → 100644
View file @ 5efefaec
function [X,Zf] = moving_average(N,X,Zi,dim)
% Like filter() for the special case of moving-average kernels.
% [Y,Zf] = moving_average(N,X,Zi,Dim)
%
% This is an overall very fast implementation whose running time does now grow with N (beyond
% N=500). The algorithm does not run into numerical problems for large data sizes unlike the usual
% cumsum-based implementations.
%
% In:
% N : filter length in samples
%
% X : data matrix
%
% Zi : initial filter conditions (default: [])
%
% Dim : dimension along which to filter (default: first non-singleton dimension)
%
% Out:
% X : the filtered data
%
% Zf : final filter conditions
%
% See also:
% filter
%
% Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
% 2012-01-10
% determine the dimension along which to filter
if nargin <= 3
if isscalar(X)
dim = 1;
else
dim = find(size(X)~=1,1);
end
end
% empty initial state
if nargin <= 2
Zi = []; end
lenx = size(X,dim);
if lenx == 0
% empty X
Zf = Zi;
else
if N < 500
% short N: use filter
[X,Zf] = filter(ones(N,1)/N,1,X,Zi,dim);
else
% we try to avoid permuting dimensions below as this would increase the running time by ~3x
if ndims(X) == 2
if dim == 1
% --- process along 1st dimension ---
if isempty(Zi)
% zero initial state
Zi = zeros(N,size(X,2));
elseif size(Zi,1) == N-1
% reverse engineer filter's initial state (assuming a moving average)
tmp = diff(Zi(end:-1:1,:),1,1);
Zi = [tmp(end:-1:1,:); Zi(end,:)]*N;
Zi = [-sum(Zi,1); Zi];
elseif ~isequal(size(Zi),[N,size(X,2)])
error('These initial conditions do not have the correct format.');
end
% pre-pend initial state & get dimensions
Y = [Zi; X]; M = size(Y,1);
% get alternating index vector (for additions & subtractions)
I = [1:M-N; 1+N:M];
% get sign vector (also alternating, and includes the scaling)
S = [-ones(1,M-N); ones(1,M-N)]/N;
% run moving average
X = cumsum(bsxfun(@times,Y(I(:),:),S(:)),1);
% read out result
X = X(2:2:end,:);
% construct final state
if nargout > 1
Zf = [-(X(end,:)*N-Y(end-N+1,:)); Y(end-N+2:end,:)]; end
else
% --- process along 2nd dimension ---
if isempty(Zi)
% zero initial state
Zi = zeros(N,size(X,1));
elseif size(Zi,1) == N-1
% reverse engineer filter's initial state (assuming a moving average)
tmp = diff(Zi(end:-1:1,:),1,1);
Zi = [tmp(end:-1:1,:); Zi(end,:)]*N;
Zi = [-sum(Zi,1); Zi];
elseif ~isequal(size(Zi),[N,size(X,1)])
error('These initial conditions do not have the correct format.');
end
% pre-pend initial state & get dimensions
Y = [Zi' X]; M = size(Y,2);
% get alternating index vector (for additions & subtractions)
I = [1:M-N; 1+N:M];
% get sign vector (also alternating, and includes the scaling)
S = [-ones(1,M-N); ones(1,M-N)]/N;
% run moving average
X = cumsum(bsxfun(@times,Y(:,I(:)),S(:)'),2);
% read out result
X = X(:,2:2:end);
% construct final state
if nargout > 1
Zf = [-(X(:,end)*N-Y(:,end-N+1)) Y(:,end-N+2:end)]'; end
end
else
% --- ND array ---
[X,nshifts] = shiftdim(X,dim-1);
shape = size(X); X = reshape(X,size(X,1),[]);
if isempty(Zi)
% zero initial state
Zi = zeros(N,size(X,2));
elseif size(Zi,1) == N-1
% reverse engineer filter's initial state (assuming a moving average)
tmp = diff(Zi(end:-1:1,:),1,1);
Zi = [tmp(end:-1:1,:); Zi(end,:)]*N;
Zi = [-sum(Zi,1); Zi];
elseif ~isequal(size(Zi),[N,size(X,2)])
error('These initial conditions do not have the correct format.');
end
% pre-pend initial state & get dimensions
Y = [Zi; X]; M = size(Y,1);
% get alternating index vector (for additions & subtractions)
I = [1:M-N; 1+N:M];
% get sign vector (also alternating, and includes the scaling)
S = [-ones(1,M-N); ones(1,M-N)]/N;
% run moving average
X = cumsum(bsxfun(@times,Y(I(:),:),S(:)),1);
% read out result
X = X(2:2:end,:);
% construct final state
if nargout > 1
Zf = [-(X(end,:)*N-Y(end-N+1,:)); Y(end-N+2:end,:)]; end
X = reshape(X,shape);
X = shiftdim(X,ndims(X)-nshifts);
end
end
end
GPSMq/mreGPSMq/BE/LpFB_GPSM_v2.m 0 → 100644
View file @ 5efefaec
function out = LpFB_GPSM_v2(in,freq,fs)
% LpFB_GPSM.m Lowpass-Filterbank applied to power features per auditory
% channel. Temporal-integrated power features are used to normalize short-time envelope power feature.
% (make them independent from the short-time intensity of the carrier)
% INPUT
% in: envelope of the corresponding auditory channel
% freq: vector that gives the modulation frequencies [1,2,4,...,256 Hz]
% fs: sampling frequency
% OUTPUT
% out: matrix containing dc-values;
% different window length (depending on the modulation filter center frequency)
% are used for computing dc-values,
%
% Usage: out = LpFB_GPSM_v2(in,freq,fs)
% author: thomas.biberger@uni-oldenburg
% date: 2013-07-22
% modified: 2018-10-12
out=zeros(length(in),length(freq)); % preallocation
factor=1.5; % 1.4
count =0;
if freq(1)==0;
start=2;
elseif freq>0;
start=1;
end
for ii=start:length(freq);
if freq(ii)==0
else
if freq(ii)<=8;
N=round(fs/(freq(ii)/factor));
% average window
b=ones(1,N)/N; % numerator of the rectangular moving average window
% out(:,ii) = dlfconvComp(b,in);
out(:,ii)=moving_average(N,in);
elseif freq(ii)>8 && count==0;
factor=round(freq(ii)/8);
count=count+1;
N=round(fs/(freq(ii)/factor));
% average window
b=ones(1,N)/N; % numerator of the rectangular moving average window
% out(:,ii) = dlfconvComp(b,in);
out(:,ii)=moving_average(N,in);
else freq(ii)>8 && count>=1;
out(:,ii)=out(:,ii-1);
end
end
end
GPSMq/mreGPSMq/FE/mreGPSMq_init.m
View file @ 5efefaec
% revision 1.00.3 beta, 13-01-2005 14:24
function [def simdef simwork]= mreGPSMq_init(def,simdef)
% sets some simwork variables depending on simdef structure
% mreGPSMq_init.m some pre-calculations of e.g., filter coefficients that have
% been calculated only one time
%
% global def
% global work
% global simdef
% global simwork
%% tdepsm specific
% INPUT:
% def: struct for various parameters, e.g.,fs
% simdef: struct for model related parameters, e.g., filter order,
% cut-off frequency
%
% OUTPUT:
% def: struct for stimulus related parameters,e.g.,fs
% simdef: struct for model related parameters, e.g., filter order,
% cut-off frequency
% simwork: struct for parameters applied during simulation, e.g.,
% filter coefficients
%
% Usage: [def simdef simwork]= mreGPSMq_init(def,simdef)
% thomas.biberger@uni-oldenburg.de;
% date: 2018-04-12
% Resampling
......@@ -28,7 +35,6 @@ simwork.analyzer = gfb_analyzer_new_tb(def.samplerate, simdef.gt_MinCF, simdef.g
simdef.gt_Dens);
% iso hearing threshold
% vsDesF=0:def.samplerate/def.intervallen:floor((def.intervallen/2-1))*def.samplerate/def.intervallen; % old: freq vector needed for isothr
vsDesF=0:def.samplerate/def.intervallen:floor((def.intervallen/2))*def.samplerate/def.intervallen; % NEW: freq vector needed for isothr
[simwork.vIsoThrDB, simwork.vsF] = iso389_7(vsDesF);
......
GPSMq/mreGPSMq/FE/mreGPSMq_init_v2.m 0 → 100644
View file @ 5efefaec
function [def simdef simwork]= mreGPSMq_init_v2(def,simdef)
% mreGPSMq_init.m some pre-calculations of e.g., filter coefficients that have
% been calculated only one time
%
% INPUT:
% def: struct for various parameters, e.g.,fs
% simdef: struct for model related parameters, e.g., filter order,
% cut-off frequency
%
% OUTPUT:
% def: struct for stimulus related parameters,e.g.,fs
% simdef: struct for model related parameters, e.g., filter order,
% cut-off frequency
% simwork: struct for parameters applied during simulation, e.g.,
% filter coefficients
%
% Usage: [def simdef simwork]= mreGPSMq_init(def,simdef)
% thomas.biberger@uni-oldenburg.de;
% date: 2018-11-27
% Resampling
if strcmp(simdef.resampling,'on')==1;
def.downsample_factor=2;
def.samplerate_down=def.samplerate/def.downsample_factor; % downsample signal by factor 2
elseif strcmp(simdef.resampling,'off')==1;
def.samplerate_down=def.samplerate; % no downsampling
else
disp('Please define if resampling should be used!')
end
% parameter (CFs etc.) of the GTFB
simwork.analyzer = gfb_analyzer_new_tb(def.samplerate, simdef.gt_MinCF, simdef.gt_align, simdef.gt_MaxCF,...
simdef.gt_Dens);
% iso hearing threshold
vsDesF=0:def.samplerate/def.intervallen:floor((def.intervallen/2))*def.samplerate/def.intervallen; % NEW: freq vector needed for isothr
[simwork.vIsoThrDB, simwork.vsF] = iso389_7(vsDesF);
% 1st order butterworth-filter for 150Hz-envelope-lp
simdef.env_lp_order=1;
simdef.env_lp_cutoff=150;
simdef.env_lp_cutoff_norm=simdef.env_lp_cutoff/(def.samplerate/2);
[simwork.b_env_lp,simwork.a_env_lp] = butter(simdef.env_lp_order,simdef.env_lp_cutoff_norm,'low');
% Lowpass-Filterbank applied to power features per auditory
% channel. Temporal-integrated power features are used to normalize short-time envelope power feature.
% (make them independent from the short-time intensity of the carrier)
mod_freq=[0 1 2 4 8 16 32 64 128 256];
factor=1.5;
simdef.LpFB.b=zeros(round(def.samplerate_down/(mod_freq(2)/factor)),length(mod_freq));
simdef.LpFB.N_per_freq=zeros(length(mod_freq),1);
for ii=1:length(mod_freq)
if mod_freq(ii)==0
else
if mod_freq(ii)>8
factor=round(mod_freq(ii)/8);
end
N=round(def.samplerate_down/(mod_freq(ii)/factor));
simdef.LpFB.N_per_freq(ii,1)=N;
% b(:,ii)=ones(1,N)/N; % numerator of the rectangular moving average window
simdef.LpFB.b(1:N,ii)=ones(N,1)/N; % numerator of the rectangular moving average window
% b=b';
end
end
end
GPSMq/mreGPSMq/FE/mreGPSMq_model_def.m
View file @ 5efefaec
% AMDepthDisc_pemo_cfg - example perception model configuration file -
%
% This matlab script is called by sim_main when starting the experiment
% 'AMDepthDisc' and model 'pemo'.
% AMDepthDisc_pemo_cfg outputs the structure 'simdef' containing the complete
% model configuration for the experiment.
%
% See also help AMDepthDisc_set, AMDepthDisc_user, sim_main
%
% Copyright (c) 1999 - 2004 Stephan Ewert.
% $Revision: 1.00.2 beta$ $Date: 30-04-2004 10:49 $
% modified version,
% 24/10/2011, C. Iben
function[simdef] = mreGPSMq_model_def()
%% model variables
simdef.MaxSample = 1.0; % maximal amplitude of input signal - DO NOT TOUCH
simdef.MaxLevel = 100; % simulated level for MaxSample - DO NOT TOUCH
simdef.store_ref = 0; % simulation either with stored references (1) or running references (0)
simdef.templ_num = 16; % number of int. rep. used for deriving the template
simdef.nstore = 16; % number of stored internal ref's (no meaning when running references)
simdef.BMfilter = 'Gamma'; % 'Gamma' or 'DRNL_CASP_2008', 'DRNL'
% mreGPSMq_model_def.m this function sets parameters of the gammatone
% filterbank and the modulation filterbank and turns resampling on/off
% OUTPUT:
% def: struct for stimulus related parameters,e.g.,fs
% simdef: struct for model related parameters, e.g., filter order,
% cut-off frequency
% simwork: struct for parameters applied during simulation, e.g.,
% filter coefficients
%
% Usage: [simdef] = mreGPSMq_model_def()
% thomas.biberger@uni-oldenburg.de;
% date: 2013-03-07
%% Resampling
simdef.resampling='on'; % set resampling 'on' or 'off'
%% Gammatone filter variables (don't mention values below, now we are using fixed values deposited in gfb_analyzer_new_tb.m)
simdef.gt_MinCF = 100; % lowest CF in Hz
simdef.gt_Dens ='terz'; % only 'terz' is a valid option
% *************************************************************************
% if simdef.gt_Dens= 'terz' is selected one ERB wide filters (as default simdef.gt_BW=1) with center frequencies of 63 80 100 125 160 200 250 315 400 500 630 800 1000
% 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500 16000 Hz are applied.
% Don't mind values given below for simdef.gt_MinCF, simdef.gt_MaxCF, simdef.gt_align = 1000;
% *************************************************************************
simdef.gt_MinCF = 63; % lowest CF in Hz
simdef.gt_MaxCF = 16000; % highest CF in Hz (mj init 4000)
simdef.gt_align = 1000; % 15/11/10 CI: specify base frequency around which first filter is aligned
simdef.gt_MultipleCFs = []; % use combinations of low/upp like [100 1000 3000 4000 9000 9000] overrides min/maxCF
simdef.gt_BW = 1.0; % bandwidth of the filter in ERB
simdef.gt_Dens ='terz'; % default:[1] filter density in 1/ERB, or 'terz', 'octave'
% if schroeder == 1
% simdef.gt_useGfbMex = 0; % use 32-bit float precision c/mex implementation of gammatone filterbank instead of matlab implementation
% simdef.gt_useHwrlpmex = 0; % use 32-bit float precision c/mex implementation of 'haircell' envelope extraction
% simdef.gt_Nlalmex = 0; % use 32-bit float precision c/mex implementation of adaptation loops
% else