1 /* Sound_to_Harmonicity_GNE.cpp
2 *
3 * Copyright (C) 1999-2012,2015-2021 Paul Boersma
4 *
5 * This code is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This code is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 * See the GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this work. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 /* a replication of:
20 D. Michaelis, T. Gramss & H.W. Strube (1997):
21 "Glottal-to-noise excitation ratio -- a new measure
22 for describing pathological voices."
23 ACUSTICA - acta acustica 83: 700-706.
24 henceforth abbreviated as "MGS".
25 */
26
27 #include "Sound_to_Harmonicity.h"
28 #include "Sound_and_LPC.h"
29 #include "Sound_and_Spectrum.h"
30
bandFilter(Spectrum me,double fmid,double bandwidth)31 static void bandFilter (Spectrum me, double fmid, double bandwidth) {
32 double *re = & my z [1] [0], *im = & my z [2] [0];
33 double fmin = fmid - bandwidth / 2.0, fmax = fmid + bandwidth / 2.0;
34 double twopibybandwidth = 2.0 * NUMpi / bandwidth;
35 for (integer col = 1; col <= my nx; col ++) {
36 double x = my x1 + (col - 1) * my dx;
37 if (x < fmin || x > fmax) {
38 re [col] = 0.0;
39 im [col] = 0.0;
40 } else {
41 double factor = 0.5 + 0.5 * cos (twopibybandwidth * (x - fmid));
42 re [col] *= factor;
43 im [col] *= factor;
44 }
45 }
46 }
47
Sound_to_Harmonicity_GNE(Sound me,double fmin,double fmax,double bandwidth,double step)48 autoMatrix Sound_to_Harmonicity_GNE (Sound me,
49 double fmin, // 500 Hz
50 double fmax, // 4500 Hz
51 double bandwidth, // 1000 Hz
52 double step) // 80 Hz
53 {
54 try {
55 autoSound envelope [1+100];
56 integer nenvelopes = Melder_ifloor ((fmax - fmin) / step);
57 for (integer ienvelope = 1; ienvelope <= 100; ienvelope ++)
58 Melder_assert (! envelope [ienvelope].get());
59
60 /*
61 * Step 1: down-sampling to 10 kHz,
62 * in order to be able to flatten the spectrum
63 * (since the human voice does not contain much above 5 kHz).
64 */
65 autoSound original10k = Sound_resample (me, 10000, 500);
66 Vector_subtractMean (original10k.get());
67 double duration = my xmax - my xmin;
68
69 /*
70 * Step 2: inverse filtering of the speech signal
71 * by 13th-order "autocorrelation method"
72 * with a Gaussian (not Hann, like MGS!) window of 30 ms length
73 * and 10 ms shift between successive frames.
74 * Since we need a spectrally flat signal (not an approximation
75 * of the source signal), we must turn the pre-emphasis off
76 * (by setting its turnover point at 1,000,000,000 Hz);
77 * otherwise, the pre-emphasis would cause an overestimation
78 * in the LPC object of the high frequencies, so that inverse
79 * filtering would yield weakened high frequencies.
80 */
81 autoLPC lpc = Sound_to_LPC_autocorrelation (original10k.get(), 13, 30e-3, 10e-3, 1e9);
82 autoSound flat = LPC_Sound_filterInverse (lpc.get(), original10k.get());
83 autoSpectrum flatSpectrum = Sound_to_Spectrum (flat.get(), true);
84 autoSpectrum hilbertSpectrum = Data_copy (flatSpectrum.get());
85 for (integer col = 1; col <= hilbertSpectrum -> nx; col ++) {
86 hilbertSpectrum -> z [1] [col] = flatSpectrum -> z [2] [col];
87 hilbertSpectrum -> z [2] [col] = - flatSpectrum -> z [1] [col];
88 }
89 double fmid = fmin;
90 integer ienvelope = 1;
91 autoMelderMonitor monitor (U"Computing Hilbert envelopes...");
92 while (fmid <= fmax) {
93 /*
94 * Step 3: calculate Hilbert envelopes of bands.
95 */
96 autoSpectrum bandSpectrum = Data_copy (flatSpectrum.get());
97 autoSpectrum hilbertBandSpectrum = Data_copy (hilbertSpectrum.get());
98 /*
99 * 3a: Filter both the spectrum of the original flat sound and its Hilbert transform.
100 */
101 bandFilter (bandSpectrum.get(), fmid, bandwidth);
102 bandFilter (hilbertBandSpectrum.get(), fmid, bandwidth);
103 /*
104 * 3b: Create both the band-filtered flat sound and its Hilbert transform.
105 */
106 autoSound band = Spectrum_to_Sound (bandSpectrum.get());
107 /*if (graphics) {
108 Graphics_beginMovieFrame (graphics, & Melder_WHITE);
109 Spectrum_draw (bandSpectrum, graphics, 0, 5000, 0, 0, true);
110 Graphics_endMovieFrame (graphics, 0.0);
111 }*/
112 Melder_monitor (ienvelope / (nenvelopes + 1.0), U"Computing Hilbert envelope ", ienvelope, U"...");
113 autoSound hilbertBand = Spectrum_to_Sound (hilbertBandSpectrum.get());
114 envelope [ienvelope] = Sound_extractPart (band.get(), 0, duration, kSound_windowShape::RECTANGULAR, 1.0, true);
115 /*
116 * 3c: Compute the Hilbert envelope of the band-passed flat signal.
117 */
118 for (integer col = 1; col <= envelope [ienvelope] -> nx; col ++) {
119 double self = envelope [ienvelope] -> z [1] [col], other = hilbertBand -> z [1] [col];
120 envelope [ienvelope] -> z [1] [col] = hypot (self, other);
121 }
122 Vector_subtractMean (envelope [ienvelope].get());
123 /*
124 * Next band.
125 */
126 fmid += step;
127 ienvelope += 1;
128 }
129
130 /*
131 * Step 4: crosscorrelation
132 */
133 nenvelopes = ienvelope - 1;
134 autoMatrix cc = Matrix_createSimple (nenvelopes, nenvelopes);
135 for (integer row = 2; row <= nenvelopes; row ++) {
136 for (integer col = 1; col <= row - 1; col ++) {
137 autoSound crossCorrelation = Sounds_crossCorrelate_short (envelope [row].get(), envelope [col].get(), -3.1e-4, 3.1e-4, true);
138 /*
139 * Step 5: the maximum of each correlation function
140 */
141 double ccmax = Vector_getMaximum (crossCorrelation.get(), 0.0, 0.0, kVector_peakInterpolation :: NONE);
142 cc -> z [row] [col] = ccmax;
143 }
144 }
145
146 /*
147 * Step 6: maximum of the maxima, ignoring those too close to the diagonal.
148 */
149 for (integer row = 2; row <= nenvelopes; row ++) {
150 for (integer col = 1; col <= row - 1; col ++) {
151 if (integer_abs (row - col) < bandwidth / 2.0 / step) {
152 cc -> z [row] [col] = 0.0;
153 }
154 }
155 }
156
157 return cc;
158 } catch (MelderError) {
159 Melder_throw (me, U": not converted to Harmonicity (GNE).");
160 }
161 }
162
163 /* End of file Sound_to_Harmonicity_GNE.cpp */
164