BAMQFrontEnd.m 4.58 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
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;