1/*
2 * Copyright (C) 2010 Google Inc. All rights reserved.
3 *
4 * Redistribution and use in source and binary forms, with or without
5 * modification, are permitted provided that the following conditions
6 * are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 * 3. Neither the name of Apple Inc. ("Apple") nor the names of
14 * its contributors may be used to endorse or promote products derived
15 * from this software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
18 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
19 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
20 * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
21 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29#include "config.h"
30
31#if ENABLE(WEB_AUDIO)
32
33#include "FFTFrame.h"
34
35#include "Logging.h"
36#include <complex>
37#include <wtf/MathExtras.h>
38
39#ifndef NDEBUG
40#include <stdio.h>
41#endif
42
43namespace WebCore {
44
45void FFTFrame::doPaddedFFT(const float* data, size_t dataSize)
46{
47 // Zero-pad the impulse response
48 AudioFloatArray paddedResponse(fftSize()); // zero-initialized
49 paddedResponse.copyToRange(data, 0, dataSize);
50
51 // Get the frequency-domain version of padded response
52 doFFT(paddedResponse.data());
53}
54
55std::unique_ptr<FFTFrame> FFTFrame::createInterpolatedFrame(const FFTFrame& frame1, const FFTFrame& frame2, double x)
56{
57 auto newFrame = std::make_unique<FFTFrame>(frame1.fftSize());
58
59 newFrame->interpolateFrequencyComponents(frame1, frame2, x);
60
61 // In the time-domain, the 2nd half of the response must be zero, to avoid circular convolution aliasing...
62 int fftSize = newFrame->fftSize();
63 AudioFloatArray buffer(fftSize);
64 newFrame->doInverseFFT(buffer.data());
65 buffer.zeroRange(fftSize / 2, fftSize);
66
67 // Put back into frequency domain.
68 newFrame->doFFT(buffer.data());
69
70 return newFrame;
71}
72
73void FFTFrame::interpolateFrequencyComponents(const FFTFrame& frame1, const FFTFrame& frame2, double interp)
74{
75 // FIXME : with some work, this method could be optimized
76
77 float* realP = realData();
78 float* imagP = imagData();
79
80 const float* realP1 = frame1.realData();
81 const float* imagP1 = frame1.imagData();
82 const float* realP2 = frame2.realData();
83 const float* imagP2 = frame2.imagData();
84
85 m_FFTSize = frame1.fftSize();
86 m_log2FFTSize = frame1.log2FFTSize();
87
88 double s1base = (1.0 - interp);
89 double s2base = interp;
90
91 double phaseAccum = 0.0;
92 double lastPhase1 = 0.0;
93 double lastPhase2 = 0.0;
94
95 realP[0] = static_cast<float>(s1base * realP1[0] + s2base * realP2[0]);
96 imagP[0] = static_cast<float>(s1base * imagP1[0] + s2base * imagP2[0]);
97
98 int n = m_FFTSize / 2;
99
100 for (int i = 1; i < n; ++i) {
101 std::complex<double> c1(realP1[i], imagP1[i]);
102 std::complex<double> c2(realP2[i], imagP2[i]);
103
104 double mag1 = abs(c1);
105 double mag2 = abs(c2);
106
107 // Interpolate magnitudes in decibels
108 double mag1db = 20.0 * log10(mag1);
109 double mag2db = 20.0 * log10(mag2);
110
111 double s1 = s1base;
112 double s2 = s2base;
113
114 double magdbdiff = mag1db - mag2db;
115
116 // Empirical tweak to retain higher-frequency zeroes
117 double threshold = (i > 16) ? 5.0 : 2.0;
118
119 if (magdbdiff < -threshold && mag1db < 0.0) {
120 s1 = pow(s1, 0.75);
121 s2 = 1.0 - s1;
122 } else if (magdbdiff > threshold && mag2db < 0.0) {
123 s2 = pow(s2, 0.75);
124 s1 = 1.0 - s2;
125 }
126
127 // Average magnitude by decibels instead of linearly
128 double magdb = s1 * mag1db + s2 * mag2db;
129 double mag = pow(10.0, 0.05 * magdb);
130
131 // Now, deal with phase
132 double phase1 = arg(c1);
133 double phase2 = arg(c2);
134
135 double deltaPhase1 = phase1 - lastPhase1;
136 double deltaPhase2 = phase2 - lastPhase2;
137 lastPhase1 = phase1;
138 lastPhase2 = phase2;
139
140 // Unwrap phase deltas
141 if (deltaPhase1 > piDouble)
142 deltaPhase1 -= 2.0 * piDouble;
143 if (deltaPhase1 < -piDouble)
144 deltaPhase1 += 2.0 * piDouble;
145 if (deltaPhase2 > piDouble)
146 deltaPhase2 -= 2.0 * piDouble;
147 if (deltaPhase2 < -piDouble)
148 deltaPhase2 += 2.0 * piDouble;
149
150 // Blend group-delays
151 double deltaPhaseBlend;
152
153 if (deltaPhase1 - deltaPhase2 > piDouble)
154 deltaPhaseBlend = s1 * deltaPhase1 + s2 * (2.0 * piDouble + deltaPhase2);
155 else if (deltaPhase2 - deltaPhase1 > piDouble)
156 deltaPhaseBlend = s1 * (2.0 * piDouble + deltaPhase1) + s2 * deltaPhase2;
157 else
158 deltaPhaseBlend = s1 * deltaPhase1 + s2 * deltaPhase2;
159
160 phaseAccum += deltaPhaseBlend;
161
162 // Unwrap
163 if (phaseAccum > piDouble)
164 phaseAccum -= 2.0 * piDouble;
165 if (phaseAccum < -piDouble)
166 phaseAccum += 2.0 * piDouble;
167
168 std::complex<double> c = std::polar(mag, phaseAccum);
169
170 realP[i] = static_cast<float>(c.real());
171 imagP[i] = static_cast<float>(c.imag());
172 }
173}
174
175double FFTFrame::extractAverageGroupDelay()
176{
177 float* realP = realData();
178 float* imagP = imagData();
179
180 double aveSum = 0.0;
181 double weightSum = 0.0;
182 double lastPhase = 0.0;
183
184 int halfSize = fftSize() / 2;
185
186 const double kSamplePhaseDelay = (2.0 * piDouble) / double(fftSize());
187
188 // Calculate weighted average group delay
189 for (int i = 0; i < halfSize; i++) {
190 std::complex<double> c(realP[i], imagP[i]);
191 double mag = abs(c);
192 double phase = arg(c);
193
194 double deltaPhase = phase - lastPhase;
195 lastPhase = phase;
196
197 // Unwrap
198 if (deltaPhase < -piDouble)
199 deltaPhase += 2.0 * piDouble;
200 if (deltaPhase > piDouble)
201 deltaPhase -= 2.0 * piDouble;
202
203 aveSum += mag * deltaPhase;
204 weightSum += mag;
205 }
206
207 // Note how we invert the phase delta wrt frequency since this is how group delay is defined
208 double ave = aveSum / weightSum;
209 double aveSampleDelay = -ave / kSamplePhaseDelay;
210
211 // Leave 20 sample headroom (for leading edge of impulse)
212 if (aveSampleDelay > 20.0)
213 aveSampleDelay -= 20.0;
214
215 // Remove average group delay (minus 20 samples for headroom)
216 addConstantGroupDelay(-aveSampleDelay);
217
218 // Remove DC offset
219 realP[0] = 0.0f;
220
221 return aveSampleDelay;
222}
223
224void FFTFrame::addConstantGroupDelay(double sampleFrameDelay)
225{
226 int halfSize = fftSize() / 2;
227
228 float* realP = realData();
229 float* imagP = imagData();
230
231 const double kSamplePhaseDelay = (2.0 * piDouble) / double(fftSize());
232
233 double phaseAdj = -sampleFrameDelay * kSamplePhaseDelay;
234
235 // Add constant group delay
236 for (int i = 1; i < halfSize; i++) {
237 std::complex<double> c(realP[i], imagP[i]);
238 double mag = abs(c);
239 double phase = arg(c);
240
241 phase += i * phaseAdj;
242
243 std::complex<double> c2 = std::polar(mag, phase);
244
245 realP[i] = static_cast<float>(c2.real());
246 imagP[i] = static_cast<float>(c2.imag());
247 }
248}
249
250#ifndef NDEBUG
251void FFTFrame::print()
252{
253 FFTFrame& frame = *this;
254 float* realP = frame.realData();
255 float* imagP = frame.imagData();
256 LOG(WebAudio, "**** \n");
257 LOG(WebAudio, "DC = %f : nyquist = %f\n", realP[0], imagP[0]);
258
259 int n = m_FFTSize / 2;
260
261 for (int i = 1; i < n; i++) {
262 double mag = sqrt(realP[i] * realP[i] + imagP[i] * imagP[i]);
263 double phase = atan2(realP[i], imagP[i]);
264
265 LOG(WebAudio, "[%d] (%f %f)\n", i, mag, phase);
266 }
267 LOG(WebAudio, "****\n");
268}
269#endif // NDEBUG
270
271} // namespace WebCore
272
273#endif // ENABLE(WEB_AUDIO)
274