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