Line data Source code
1 : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : /* vim:set ts=4 sw=4 sts=4 et cindent: */
3 : /*
4 : * Copyright (C) 2010 Google Inc. All rights reserved.
5 : *
6 : * Redistribution and use in source and binary forms, with or without
7 : * modification, are permitted provided that the following conditions
8 : * are met:
9 : *
10 : * 1. Redistributions of source code must retain the above copyright
11 : * notice, this list of conditions and the following disclaimer.
12 : * 2. Redistributions in binary form must reproduce the above copyright
13 : * notice, this list of conditions and the following disclaimer in the
14 : * documentation and/or other materials provided with the distribution.
15 : * 3. Neither the name of Apple Computer, Inc. ("Apple") nor the names of
16 : * its contributors may be used to endorse or promote products derived
17 : * from this software without specific prior written permission.
18 : *
19 : * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
20 : * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21 : * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22 : * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
23 : * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24 : * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25 : * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
26 : * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27 : * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
28 : * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 : */
30 :
31 : #include "FFTBlock.h"
32 :
33 : #include <complex>
34 :
35 : namespace mozilla {
36 :
37 : typedef std::complex<double> Complex;
38 :
39 0 : FFTBlock* FFTBlock::CreateInterpolatedBlock(const FFTBlock& block0, const FFTBlock& block1, double interp)
40 : {
41 0 : FFTBlock* newBlock = new FFTBlock(block0.FFTSize());
42 :
43 0 : newBlock->InterpolateFrequencyComponents(block0, block1, interp);
44 :
45 : // In the time-domain, the 2nd half of the response must be zero, to avoid circular convolution aliasing...
46 0 : int fftSize = newBlock->FFTSize();
47 0 : AlignedTArray<float> buffer(fftSize);
48 0 : newBlock->GetInverseWithoutScaling(buffer.Elements());
49 0 : AudioBufferInPlaceScale(buffer.Elements(), 1.0f / fftSize, fftSize / 2);
50 0 : PodZero(buffer.Elements() + fftSize / 2, fftSize / 2);
51 :
52 : // Put back into frequency domain.
53 0 : newBlock->PerformFFT(buffer.Elements());
54 :
55 0 : return newBlock;
56 : }
57 :
58 0 : void FFTBlock::InterpolateFrequencyComponents(const FFTBlock& block0, const FFTBlock& block1, double interp)
59 : {
60 : // FIXME : with some work, this method could be optimized
61 :
62 0 : ComplexU* dft = mOutputBuffer.Elements();
63 :
64 0 : const ComplexU* dft1 = block0.mOutputBuffer.Elements();
65 0 : const ComplexU* dft2 = block1.mOutputBuffer.Elements();
66 :
67 0 : MOZ_ASSERT(mFFTSize == block0.FFTSize());
68 0 : MOZ_ASSERT(mFFTSize == block1.FFTSize());
69 0 : double s1base = (1.0 - interp);
70 0 : double s2base = interp;
71 :
72 0 : double phaseAccum = 0.0;
73 0 : double lastPhase1 = 0.0;
74 0 : double lastPhase2 = 0.0;
75 :
76 0 : int n = mFFTSize / 2;
77 :
78 0 : dft[0].r = static_cast<float>(s1base * dft1[0].r + s2base * dft2[0].r);
79 0 : dft[n].r = static_cast<float>(s1base * dft1[n].r + s2base * dft2[n].r);
80 :
81 0 : for (int i = 1; i < n; ++i) {
82 0 : Complex c1(dft1[i].r, dft1[i].i);
83 0 : Complex c2(dft2[i].r, dft2[i].i);
84 :
85 0 : double mag1 = abs(c1);
86 0 : double mag2 = abs(c2);
87 :
88 : // Interpolate magnitudes in decibels
89 0 : double mag1db = 20.0 * log10(mag1);
90 0 : double mag2db = 20.0 * log10(mag2);
91 :
92 0 : double s1 = s1base;
93 0 : double s2 = s2base;
94 :
95 0 : double magdbdiff = mag1db - mag2db;
96 :
97 : // Empirical tweak to retain higher-frequency zeroes
98 0 : double threshold = (i > 16) ? 5.0 : 2.0;
99 :
100 0 : if (magdbdiff < -threshold && mag1db < 0.0) {
101 0 : s1 = pow(s1, 0.75);
102 0 : s2 = 1.0 - s1;
103 0 : } else if (magdbdiff > threshold && mag2db < 0.0) {
104 0 : s2 = pow(s2, 0.75);
105 0 : s1 = 1.0 - s2;
106 : }
107 :
108 : // Average magnitude by decibels instead of linearly
109 0 : double magdb = s1 * mag1db + s2 * mag2db;
110 0 : double mag = pow(10.0, 0.05 * magdb);
111 :
112 : // Now, deal with phase
113 0 : double phase1 = arg(c1);
114 0 : double phase2 = arg(c2);
115 :
116 0 : double deltaPhase1 = phase1 - lastPhase1;
117 0 : double deltaPhase2 = phase2 - lastPhase2;
118 0 : lastPhase1 = phase1;
119 0 : lastPhase2 = phase2;
120 :
121 : // Unwrap phase deltas
122 0 : if (deltaPhase1 > M_PI)
123 0 : deltaPhase1 -= 2.0 * M_PI;
124 0 : if (deltaPhase1 < -M_PI)
125 0 : deltaPhase1 += 2.0 * M_PI;
126 0 : if (deltaPhase2 > M_PI)
127 0 : deltaPhase2 -= 2.0 * M_PI;
128 0 : if (deltaPhase2 < -M_PI)
129 0 : deltaPhase2 += 2.0 * M_PI;
130 :
131 : // Blend group-delays
132 : double deltaPhaseBlend;
133 :
134 0 : if (deltaPhase1 - deltaPhase2 > M_PI)
135 0 : deltaPhaseBlend = s1 * deltaPhase1 + s2 * (2.0 * M_PI + deltaPhase2);
136 0 : else if (deltaPhase2 - deltaPhase1 > M_PI)
137 0 : deltaPhaseBlend = s1 * (2.0 * M_PI + deltaPhase1) + s2 * deltaPhase2;
138 : else
139 0 : deltaPhaseBlend = s1 * deltaPhase1 + s2 * deltaPhase2;
140 :
141 0 : phaseAccum += deltaPhaseBlend;
142 :
143 : // Unwrap
144 0 : if (phaseAccum > M_PI)
145 0 : phaseAccum -= 2.0 * M_PI;
146 0 : if (phaseAccum < -M_PI)
147 0 : phaseAccum += 2.0 * M_PI;
148 :
149 0 : dft[i].r = static_cast<float>(mag * cos(phaseAccum));
150 0 : dft[i].i = static_cast<float>(mag * sin(phaseAccum));
151 : }
152 0 : }
153 :
154 0 : double FFTBlock::ExtractAverageGroupDelay()
155 : {
156 0 : ComplexU* dft = mOutputBuffer.Elements();
157 :
158 0 : double aveSum = 0.0;
159 0 : double weightSum = 0.0;
160 0 : double lastPhase = 0.0;
161 :
162 0 : int halfSize = FFTSize() / 2;
163 :
164 0 : const double kSamplePhaseDelay = (2.0 * M_PI) / double(FFTSize());
165 :
166 : // Remove DC offset
167 0 : dft[0].r = 0.0f;
168 :
169 : // Calculate weighted average group delay
170 0 : for (int i = 1; i < halfSize; i++) {
171 0 : Complex c(dft[i].r, dft[i].i);
172 0 : double mag = abs(c);
173 0 : double phase = arg(c);
174 :
175 0 : double deltaPhase = phase - lastPhase;
176 0 : lastPhase = phase;
177 :
178 : // Unwrap
179 0 : if (deltaPhase < -M_PI)
180 0 : deltaPhase += 2.0 * M_PI;
181 0 : if (deltaPhase > M_PI)
182 0 : deltaPhase -= 2.0 * M_PI;
183 :
184 0 : aveSum += mag * deltaPhase;
185 0 : weightSum += mag;
186 : }
187 :
188 : // Note how we invert the phase delta wrt frequency since this is how group delay is defined
189 0 : double ave = aveSum / weightSum;
190 0 : double aveSampleDelay = -ave / kSamplePhaseDelay;
191 :
192 : // Leave 20 sample headroom (for leading edge of impulse)
193 0 : aveSampleDelay -= 20.0;
194 0 : if (aveSampleDelay <= 0.0)
195 0 : return 0.0;
196 :
197 : // Remove average group delay (minus 20 samples for headroom)
198 0 : AddConstantGroupDelay(-aveSampleDelay);
199 :
200 0 : return aveSampleDelay;
201 : }
202 :
203 0 : void FFTBlock::AddConstantGroupDelay(double sampleFrameDelay)
204 : {
205 0 : int halfSize = FFTSize() / 2;
206 :
207 0 : ComplexU* dft = mOutputBuffer.Elements();
208 :
209 0 : const double kSamplePhaseDelay = (2.0 * M_PI) / double(FFTSize());
210 :
211 0 : double phaseAdj = -sampleFrameDelay * kSamplePhaseDelay;
212 :
213 : // Add constant group delay
214 0 : for (int i = 1; i < halfSize; i++) {
215 0 : Complex c(dft[i].r, dft[i].i);
216 0 : double mag = abs(c);
217 0 : double phase = arg(c);
218 :
219 0 : phase += i * phaseAdj;
220 :
221 0 : dft[i].r = static_cast<float>(mag * cos(phase));
222 0 : dft[i].i = static_cast<float>(mag * sin(phase));
223 : }
224 0 : }
225 :
226 : } // namespace mozilla
|
