1 /*-------------- Telecommunications & Signal Processing Lab ---------------
2 McGill University
3
4 Routine:
5 struct MOVB_PD PQmovPD (const double *Ehs[])
6 void PQchanPD (int Nc, int Nchan, const struct PQ_MOVBI MOVBI[2],
7 double *Pc[], double *Qc)
8 void PQavgPD (const struct MOVC_PD *PD, int Np,
9 double *ADBB, double *MFPDB)
10
11 Purpose:
12 Calculate the probability of detection values
13 Average the probability of detection values
14
15 Description:
16 PQmovPD:
17 This routine calculates the probability of detection for each frequency for a
18 frame of data.
19
20 PQchanPD:
21 This routine compares probability of detection values across channels.
22
23 PQavgPD:
24 This routine average the probability of detection values to form the model
25 output variables.
26
27 Parameters:
28 PQmovPD:
29 <- struct MOVB_PD PQmovPD
30 Probability of detection values
31 -> const double *Ehs[]
32 Excitation patterns
33
34 PQchanPD:
35 -> int Nc
36 Number of bands
37 -> int Nchan
38 Number of channels
39 -> const struct PQ_MOVBI MOVBI[2]
40 Input values for each channel
41 <- double *Pc
42 Output probability of detection value (max. over channels)
43 <- double *Qc
44 Output total number of steps above the threshold
45
46 PQavgPD:
47 -> const struct MOVC_PD *PD
48 Probability of detection model output variable precursors
49 -> int Np
50 Number of PEAQ frames
51 <- double *ADBB
52 Output average distorted block MOV
53 <- double *MFPDB
54 Maximum filtered probability of detection MOV
55
56 Author / revision:
57 P. Kabal Copyright (C) 2003
58 $Revision: 1.7 $ $Date: 2003/05/13 01:15:13 $
59
60 -------------------------------------------------------------------------*/
61
62 #include <math.h>
63
64 #include "PQevalAudio.h"
65
66 #define MAXV(a, b) (((a) > (b)) ? (a) : (b))
67 #define ABSV(x) (((x) < 0) ? -(x) : (x))
68 #define DBP(x) (10.0 * log10 (x))
69
70
71 struct MOVB_PD
PQmovPD(const double * Ehs[])72 PQmovPD (const double *Ehs[])
73
74 {
75 /* Polynomial approximation parameters */
76 static const double c[] =
77 {-0.198719, 0.0550197, -0.00102438, 5.05622e-6, 9.01033e-11};
78 static const double d1 = 5.95072;
79 static const double d2 = 6.39468;
80 static const double g = 1.71332;
81
82 /* Slope parameters */
83 static const double bP = 4;
84 static const double bM = 6;
85
86 static const int Nc = PQ_NC_B;
87
88 int m;
89 double EdBR, EdBT, edB;
90 double L, b, s;
91 struct MOVB_PD PDV;
92
93 for (m = 0; m < Nc; ++m) {
94 EdBR = DBP (Ehs[0][m]);
95 EdBT = DBP (Ehs[1][m]);
96 edB = EdBR - EdBT;
97 if (edB > 0) {
98 L = 0.3 * EdBR + 0.7 * EdBT;
99 b = bP;
100 }
101 else {
102 L = EdBT;
103 b = bM;
104 }
105 if (L > 0)
106 s = d1 * pow (d2 / L, g)
107 + c[0] + L * (c[1] + L * (c[2] + L * (c[3] + L * c[4])));
108 else
109 s = 1e30;
110 PDV.p[m] = 1 - pow (0.5, pow (edB / s, b)); /* Detection probability */
111 PDV.q[m] = ABSV ((int) (edB)) / s; /* Steps above threshold */
112 }
113
114 return PDV;
115 }
116
117 /* Compare prob. detection values across channels */
118
119
120 void
PQchanPD(int Nc,int Nchan,const struct PQ_MOVBI MOVBI[2],double * Pc,double * Qc)121 PQchanPD (int Nc, int Nchan, const struct PQ_MOVBI MOVBI[2],
122 double *Pc, double *Qc)
123
124 {
125 int m;
126 double Pr, Qcx, pbin, qbin;
127 const struct MOVB_PD *PD1;
128 const struct MOVB_PD *PD2;
129
130
131 Pr = 1;
132 Qcx = 0;
133 if (Nchan > 1) {
134 PD1 = &MOVBI[0].PDV;
135 PD2 = &MOVBI[1].PDV;
136 for (m = 0; m < Nc; ++m) {
137 pbin = MAXV (PD1->p[m], PD2->p[m]);
138 qbin = MAXV (PD1->q[m], PD2->q[m]);
139 Pr = Pr * (1 - pbin);
140 Qcx += qbin;
141 }
142 }
143 else {
144 PD1 = &MOVBI[0].PDV;
145 for (m = 0; m < Nc; ++m) {
146 Pr = Pr * (1 - PD1->p[m]);
147 Qcx += PD1->q[m];
148 }
149 }
150
151 *Pc = 1 - Pr;
152 *Qc = Qcx;
153
154 return;
155 }
156
157 /* Average probability of detection values */
158
159
160 void
PQavgPD(const struct MOVC_PD * PD,int Np,double * ADBB,double * MFPDB)161 PQavgPD (const struct MOVC_PD *PD, int Np, double *ADBB, double *MFPDB)
162
163 {
164 int i, Nd;
165 double Phc, Pcmax, Qsum;
166
167 static const double c0 = 0.9;
168 const double c1 = PD->c1;
169
170 Phc = 0;
171 Pcmax = 0;
172 Qsum = 0;
173 Nd = 0;
174 for (i = 0; i < Np; ++i) {
175 Phc = c0 * Phc + (1 - c0) * PD->Pc[i];
176 Pcmax = MAXV (Pcmax * c1, Phc);
177
178 if (PD->Pc[i] > 0.5) {
179 ++Nd;
180 Qsum += PD->Qc[i];
181 }
182 }
183
184 if (Nd == 0)
185 *ADBB = 0;
186 else if (Qsum > 0)
187 *ADBB = log10 (Qsum / Nd);
188 else
189 *ADBB = -0.5;
190
191 *MFPDB = Pcmax;
192
193 return;
194 }
195