1 /*-------------- Telecommunications & Signal Processing Lab ---------------
2 McGill University
3
4 Routine:
5 struct PQ_MOVBI PQeval (const double xR[], const double xT[],
6 const struct PQ_Par *PQpar, struct PQ_FiltMem *Fmem)
7
8 Purpose:
9 Calculate PEAQ FFT model values for one frame (PEAQ Basic version)
10
11 Description:
12 This routine calculates the PEAQ model output variable precursors for one
13 frame (PEAQ Basic version).
14
15 Parameters:
16 <- struct PQ_MOVBI PQeval
17 Model output variable precursors
18 -> const double xR[]
19 Reference signal
20 -> const double xT[]
21 Test signal
22 -> const struct PQ_Par *PQpar
23 Tables and parameters
24 <-> struct PQ_FiltMem *Fmem
25 Filter parameters and memory
26
27 Author / revision:
28 P. Kabal Copyright (C) 2003
29 $Revision: 1.14 $ $Date: 2003/05/13 01:10:26 $
30
31 -------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <math.h>
35 #include <stdio.h>
36
37 #include <libtsp.h>
38 #include "PQevalAudio.h"
39
40 static void
41 PQ_DFTFrame (const double x[], const double hw[], double X2[], int NF);
42 static void
43 PQ_excitCB (const double *X2[2], double EbN[], double *Es[2],
44 const struct Par_CB *CB);
45
46
47 struct PQ_MOVBI
PQeval(const double xR[],const double xT[],const struct PQ_Par * PQpar,struct PQ_FiltMem * Fmem)48 PQeval (const double xR[], const double xT[], const struct PQ_Par *PQpar,
49 struct PQ_FiltMem *Fmem)
50
51 {
52 double xFFT[2 * PQ_NF];
53 double *X2[2] = {&xFFT[0], &xFFT[PQ_NF]};
54 double xCB[6 * PQ_NC_B];
55 double *Es[2] = {&xCB[0*PQ_NC_B], &xCB[PQ_NC_B]};
56 double *Ehs[2] = {&xCB[2*PQ_NC_B], &xCB[3*PQ_NC_B]};
57 double *EP[2] = {&xCB[4*PQ_NC_B], &xCB[5*PQ_NC_B]};
58 double xMod[2 * PQ_MAXNC];
59 double *M[2] = {&xMod[0], &xMod[PQ_MAXNC]};
60 double ERavg[PQ_MAXNC];
61
62 double EbN[PQ_NC_B];
63 struct PQ_MOVBI MOVB;
64
65 static const int NF = PQ_NF;
66
67 /* Windowed DFT */
68 PQ_DFTFrame (xR, PQpar->FFT.hw, X2[0], NF);
69 PQ_DFTFrame (xT, PQpar->FFT.hw, X2[1], NF);
70
71 /* Critical band grouping and frequency spreading */
72 PQ_excitCB ((const double **) X2, EbN, Es, &PQpar->CB);
73
74 /* Time domain smoothing => "Excitation patterns" */
75 PQtimeSpread (Es[0], &PQpar->TDS, Ehs[0], Fmem->TDS.x[0]);
76 PQtimeSpread (Es[1], &PQpar->TDS, Ehs[1], Fmem->TDS.x[1]);
77
78 /* Level and pattern adaptation => "Spectrally adapted patterns" */
79 PQadapt ((const double **) Ehs, &PQpar->Patt, EP, &Fmem->Adap);
80
81 /* Modulation patterns */
82 PQmodPatt ((const double **) Es, &PQpar->MPatt, M, ERavg, &Fmem->Env);
83
84 /* Loudness calculation */
85 MOVB.LoudV.NRef = PQloud (Ehs[0], &PQpar->Loud);
86 MOVB.LoudV.NTest = PQloud (Ehs[1], &PQpar->Loud);
87
88 /* Modulation differences */
89 MOVB.MDiffV = PQmovModDiffB ((const double **) M, ERavg, &PQpar->MDiff);
90
91 /* Noise Loudness */
92 MOVB.NLoudV.NL = PQmovNLoudB ((const double **) M, (const double **) EP,
93 &PQpar->NLoud);
94
95 /* Bandwidth */
96 MOVB.BWV = PQmovBW ((const double **) X2, &PQpar->BW);
97
98 /* Noise-to-mask ratios */
99 MOVB.NMRV = PQmovNMRB (EbN, Ehs[0], &PQpar->NMR);
100
101 /* Probability of detection */
102 MOVB.PDV = PQmovPD ((const double **) Ehs);
103
104 /* Error harmonic structure */
105 MOVB.EHSV.EHS = PQmovEHS (xR, xT, (const double **) X2, &PQpar->EHS);
106
107 return MOVB;
108 }
109 /* Process a frame of time samples, returning a squared magnitude spectrum */
110
111 static void
PQ_DFTFrame(const double x[],const double hw[],double X2[],int NF)112 PQ_DFTFrame (const double x[], const double hw[], double X2[], int NF)
113
114 {
115 double X[PQ_NF];
116
117 assert (NF == PQ_NF);
118
119 /* Window the data */
120 VRdMult (x, hw, X, NF);
121
122 /* Real FFT */
123 SPdRFFT (X, NF, +1);
124
125 /* Squared magnitude */
126 VRdRFFTMSq (X, X2, NF);
127
128 return;
129 }
130
131 /* DFT values to critical band values */
132
133 static void
PQ_excitCB(const double * X2[2],double EbN[],double * Es[2],const struct Par_CB * CB)134 PQ_excitCB (const double *X2[2], double EbN[], double *Es[2],
135 const struct Par_CB *CB)
136
137 {
138 int k;
139 double Xw2[2][PQ_NF/2+1];
140 double XwN2[PQ_NF/2+1];
141 double E[2][PQ_NC_B];
142 double Eb[2][PQ_NC_B];
143
144 const int Nc = CB->Nc;
145 static const int NF = PQ_NF;
146
147 /* Outer and middle ear filtering */
148 VRdMult (X2[0], CB->W2, Xw2[0], NF/2+1);
149 VRdMult (X2[1], CB->W2, Xw2[1], NF/2+1);
150
151 /* Form the difference magnitude signal */
152 for (k = 0; k <= NF/2; ++k)
153 XwN2[k] = Xw2[0][k] - 2 * sqrt (Xw2[0][k] * Xw2[1][k]) + Xw2[1][k];
154
155 /* Group into partial critical bands */
156 PQgroupCB (Xw2[0], CB, Eb[0]);
157 PQgroupCB (Xw2[1], CB, Eb[1]);
158 PQgroupCB (XwN2, CB, EbN); /* Noise patterns */
159
160 /* Add internal noise => "Pitch patterns" */
161 VRdAdd (Eb[0], CB->EIN, E[0], Nc);
162 VRdAdd (Eb[1], CB->EIN, E[1], Nc);
163
164 /* Critical band spreading => "Unsmeared excitation patterns" */
165 PQspreadCB (E[0], CB, Es[0]);
166 PQspreadCB (E[1], CB, Es[1]);
167
168 return;
169 }
170