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 Computer, 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#include <wtf/OwnPtr.h>
39
40#ifndef NDEBUG
41#include <stdio.h>
42#endif
43
44namespace WebCore {
45
46void FFTFrame::doPaddedFFT(const float* data, size_t dataSize)
47{
48    // Zero-pad the impulse response
49    AudioFloatArray paddedResponse(fftSize()); // zero-initialized
50    paddedResponse.copyToRange(data, 0, dataSize);
51
52    // Get the frequency-domain version of padded response
53    doFFT(paddedResponse.data());
54}
55
56PassOwnPtr<FFTFrame> FFTFrame::createInterpolatedFrame(const FFTFrame& frame1, const FFTFrame& frame2, double x)
57{
58    OwnPtr<FFTFrame> newFrame = adoptPtr(new FFTFrame(frame1.fftSize()));
59
60    newFrame->interpolateFrequencyComponents(frame1, frame2, x);
61
62    // In the time-domain, the 2nd half of the response must be zero, to avoid circular convolution aliasing...
63    int fftSize = newFrame->fftSize();
64    AudioFloatArray buffer(fftSize);
65    newFrame->doInverseFFT(buffer.data());
66    buffer.zeroRange(fftSize / 2, fftSize);
67
68    // Put back into frequency domain.
69    newFrame->doFFT(buffer.data());
70
71    return newFrame.release();
72}
73
74void FFTFrame::interpolateFrequencyComponents(const FFTFrame& frame1, const FFTFrame& frame2, double interp)
75{
76    // FIXME : with some work, this method could be optimized
77
78    float* realP = realData();
79    float* imagP = imagData();
80
81    const float* realP1 = frame1.realData();
82    const float* imagP1 = frame1.imagData();
83    const float* realP2 = frame2.realData();
84    const float* imagP2 = frame2.imagData();
85
86    m_FFTSize = frame1.fftSize();
87    m_log2FFTSize = frame1.log2FFTSize();
88
89    double s1base = (1.0 - interp);
90    double s2base = interp;
91
92    double phaseAccum = 0.0;
93    double lastPhase1 = 0.0;
94    double lastPhase2 = 0.0;
95
96    realP[0] = static_cast<float>(s1base * realP1[0] + s2base * realP2[0]);
97    imagP[0] = static_cast<float>(s1base * imagP1[0] + s2base * imagP2[0]);
98
99    int n = m_FFTSize / 2;
100
101    for (int i = 1; i < n; ++i) {
102        std::complex<double> c1(realP1[i], imagP1[i]);
103        std::complex<double> c2(realP2[i], imagP2[i]);
104
105        double mag1 = abs(c1);
106        double mag2 = abs(c2);
107
108        // Interpolate magnitudes in decibels
109        double mag1db = 20.0 * log10(mag1);
110        double mag2db = 20.0 * log10(mag2);
111
112        double s1 = s1base;
113        double s2 = s2base;
114
115        double magdbdiff = mag1db - mag2db;
116
117        // Empirical tweak to retain higher-frequency zeroes
118        double threshold =  (i > 16) ? 5.0 : 2.0;
119
120        if (magdbdiff < -threshold && mag1db < 0.0) {
121            s1 = pow(s1, 0.75);
122            s2 = 1.0 - s1;
123        } else if (magdbdiff > threshold && mag2db < 0.0) {
124            s2 = pow(s2, 0.75);
125            s1 = 1.0 - s2;
126        }
127
128        // Average magnitude by decibels instead of linearly
129        double magdb = s1 * mag1db + s2 * mag2db;
130        double mag = pow(10.0, 0.05 * magdb);
131
132        // Now, deal with phase
133        double phase1 = arg(c1);
134        double phase2 = arg(c2);
135
136        double deltaPhase1 = phase1 - lastPhase1;
137        double deltaPhase2 = phase2 - lastPhase2;
138        lastPhase1 = phase1;
139        lastPhase2 = phase2;
140
141        // Unwrap phase deltas
142        if (deltaPhase1 > piDouble)
143            deltaPhase1 -= 2.0 * piDouble;
144        if (deltaPhase1 < -piDouble)
145            deltaPhase1 += 2.0 * piDouble;
146        if (deltaPhase2 > piDouble)
147            deltaPhase2 -= 2.0 * piDouble;
148        if (deltaPhase2 < -piDouble)
149            deltaPhase2 += 2.0 * piDouble;
150
151        // Blend group-delays
152        double deltaPhaseBlend;
153
154        if (deltaPhase1 - deltaPhase2 > piDouble)
155            deltaPhaseBlend = s1 * deltaPhase1 + s2 * (2.0 * piDouble + deltaPhase2);
156        else if (deltaPhase2 - deltaPhase1 > piDouble)
157            deltaPhaseBlend = s1 * (2.0 * piDouble + deltaPhase1) + s2 * deltaPhase2;
158        else
159            deltaPhaseBlend = s1 * deltaPhase1 + s2 * deltaPhase2;
160
161        phaseAccum += deltaPhaseBlend;
162
163        // Unwrap
164        if (phaseAccum > piDouble)
165            phaseAccum -= 2.0 * piDouble;
166        if (phaseAccum < -piDouble)
167            phaseAccum += 2.0 * piDouble;
168
169        std::complex<double> c = std::polar(mag, phaseAccum);
170
171        realP[i] = static_cast<float>(c.real());
172        imagP[i] = static_cast<float>(c.imag());
173    }
174}
175
176double FFTFrame::extractAverageGroupDelay()
177{
178    float* realP = realData();
179    float* imagP = imagData();
180
181    double aveSum = 0.0;
182    double weightSum = 0.0;
183    double lastPhase = 0.0;
184
185    int halfSize = fftSize() / 2;
186
187    const double kSamplePhaseDelay = (2.0 * piDouble) / double(fftSize());
188
189    // Calculate weighted average group delay
190    for (int i = 0; i < halfSize; i++) {
191        std::complex<double> c(realP[i], imagP[i]);
192        double mag = abs(c);
193        double phase = arg(c);
194
195        double deltaPhase = phase - lastPhase;
196        lastPhase = phase;
197
198        // Unwrap
199        if (deltaPhase < -piDouble)
200            deltaPhase += 2.0 * piDouble;
201        if (deltaPhase > piDouble)
202            deltaPhase -= 2.0 * piDouble;
203
204        aveSum += mag * deltaPhase;
205        weightSum += mag;
206    }
207
208    // Note how we invert the phase delta wrt frequency since this is how group delay is defined
209    double ave = aveSum / weightSum;
210    double aveSampleDelay = -ave / kSamplePhaseDelay;
211
212    // Leave 20 sample headroom (for leading edge of impulse)
213    if (aveSampleDelay > 20.0)
214        aveSampleDelay -= 20.0;
215
216    // Remove average group delay (minus 20 samples for headroom)
217    addConstantGroupDelay(-aveSampleDelay);
218
219    // Remove DC offset
220    realP[0] = 0.0f;
221
222    return aveSampleDelay;
223}
224
225void FFTFrame::addConstantGroupDelay(double sampleFrameDelay)
226{
227    int halfSize = fftSize() / 2;
228
229    float* realP = realData();
230    float* imagP = imagData();
231
232    const double kSamplePhaseDelay = (2.0 * piDouble) / double(fftSize());
233
234    double phaseAdj = -sampleFrameDelay * kSamplePhaseDelay;
235
236    // Add constant group delay
237    for (int i = 1; i < halfSize; i++) {
238        std::complex<double> c(realP[i], imagP[i]);
239        double mag = abs(c);
240        double phase = arg(c);
241
242        phase += i * phaseAdj;
243
244        std::complex<double> c2 = std::polar(mag, phase);
245
246        realP[i] = static_cast<float>(c2.real());
247        imagP[i] = static_cast<float>(c2.imag());
248    }
249}
250
251#ifndef NDEBUG
252void FFTFrame::print()
253{
254    FFTFrame& frame = *this;
255    float* realP = frame.realData();
256    float* imagP = frame.imagData();
257    LOG(WebAudio, "**** \n");
258    LOG(WebAudio, "DC = %f : nyquist = %f\n", realP[0], imagP[0]);
259
260    int n = m_FFTSize / 2;
261
262    for (int i = 1; i < n; i++) {
263        double mag = sqrt(realP[i] * realP[i] + imagP[i] * imagP[i]);
264        double phase = atan2(realP[i], imagP[i]);
265
266        LOG(WebAudio, "[%d] (%f %f)\n", i, mag, phase);
267    }
268    LOG(WebAudio, "****\n");
269}
270#endif // NDEBUG
271
272} // namespace WebCore
273
274#endif // ENABLE(WEB_AUDIO)
275