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