1 /* Sound_to_Intensity.cpp
2 *
3 * Copyright (C) 1992-2012,2014-2020 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 /*
20 * pb 2002/07/16 GPL
21 * pb 2003/05/20 default time step is four times oversampling
22 * pb 2003/07/10 NUMbessel_i0_f
23 * pb 2003/11/19 Sound_to_Intensity veryAccurate
24 * pb 2003/12/15 removed bug introduced by previous change
25 * pb 2004/10/27 subtractMean
26 * pb 2006/12/31 compatible with stereo sounds
27 * pb 2007/01/27 for stereo sounds, add channel energies
28 * pb 2007/02/14 honoured precondition of Sampled_shortTermAnalysis (by checking whether minimumPitch is defined)
29 * pb 2008/01/19 double
30 * pb 2011/03/04 C++
31 * pb 2011/03/28 C++
32 */
33
34 #include "Sound_to_Intensity.h"
35
Sound_to_Intensity_(Sound me,double minimumPitch,double timeStep,bool subtractMeanPressure)36 static autoIntensity Sound_to_Intensity_ (Sound me, double minimumPitch, double timeStep, bool subtractMeanPressure) {
37 try {
38 /*
39 Preconditions.
40 */
41 Melder_require (isdefined (minimumPitch),
42 U"Minimum pitch is undefined.");
43 Melder_require (isdefined (timeStep),
44 U"Time step is undefined.");
45 Melder_require (timeStep >= 0.0,
46 U"Time step should be zero (= automatic) or positive, instead of ", timeStep, U" seconds.");
47 Melder_require (my dx > 0.0,
48 U"The Sound's time step should be positive, instead of ", my dx, U" seconds.");
49 Melder_require (minimumPitch > 0.0,
50 U"Minimum pitch should be positive, instead of ", minimumPitch, U" Hz.");
51 /*
52 Defaults.
53 */
54 constexpr double minimumNumberOfPeriodsNeededForReliablePitchMeasurement = 3.2;
55 const double maximumPeriod = 1.0 / minimumPitch;
56 const double logicalWindowDuration = minimumNumberOfPeriodsNeededForReliablePitchMeasurement * maximumPeriod; // == 3.2 / minimumPitch
57 if (timeStep == 0.0) {
58 constexpr double defaultOversampling = 4.0;
59 timeStep = logicalWindowDuration / defaultOversampling; // == 0.8 / minimumPitch
60 }
61
62 const double physicalWindowDuration = 2.0 * logicalWindowDuration; // == 6.4 / minimumPitch
63 Melder_assert (physicalWindowDuration > 0.0);
64 const double halfWindowDuration = 0.5 * physicalWindowDuration;
65 const integer halfWindowSamples = Melder_ifloor (halfWindowDuration / my dx);
66 const integer windowNumberOfSamples = 2 * halfWindowSamples + 1;
67 autoVEC amplitude = zero_VEC (windowNumberOfSamples);
68 autoVEC window = zero_VEC (windowNumberOfSamples);
69 const integer windowCentreSampleNumber = halfWindowSamples + 1;
70
71 for (integer i = 1; i <= windowNumberOfSamples; i ++) {
72 const double x = (i - windowCentreSampleNumber) * my dx / halfWindowDuration;
73 const double root = sqrt (Melder_clippedLeft (0.0, 1.0 - sqr (x))); // clipping should be rare
74 window [i] = NUMbessel_i0_f ((2.0 * NUMpi * NUMpi + 0.5) * root);
75 }
76
77 integer numberOfFrames;
78 double thyFirstTime;
79 try {
80 Sampled_shortTermAnalysis (me, physicalWindowDuration, timeStep, & numberOfFrames, & thyFirstTime);
81 } catch (MelderError) {
82 const double physicalSoundDuration = my nx * my dx;
83 Melder_throw (U"The physical duration of the sound (the number of samples times the sampling period) in an intensity analysis "
84 "should be at least 6.4 divided by the minimum pitch (", minimumPitch, U" Hz), "
85 U"i.e. at least ", physicalWindowDuration, U" s, instead of ", physicalSoundDuration, U" s.");
86 }
87 autoIntensity thee = Intensity_create (my xmin, my xmax, numberOfFrames, timeStep, thyFirstTime);
88 for (integer iframe = 1; iframe <= numberOfFrames; iframe ++) {
89 const double midTime = Sampled_indexToX (thee.get(), iframe);
90 const integer soundCentreSampleNumber = Sampled_xToNearestIndex (me, midTime); // time accuracy is half a sampling period
91
92 integer leftSample = soundCentreSampleNumber - halfWindowSamples;
93 integer rightSample = soundCentreSampleNumber + halfWindowSamples;
94 /*
95 Catch some edge cases, which are uncommon because Sampled_shortTermAnalysis() filtered out most problems.
96 */
97 Melder_clipLeft (1_integer, & leftSample);
98 Melder_clipRight (& rightSample, my nx);
99 Melder_require (rightSample >= leftSample,
100 U"Unexpected edge case: right sample (", rightSample, U") less than left sample (", leftSample, U").");
101
102 const integer windowFromSoundOffset = windowCentreSampleNumber - soundCentreSampleNumber;
103 VEC amplitudePart = amplitude.part (windowFromSoundOffset + leftSample, windowFromSoundOffset + rightSample);
104 constVEC windowPart = window.part (windowFromSoundOffset + leftSample, windowFromSoundOffset + rightSample);
105 longdouble sumxw = 0.0, sumw = 0.0;
106 for (integer ichan = 1; ichan <= my ny; ichan ++) {
107 amplitudePart <<= my z [ichan].part (leftSample, rightSample);
108 if (subtractMeanPressure)
109 centre_VEC_inout (amplitudePart);
110 for (integer isamp = 1; isamp <= amplitudePart.size; isamp ++) {
111 sumxw += sqr (amplitudePart [isamp]) * windowPart [isamp];
112 sumw += windowPart [isamp];
113 }
114 }
115 const double intensity_in_Pa2 = double (sumxw / sumw);
116 constexpr double hearingThreshold_in_Pa = 2.0e-5;
117 constexpr double hearingThreshold_in_Pa2 = sqr (hearingThreshold_in_Pa);
118 const double intensity_re_hearingThreshold = intensity_in_Pa2 / hearingThreshold_in_Pa2;
119 const double intensity_in_dB_re_hearingThreshold = ( intensity_re_hearingThreshold < 1.0e-30 ? -300.0 :
120 10.0 * log10 (intensity_re_hearingThreshold) );
121 thy z [1] [iframe] = intensity_in_dB_re_hearingThreshold;
122 }
123 return thee;
124 } catch (MelderError) {
125 Melder_throw (me, U": intensity analysis not performed.");
126 }
127 }
128
Sound_to_Intensity(Sound me,double minimumPitch,double timeStep,bool subtractMeanPressure)129 autoIntensity Sound_to_Intensity (Sound me, double minimumPitch, double timeStep, bool subtractMeanPressure) {
130 const bool veryAccurate = false;
131 if (veryAccurate) {
132 autoSound up = Sound_upsample (me); // because squaring doubles the frequency content, i.e. you get super-Nyquist components
133 return Sound_to_Intensity_ (up.get(), minimumPitch, timeStep, subtractMeanPressure);
134 } else {
135 return Sound_to_Intensity_ (me, minimumPitch, timeStep, subtractMeanPressure);
136 }
137 }
138
Sound_to_IntensityTier(Sound me,double minimumPitch,double timeStep,bool subtractMean)139 autoIntensityTier Sound_to_IntensityTier (Sound me, double minimumPitch, double timeStep, bool subtractMean) {
140 try {
141 autoIntensity intensity = Sound_to_Intensity (me, minimumPitch, timeStep, subtractMean);
142 return Intensity_downto_IntensityTier (intensity.get());
143 } catch (MelderError) {
144 Melder_throw (me, U": no IntensityTier created.");
145 }
146 }
147
148 /* End of file Sound_to_Intensity.cpp */
149