1 /*-------------- Telecommunications & Signal Processing Lab ---------------
2 McGill University
3
4 Routine:
5 double PQmovEHS (const double *X2[], const struct Par_EHS *EHS)
6 double PQavgEHS (const struct MOVC_EHS *EHS, int Nchan, int Np)
7
8 Purpose:
9 Calculate the error harmonic structure parameter
10 Average the error harmonic structure values
11
12 Description:
13 PQmovEHS:
14 This routine calculates the error harmonic structure parameter for a frame
15 of data.
16
17 PQavbEHS:
18 Average the per frame EHS values to give the final model output variable.
19
20 Parameters:
21 PQmovEHS:
22 <- double PQmovEHS
23 Maximum value of the correlation power spectrum
24 -> const double xR[]
25 Reference time samples
26 -> const double xT[]
27 Test time samples
28 -> const double *X2[]
29 DFT energies
30 -> const struct Par_EHS *EHS
31 Structure with EHS parameters
32
33 PQavgEHS:
34 <- double PQavgEHS
35 EHS model output variable
36 -> const struct MOVC_EHS *EHS
37 Input EHS values
38 -> int Nchan
39 Number of channels
40 -> int Np
41 Number of PEAQ frames
42
43 Author / revision:
44 P. Kabal Copyright (C) 2003
45 $Revision: 1.9 $ $Date: 2003/05/13 01:15:12 $
46
47 -------------------------------------------------------------------------*/
48
49 #include <assert.h>
50 #include <math.h>
51
52 #include <libtsp.h>
53 #include "PQevalAudio.h"
54
55 #define SQRV(x) ((x) * (x))
56 #define KST_MAX 1
57 #define NLAG 256
58 #define NSUM 256
59 #define ND (KST_MAX + NLAG + NSUM - 1)
60 #define SUMSQ(x,N) VRdDotProd(x,x,N)
61
62 static void
63 PQ_Corr (const double D0[], const double D1[], double C[], int NL, int M);
64 static void
65 PQ_Demean (const double x[], double y[], int N);
66 static double
67 PQ_FindPeak (const double c2[], int N);
68 static double
69 PQ_LinPosAvg (const double x[], int N);
70 static void
71 PQ_NCorr (const double C[], const double D[], double Cn[], int NL, int M);
72
73
74 double
PQmovEHS(const double xR[],const double xT[],const double * X2[],const struct Par_EHS * EHS)75 PQmovEHS (const double xR[], const double xT[], const double *X2[],
76 const struct Par_EHS *EHS)
77
78 {
79 int k, kmax;
80 double D[ND];
81 double C[KST_MAX + NLAG];
82 double Cm[NLAG];
83 double c2[NLAG/2 + 1];
84 double EnRef, EnTest;
85 double EHSV;
86
87 const int kst = EHS->kst;
88 const int NL = EHS->NL;
89 const int M = EHS->M;
90 const double *Hw = EHS->Hw;
91
92 /* Energy in the second half frame */
93 EnRef = SUMSQ (&xR[PQ_CB_ADV], PQ_NF - PQ_CB_ADV);
94 EnTest = SUMSQ (&xT[PQ_CB_ADV], PQ_NF - PQ_CB_ADV);
95
96 /* Set the return value to be negative for small energy frames */
97 if (EnRef < EHS->EnThr && EnTest < EHS->EnThr)
98 return -1.;
99
100 kmax = kst + NL - 1 + M;
101 assert (ND >= kmax && NLAG >= NL && KST_MAX +NLAG >= kst + NL);
102
103 /* Differences of log values */
104 for (k = 0; k < kmax; ++k)
105 D[k] = log (X2[1][k] / X2[0][k]);
106
107 /* Compute correlation */
108 PQ_Corr (D, D, C, kst+NL, M);
109
110 /* Normalize the correlations, remove the mean */
111 PQ_NCorr (C, D, C, kst+NL, M);
112 PQ_Demean (&C[kst], Cm, NL);
113
114 /* Window the correlation */
115 VRdMult (Cm, Hw, Cm, NL);
116
117 /* DFT */
118 SPdRFFT (Cm, NL, 1);
119
120 /* Squared magnitude */
121 VRdRFFTMSq (Cm, c2, NL);
122
123 /* Search for a peak after a valley */
124 EHSV = PQ_FindPeak (c2, NL/2+1);
125
126 return EHSV;
127 }
128
129 /* Average EHS values (negative values flag low energy frames) */
130
131
132 double
PQavgEHS(const struct MOVC_EHS * EHS,int Nchan,int Np)133 PQavgEHS (const struct MOVC_EHS *EHS, int Nchan, int Np)
134
135 {
136 int j;
137 double s, EHSB;
138
139 s = 0;
140 for (j = 0; j < Nchan; ++j)
141 s += PQ_LinPosAvg (EHS->EHS[j], Np);
142 EHSB = 1000 * s / Nchan;
143
144 return EHSB;
145 }
146 static void
PQ_Demean(const double x[],double y[],int N)147 PQ_Demean (const double x[], double y[], int N)
148
149 {
150 int i;
151 double s;
152
153 s = 0;
154 for (i = 0; i < N; ++i)
155 s += x[i];
156 s = s / N;
157
158 for (i = 0; i < N; ++i)
159 y[i] = x[i] - s;
160
161 return;
162 }
163
164 /* Correlation calculation */
165
166 #define NFFT 512
167
168
169 static void
PQ_Corr(const double D0[],const double D1[],double C[],int NL,int M)170 PQ_Corr (const double D0[], const double D1[], double C[], int NL, int M)
171
172 {
173 int n, m;
174 double D0x[NFFT];
175 double D1x[NFFT];
176 double dx[NFFT];
177
178 assert (NFFT >= NL + M - 1);
179
180 /* Calculate the correlation using DFT's */
181 VRdCopy (D0, D0x, M);
182 VRdZero (&D0x[M], NFFT-M); /* Zero padding */
183 VRdCopy (D1, D1x, M+NL-1);
184 VRdZero (&D1x[M+NL-1], NFFT-(M+NL-1));
185
186 /* DFTs of the zero-padded sequences */
187 SPdRFFT (D0x, NFFT, 1);
188 SPdRFFT (D1x, NFFT, 1);
189
190 /* Multiply the (complex) DFT sequences (D0x is conjugated) */
191 dx[0] = D0x[0] * D1x[0];
192 for (n = 1; n < NFFT/2; ++n) {
193 m = NFFT/2 + n; /* n => real part, m => imaginary part */
194 dx[n] = D0x[n] * D1x[n] + D0x[m] * D1x[m];
195 dx[m] = D0x[n] * D1x[m] - D0x[m] * D1x[n];
196 }
197 dx[NFFT/2] = D0x[NFFT/2] * D1x[NFFT/2];
198
199 /* Inverse DFT */
200 SPdRFFT (dx, NFFT, -1);
201
202 /* Return the first NL elements */
203 VRdCopy (dx, C, NL);
204
205 return;
206 }
207
208 /* Normalize the correlation values */
209
210 static void
PQ_NCorr(const double C[],const double D[],double Cn[],int NL,int M)211 PQ_NCorr (const double C[], const double D[], double Cn[], int NL, int M)
212
213 {
214 int i;
215 double s0, sj, d;
216
217 s0 = C[0];
218 sj = s0;
219 Cn[0] = 1;
220
221 for (i = 1; i < NL; ++i) {
222 sj += (SQRV (D[i+M-1]) - SQRV (D[i-1]));
223 d = s0 * sj;
224 if (d <= 0)
225 Cn[i] = 1;
226 else
227 Cn[i] = C[i] / sqrt (d);
228 }
229
230 return;
231 }
232
233 /* Find the peak of the correlation power spectrum */
234
235 static double
PQ_FindPeak(const double c2[],int N)236 PQ_FindPeak (const double c2[], int N)
237
238 {
239 int n;
240 double cprev, cmax;
241
242 assert (N >= 2);
243
244 cprev = c2[0];
245 cmax = 0;
246
247 /* Search for a peak after a valley */
248 for (n = 1; n < N; ++n) {
249 if (c2[n] > cprev) { /* Rising from a valley */
250 if (c2[n] > cmax)
251 cmax = c2[n];
252 }
253 }
254
255 return cmax;
256 }
257
258 /* Average values, omitting values which are negative */
259
260 static double
PQ_LinPosAvg(const double x[],int N)261 PQ_LinPosAvg (const double x[], int N)
262
263 {
264 int i, Nv;
265 double s;
266
267 Nv = 0;
268 s = 0;
269 for (i = 0; i < N; ++i) {
270 if (x[i] >= 0) {
271 s += x[i];
272 ++Nv;
273 }
274 }
275
276 if (Nv > 0)
277 s = s / Nv;
278
279 return s;
280 }
281