forked from m-r-s/reference-feature-extraction
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathepsi.m
206 lines (183 loc) · 7.57 KB
/
epsi.m
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
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
function [epsi epsi_std] = epsi(snr1, correct1, total1, snr2, correct2, total2)
% usage: [epsi epsi_std] = epsi(snr1, correct1, total1, snr2, correct2, total2)
% snr Vector of signal-to-noise ratios (SNR) in dB
% correct Vector of number of correct decisions at the corresponding SNRs
% total Total number of decisions (may also be a vector)
% epsi Equal-performance SNR increase (EPSI)
% epsi_std Estimated standard deviation of the EPSI
%
% - Equal-performance SNR increase v1.0 -
%
% This script calculates the increase in signal-to-noise ratio (SNR) which
% is needed on average in order to make System2 perform as well as System1.
% If the equal-performance SNR increase (EPSI) is positive the SNR needs
% to be increased, which would indicate that System2 is less robust than
% System1. A negative EPSI would inidcate that System2 is more robust than
% System1 because the SNR for System2 needs to be decrased to match the
% average recognition performance of System1 and System2.
% A detailed explanation is given in [1].
%
% Copyright (C) 2015 Marc René Schädler
% E-mail [email protected]
% Institute Carl-von-Ossietzky University Oldenburg, Germany
%
%-----------------------------------------------------------------------------
%
% Release Notes:
% v1.0 - Inital release
%
%%%%%%%%%%%%%%%%%%%%%%%%% Example with data from [1] %%%%%%%%%%%%%%%%%%%%%%%%%
% SNR = [ -6 -3 +0 +3 +6 +9]; % SNR in dB
% HSR = [90.30 93.00 93.80 95.30 96.80 98.80]; % percent correct
% MFCC = [68.66 74.58 82.25 87.50 89.08 92.00]; % percent correct
% GBFB = [71.41 77.75 84.25 88.91 92.16 92.66]; % percent correct
% ASR_decisions = 1200; % Number of independent decisions
% HSR_decisions = 1200; % Number of independent decisions
% HSR_correct = HSR./100.*HSR_decisions;
% MFCC_correct = MFCC./100.*ASR_decisions;
% GBFB_correct = GBFB./100.*ASR_decisions;
%
% % Calculate the EPSI of GBFB over HSR
% [HSR_GBFB_epsi HSR_GBFB_epsi_std] = epsi( ...
% SNR, HSR_correct, HSR_decisions, ...
% SNR, GBFB_correct, ASR_decisions);
% [HSR_MFCC_epsi HSR_MFCC_epsi_std] = epsi( ...
% SNR, HSR_correct, HSR_decisions, ...
% SNR, MFCC_correct, ASR_decisions);
% [GBFB_MFCC_epsi GBFB_MFCC_epsi_std] = epsi( ...
% SNR, GBFB_correct, ASR_decisions, ...
% SNR, MFCC_correct, ASR_decisions);
%
% disp(['The SNR needs to be increased by ' ...
% num2str(HSR_GBFB_epsi) '+-' num2str(HSR_GBFB_epsi_std) ' dB' ...
% ' for GBFB in order to make GBFB and HSR perform equally.']);
% disp(['The SNR needs to be increased by ' ...
% num2str(HSR_MFCC_epsi) '+-' num2str(HSR_MFCC_epsi_std) ' dB' ...
% ' for MFCC in order to make MFCC and HSR perform equally.']);
% disp(['The SNR needs to be increased by ' ...
% num2str(GBFB_MFCC_epsi) '+-' num2str(GBFB_MFCC_epsi_std) ' dB' ...
% ' for MFCC in order to make MFCC and GBFB perform equally.']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EPSI implementation %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Number of bootstrap samples for error estimation
num_bootstrap = 1000;
% Make sure all data is stored in column vectors
snr1 = snr1(:);
correct1 = correct1(:);
total1 = total1(:);
snr2 = snr2(:);
correct2 = correct2(:);
total2 = total2(:);
% Expand total numbers if no SNR-dependent vectors are supplied
if numel(total1) == 1
total1 = total1.*ones(size(snr1));
end
if numel(total2) == 1
total2 = total2.*ones(size(snr2));
end
% Sort results by SNR and estimate the uncertainties
[snr1 sortidx1] = sort(snr1);
correct1 = correct1(sortidx1);
total1 = total1(sortidx1);
performance1 = correct1./total1;
uncertainty1 = estd(round(total1), round(correct1), num_bootstrap);
[snr2 sortidx2] = sort(snr2);
correct2 = correct2(sortidx2);
total2 = total2(sortidx2);
performance2 = correct2./total2;
uncertainty2 = estd(round(total2), round(correct2), num_bootstrap);
% Calculate the shift in SNR to match the performance of the two systems
epsi = shiftmetric(snr1, performance1, snr2, performance2);
% Estimate the EPSI standard deviation by bootstrapping and assuming normally
% distributed deviations on the performance scale.
% Please note: The shift is independent of the scale of the y-axis but its
% estimated deviations are not. Ideally the performance would be reported on
% a scale on which the deviations could be expected to be normally distributed
% in order to better match this assumtion.
epsi_std = NaN;
if ~isnan(epsi)
bootstrap = zeros(num_bootstrap,1);
for k=1:num_bootstrap
performance1_tmp = performance1 + randn(size(performance1)).*uncertainty1;
performance2_tmp = performance2 + randn(size(performance2)).*uncertainty2;
bootstrap(k) = shiftmetric(snr1, performance1_tmp, snr2, performance2_tmp);
end
if any(~isfinite(bootstrap))
warning('Ignoring infinite numbers during error estimation!');
end
epsi_std = std(bootstrap(isfinite(bootstrap)));
end
end
function shift = shiftmetric(x1, y1, x2, y2)
% Calulcates the average shift of two graphs, y1(x1) and y2(x2), on the x-axis.
% The average is taken in steps of 0.5 over the x-range in which both graphs
% are defined. Hence, the calculated shift does not depend on the scale of the
% y-axis.
% Make sure that data comes in column vectors
x1 = x1(:);
y1 = y1(:);
x2 = x2(:);
y2 = y2(:);
% Make sure the performance graphs are monotonic (c.f. Equation 5 in [1])
y1 = monotonify(-y1(end:-1:1),0.0001);
y2 = monotonify(-y2(end:-1:1),0.0001);
y1 = -y1(end:-1:1);
y2 = -y2(end:-1:1);
shift = NaN;
assert(all(size(x1) == size(y1)) && all(size(x2) == size(y2)), ...
'x and y vectors must have same lengths');
% Get the overlapping window on the y-axis.
% This is the region in which the two graphs can be compared.
yrange = [max(min(y1),min(y2)); min(max(y1),max(y2))];
% Interpolate the corrsponding ranges on the x-axis.
xrange1 = interp1q(y1, x1, yrange);
xrange2 = interp1q(y2, x2, yrange);
% Sample the x-ranges in 0.5 dB steps
xrange1 = [ceil(xrange1(1)*2)/2; floor(xrange1(2)*2)/2];
xpoints1 = (xrange1(1):0.5:xrange1(2)).';
xrange2 = [ceil(xrange2(1)*2)/2; floor(xrange2(2)*2)/2];
xpoints2 = (xrange2(1):0.5:xrange2(2)).';
% Integrate the shifts in both directions if sampling points exist
if ~isempty(xpoints1) && ~isempty(xpoints2)
% Look up the y-values at the sample points
ypoints1 = interp1q(x1, y1, xpoints1);
ypoints2 = interp1q(x2, y2, xpoints2);
% Get the x-values of the "other" system at the same y-value
xpoints12 = interp1q(y2, x2, ypoints1);
xpoints21 = interp1q(y1, x1, ypoints2);
% Average over the difference on the x-values
shift1 = mean(xpoints12-xpoints1);
shift2 = mean(xpoints21-xpoints2);
% Weight the calculated differences from both graphs equally
shift = 0.5.*(shift1-shift2);
end
end
function y = monotonify(x, epsilon)
% Guarantees monotonically increasing output
y = x;
if nargin < 2
epsilon = eps;
end
for i=2:length(x)
y(i) = max(x(i-1)+epsilon, y(i));
end
end
function rate_std = estd(total, correct, n)
% Estimates the standard deviation of the correct-rate of a finite number
% of binary decisions by bootstrapping
if numel(total) > 1
% Handle the situation when total is a vector
if numel(n) == 1
n = n.*ones(size(total));
end
rate_std = arrayfun(@estd, total, correct, n);
else
% Generate a vector of ones and zeros reflecting correct and false decisions
x = [zeros(1,total-correct) ones(1,correct)];
% Generate n random selections drawing total samples from x
selections = ceil(rand(total,n)*total);
% Get the rates for these selections
rates = mean(x(selections));
% Get the standard deviation of the artificially generated samples
rate_std = std(rates);
end
end