Release 1.7
This commit is contained in:
parent
c934ef2c56
commit
80ef6d4c87
SSH key fingerprint:
SHA256:OLTRf6Fl87G52SiR7sWLGNzlJt4WOX+tfI2yxo0z7xk
13
CHANGES
13
CHANGES
|
|
@ -1,4 +1,17 @@
|
|||
|
||||
1.7
|
||||
|
||||
Miscellaneous bug fixes.
|
||||
|
||||
Support lower numbers of bands per octave, down to 4.
|
||||
|
||||
Further improve the performance of analyzing short signal blocks.
|
||||
|
||||
The "Frequency-Domain Filtering" and "Streaming" examples now use
|
||||
a white noise and impulse signal, respectively.
|
||||
|
||||
1.6
|
||||
|
||||
Add "API Introduction" documentation section that was missing
|
||||
from version 1.5, causing broken links.
|
||||
|
||||
|
|
|
|||
Binary file not shown.
|
Before Width: | Height: | Size: 82 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 73 KiB |
|
|
@ -111,7 +111,7 @@ sampling frequency. The gain is normalized to unity at 20 Hz.
|
|||
</p>
|
||||
<pre>
|
||||
for (int band = analyzer.bandpass_bands_begin(); band < analyzer.bandpass_bands_end(); band++) {
|
||||
float f_hz = analyzer.band_ff(band) * fs;
|
||||
double f_hz = analyzer.band_ff(band) * fs;
|
||||
band_gains[band] = 1.0 / sqrt(f_hz / 20.0);
|
||||
}
|
||||
</pre>
|
||||
|
|
@ -233,7 +233,7 @@ complete program:
|
|||
<p>Like <a href="render.html#compiling">Example 1</a>, this example
|
||||
can be built using a one-line build command:
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD">
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` filter.cc `pkg-config --libs sndfile` -o filter
|
||||
</pre>
|
||||
<p>Or using the vDSP FFT on macOS:</p>
|
||||
|
|
@ -241,19 +241,18 @@ c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` filter.cc `pkg
|
|||
c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` filter.cc `pkg-config --libs sndfile` -framework Accelerate -o filter
|
||||
</pre>
|
||||
<p>Or using PFFFT (see <a href="render.html#compiling">Example 1</a> for how to download and build PFFFT):</p>
|
||||
<pre class="build Linux NetBSD">
|
||||
<pre class="build">
|
||||
c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` filter.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o filter
|
||||
</pre>
|
||||
|
||||
<h2>Running</h2>
|
||||
<p>To filter the file <code>guitar.wav</code> that was downloaded in
|
||||
Example 1, simply run</p>
|
||||
<p>Running the following shell commands will download an example
|
||||
audio file containing five seconds of white noise and filter it,
|
||||
producing pink noise.</p>
|
||||
<pre class="run">
|
||||
./filter guitar.wav guitar_filtered.wav
|
||||
wget http://download.gaborator.com/audio/white_noise.wav
|
||||
./filter white_noise.wav pink_noise.wav
|
||||
</pre>
|
||||
<p>The resulting lowpass filtered audio in <code>guitar_filtered.wav</code> will
|
||||
sound muffled compared to the original, but less so than it would with a
|
||||
6 dB/octave filter.</p>
|
||||
|
||||
<h2>Frequency response</h2>
|
||||
<p>The following plot shows the actual measured frequency response of the
|
||||
|
|
|
|||
BIN
doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png
Normal file
BIN
doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 122 KiB |
BIN
doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png
Normal file
BIN
doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 107 KiB |
BIN
doc/gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png
Normal file
BIN
doc/gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 14 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 9.6 KiB |
|
|
@ -1,20 +1,21 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-2020 Andreas Gustafsson. This file is part of
|
||||
Copyright (C) 2018-2021 Andreas Gustafsson. This file is part of
|
||||
the Gaborator library source distribution. See the file LICENSE at
|
||||
the top level of the distribution for license information.
|
||||
-->
|
||||
<html>
|
||||
<head>
|
||||
<link rel="stylesheet" href="doc.css" />
|
||||
<link rel="icon" type="image/png" href="favicon64.png" />
|
||||
<title>The Gaborator</title>
|
||||
</head>
|
||||
<body>
|
||||
<h1>The Gaborator</h1>
|
||||
<p>The Gaborator is a library that generates constant-Q spectrograms
|
||||
for visualization and analysis of audio signals. It also supports an
|
||||
accurate inverse transformation of the spectrogram coefficients back into
|
||||
audio for spectral effects and editing.</p>
|
||||
for visualization and analysis of audio signals. It also supports a
|
||||
fast and accurate inverse transformation of the spectrogram coefficients
|
||||
back into audio for spectral effects and editing.</p>
|
||||
|
||||
<p>The Gaborator implements the invertible constant-Q transform of
|
||||
Velasco, Holighaus, Dörfler, and Grill, described in the papers
|
||||
|
|
@ -26,8 +27,8 @@ using Gaussian bandpass filters and an efficient multi-rate architecture.
|
|||
</p>
|
||||
|
||||
<p>The Gaborator is written in C++11 and compatible with C++14 and C++17.
|
||||
It has been tested on macOS, Linux, NetBSD, and iOS, on Intel x86_64 and ARM
|
||||
processors.</p>
|
||||
It has been tested on macOS, Linux, NetBSD, FreeBSD, and iOS, on Intel
|
||||
x86_64 and ARM processors.</p>
|
||||
|
||||
<p>The Gaborator is open source under the GNU Affero General Public
|
||||
License, version 3, and is also available for commercial licensing.
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-2020 Andreas Gustafsson. This file is part of
|
||||
Copyright (C) 2018-2021 Andreas Gustafsson. This file is part of
|
||||
the Gaborator library source distribution. See the file LICENSE at
|
||||
the top level of the distribution for license information.
|
||||
-->
|
||||
|
|
@ -67,7 +67,7 @@ typical f<sub>min</sub> for audio work would be 0.00045, but
|
|||
that would make the plot hard to read because both the lowpass filter
|
||||
and the lowest-frequency bandpass filters would be extremely narrow.</p>
|
||||
|
||||
<img src="allkernels_bpo12_ffmin0.03_ffref0.5_anl_wob.png" alt="Analysis filters">
|
||||
<img src="gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png" alt="Analysis filters">
|
||||
|
||||
<p>The output of each bandpass filter is shifted down in frequency to
|
||||
a complex quadrature baseband. The baseband signal is then resampled
|
||||
|
|
@ -95,7 +95,7 @@ here. The plot covers a time range of 128 signal samples, but
|
|||
conceptually, the grid extends arbitrarily far in time, in both the
|
||||
positive and the negative direction.</p>
|
||||
|
||||
<img src="grid_bpo12_ffmin0.03_ffref0.5_wob.png" alt="Sampling grid">
|
||||
<img src="gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png" alt="Sampling grid">
|
||||
|
||||
<h2>Resynthesis</h2>
|
||||
|
||||
|
|
@ -107,7 +107,7 @@ that is a <i>dual</i> of the analysis filter bank. The following
|
|||
plot shows the frequency responses of the reconstruction filters
|
||||
corresponding to the analysis filters shown earlier.</p>
|
||||
|
||||
<img src="allkernels_bpo12_ffmin0.03_ffref0.5_syn_wob.png" alt="Reconstruction filters">
|
||||
<img src="gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png" alt="Reconstruction filters">
|
||||
|
||||
<p>Although the bandpass filters may look similar to the Gaussian
|
||||
filters of the analysis filter bank, their shapes are actually subtly
|
||||
|
|
|
|||
|
|
@ -71,15 +71,17 @@ applying the same effect to an entire one-minute recording runs in a
|
|||
fraction of a second.</p>
|
||||
|
||||
<p>In analysis and visualization applications that don't need to
|
||||
perform resynthesis, it may be possible to partly hide the latency by
|
||||
perform resynthesis, it is possible to partly hide the latency by
|
||||
taking advantage of the fact that the coefficients for the higher
|
||||
frequencies exhibit lower latency than those for low frequencies.
|
||||
For example, a live spectrogram display could update the
|
||||
high-frequency parts of the display before the corresponding
|
||||
low-frequency parts. Alternatively, low-frequency parts of the
|
||||
spectrogram could be drawn multiple times, effectively animating
|
||||
spectrogram may be drawn multiple times, effectively animating
|
||||
the display of the low-frequency coefficients as they converge to
|
||||
their final values.</p>
|
||||
their final values. This approach can be seen in action in
|
||||
the <a href="https://waxingwave.com/spectrolite/">Spectrolite</a>
|
||||
iOS app.</p>
|
||||
|
||||
<h2>Does it support small blocks sizes?</h2>
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
<!DOCTYPE html>
|
||||
<!--
|
||||
Copyright (C) 2018-2020 Andreas Gustafsson. This file is part of
|
||||
Copyright (C) 2018-2021 Andreas Gustafsson. This file is part of
|
||||
the Gaborator library source distribution. See the file LICENSE at
|
||||
the top level of the distribution for license information.
|
||||
-->
|
||||
|
|
@ -32,8 +32,8 @@ parameters(unsigned int bands_per_octave,
|
|||
<dl>
|
||||
<dt><code>bands_per_octave</code></dt>
|
||||
<dd>The number of frequency bands per octave.
|
||||
Values from 6 to 384 (inclusive) are supported. Values outside
|
||||
this range may not work, or may cause degraded performance.</dd>
|
||||
Values from 4 to 384 (inclusive) are supported.
|
||||
</dd>
|
||||
<dt><code>ff_min</code></dt>
|
||||
<dd>The lower limit of the analysis frequency range, in units of the
|
||||
sample rate. The analysis filter bank will extend low enough in
|
||||
|
|
|
|||
|
|
@ -331,7 +331,7 @@ subdirectory.
|
|||
You need to have <i>libsndfile</i> is installed and supported by
|
||||
<code>pkg-config</code>.
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD">
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` render.cc `pkg-config --libs sndfile` -o render
|
||||
</pre>
|
||||
|
||||
|
|
@ -353,24 +353,24 @@ You can get the latest version from
|
|||
or the exact version that was used for testing from gaborator.com:
|
||||
</p>
|
||||
<!-- ftp https://bitbucket.org/jpommier/pffft/get/29e4f76ac53b.zip -->
|
||||
<pre class="build Linux NetBSD">
|
||||
<pre class="build Linux NetBSD FreeBSD">
|
||||
wget http://download.gaborator.com/mirror/pffft/29e4f76ac53b.zip
|
||||
unzip 29e4f76ac53b.zip
|
||||
mv jpommier-pffft-29e4f76ac53b pffft
|
||||
</pre>
|
||||
<p>Then, compile it:</p>
|
||||
<pre class="build Linux NetBSD">
|
||||
<pre class="build Linux NetBSD FreeBSD">
|
||||
cc -c -O3 -ffast-math pffft/pffft.c -o pffft/pffft.o
|
||||
</pre>
|
||||
<p>(If you are building for ARM, you will need to add <code>-mfpu=neon</code> to
|
||||
both the above compilation command and the ones below.)</p>
|
||||
<p>PFFFT is single precision only, but it comes with a copy of FFTPACK which can
|
||||
be used for double-precision FFTs. Let's compile that, too:</p>
|
||||
<pre class="build Linux NetBSD">
|
||||
<pre class="build Linux NetBSD FreeBSD">
|
||||
cc -c -O3 -ffast-math -DFFTPACK_DOUBLE_PRECISION pffft/fftpack.c -o pffft/fftpack.o
|
||||
</pre>
|
||||
<p>Then build the example and link it with both PFFFT and FFTPACK:</p>
|
||||
<pre class="build Linux NetBSD">
|
||||
<pre class="build Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` render.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o render
|
||||
</pre>
|
||||
|
||||
|
|
|
|||
|
|
@ -96,7 +96,7 @@ and a frequency range of 3 decades (0.0005 to 0.5 times the sample rate):</p>
|
|||
<p>Like <a href="render.html#compiling">Example 1</a>, this example
|
||||
can be built using a one-line build command:
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD">
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` snr.cc `pkg-config --libs sndfile` -o snr
|
||||
</pre>
|
||||
<p>Or using the vDSP FFT on macOS:</p>
|
||||
|
|
@ -104,7 +104,7 @@ c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` snr.cc `pkg-co
|
|||
c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` snr.cc `pkg-config --libs sndfile` -framework Accelerate -o snr
|
||||
</pre>
|
||||
<p>Or using PFFFT (see <a href="render.html#compiling">Example 1</a> for how to download and build PFFFT):</p>
|
||||
<pre class="build Linux NetBSD">
|
||||
<pre class="build">
|
||||
c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` snr.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o snr
|
||||
</pre>
|
||||
|
||||
|
|
|
|||
|
|
@ -247,7 +247,7 @@ it's now OK to forget them and free the memory they used:
|
|||
<h2>Compiling</h2>
|
||||
<p>Like the previous ones, this example can also be built using a one-line build command:
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD">
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` stream.cc `pkg-config --libs sndfile` -o stream
|
||||
</pre>
|
||||
<p>Or using the vDSP FFT on macOS:</p>
|
||||
|
|
@ -255,18 +255,22 @@ c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` stream.cc `pkg
|
|||
c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` stream.cc `pkg-config --libs sndfile` -framework Accelerate -o stream
|
||||
</pre>
|
||||
<p>Or using PFFFT (see <a href="render.html#compiling">Example 1</a> for how to download and build PFFFT):</p>
|
||||
<pre class="build Linux NetBSD">
|
||||
<pre class="build">
|
||||
c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` stream.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o stream
|
||||
</pre>
|
||||
|
||||
<h2>Running</h2>
|
||||
<p>To process the file <code>guitar.wav</code> that was downloaded in
|
||||
Example 1, run</p>
|
||||
<p>Running the following shell commands will download an example
|
||||
audio file containing an impulse (a single sample of maximum amplitude)
|
||||
padded with silence to a total of 65536 samples, and process it.</p>
|
||||
<pre class="run">
|
||||
./stream guitar.wav guitar_streamed.wav
|
||||
wget http://download.gaborator.com/audio/impulse.wav
|
||||
./stream impulse.wav impulse_streamed.wav
|
||||
</pre>
|
||||
<p>The file <code>guitar_streamed.wav</code> will be identical to <code>guitar.wav</code>
|
||||
except that the signal will be of opposite polarity, and delayed by the latency of
|
||||
|
||||
<p>The file <code>impulse_streamed.wav</code> will be identical to
|
||||
<code>impulse.wav</code> except that the impulse will be of
|
||||
opposite polarity, and delayed by the latency of
|
||||
<code>analysis_support + synthesis_support</code> samples.</p>
|
||||
|
||||
<div class="nav"><span class="prev"><a href="filter.html">Previous: Example 2: Frequency-Domain Filtering</a></span><span class="next"><a href="snr.html">Next: Example 4: Measuring the Signal-to-Noise Ratio</a></span></div>
|
||||
|
|
|
|||
|
|
@ -188,7 +188,7 @@ we need to choose a file format for the output file by filling in the
|
|||
<p>Like <a href="render.html#compiling">Example 1</a>, this example
|
||||
can be built using a one-line build command:
|
||||
</p>
|
||||
<pre class="build Darwin Linux NetBSD">
|
||||
<pre class="build Darwin Linux NetBSD FreeBSD">
|
||||
c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` synth.cc `pkg-config --libs sndfile` -o synth
|
||||
</pre>
|
||||
<p>Or using the vDSP FFT on macOS:</p>
|
||||
|
|
@ -196,7 +196,7 @@ c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` synth.cc `pkg-
|
|||
c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` synth.cc `pkg-config --libs sndfile` -framework Accelerate -o synth
|
||||
</pre>
|
||||
<p>Or using PFFFT (see <a href="render.html">Example 1</a> for how to download and build PFFFT):</p>
|
||||
<pre class="build Linux NetBSD">
|
||||
<pre class="build">
|
||||
c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --cflags sndfile` synth.cc pffft/pffft.o pffft/fftpack.o `pkg-config --libs sndfile` -o synth
|
||||
</pre>
|
||||
|
||||
|
|
|
|||
|
|
@ -32,7 +32,7 @@ int main(int argc, char **argv) {
|
|||
gaborator::analyzer<float> analyzer(params);
|
||||
std::vector<float> band_gains(analyzer.bands_end());
|
||||
for (int band = analyzer.bandpass_bands_begin(); band < analyzer.bandpass_bands_end(); band++) {
|
||||
float f_hz = analyzer.band_ff(band) * fs;
|
||||
double f_hz = analyzer.band_ff(band) * fs;
|
||||
band_gains[band] = 1.0 / sqrt(f_hz / 20.0);
|
||||
}
|
||||
band_gains[analyzer.band_lowpass()] = band_gains[analyzer.bandpass_bands_end() - 1];
|
||||
|
|
|
|||
42
gaborator/affine_transform.h
Normal file
42
gaborator/affine_transform.h
Normal file
|
|
@ -0,0 +1,42 @@
|
|||
//
|
||||
// A class for affine transforms (ax + b) of scalar values
|
||||
//
|
||||
// Copyright (C) 2020-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
||||
#ifndef _GABORATOR_AFFINE_TRANSFORM_H
|
||||
#define _GABORATOR_AFFINE_TRANSFORM_H
|
||||
|
||||
namespace gaborator {
|
||||
|
||||
struct affine_transform {
|
||||
affine_transform(): a(0), b(0) { }
|
||||
affine_transform(double a_, double b_): a(a_), b(b_) { }
|
||||
affine_transform(const affine_transform &rhs): a(rhs.a), b(rhs.b) { }
|
||||
double operator()(double x) const { return a * x + b; }
|
||||
affine_transform inverse() const {
|
||||
return affine_transform(1.0 / a, -b / a);
|
||||
}
|
||||
static affine_transform identity() { return affine_transform(1, 0); }
|
||||
double a, b;
|
||||
};
|
||||
|
||||
// Composition
|
||||
|
||||
static inline affine_transform
|
||||
operator *(const affine_transform &a, const affine_transform &b) {
|
||||
return affine_transform(a.a * b.a, a.a * b.b + a.b);
|
||||
}
|
||||
|
||||
// Equality
|
||||
|
||||
static inline bool
|
||||
operator ==(const affine_transform &a, const affine_transform &b) {
|
||||
return a.a == b.a && a.b == b.b;
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
#endif
|
||||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// Fast Fourier transform
|
||||
//
|
||||
// Copyright (C) 2016-2020 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2016-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -14,13 +14,16 @@
|
|||
#if GABORATOR_USE_VDSP
|
||||
#include "gaborator/fft_vdsp.h"
|
||||
#define GABORATOR_USE_REAL_FFT 1
|
||||
#define GABORATOR_MIN_FFT_SIZE 1
|
||||
#elif GABORATOR_USE_PFFFT
|
||||
#include "gaborator/fft_pffft.h"
|
||||
#define GABORATOR_USE_REAL_FFT 1
|
||||
#define GABORATOR_MIN_FFT_SIZE 32
|
||||
#else
|
||||
// Use the naive FFT
|
||||
// Do not define GABORATOR_USE_REAL_FFT as it is slower than
|
||||
// using the complex code.
|
||||
#define GABORATOR_MIN_FFT_SIZE 1
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// Fast Fourier transform, naive reference implementations
|
||||
//
|
||||
// Copyright (C) 1992-2020 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 1992-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// Fast Fourier transform using PFFFT
|
||||
//
|
||||
// Copyright (C) 2017-2020 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2017-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -76,14 +76,14 @@ private:
|
|||
PFFFT_Setup *setup;
|
||||
};
|
||||
|
||||
// Support transforming std::vector<std::complex<float> >::iterator
|
||||
// Support transforming std::vector<std::complex<float>>::iterator
|
||||
|
||||
template <>
|
||||
struct fft<std::vector<std::complex<float> >::iterator>:
|
||||
struct fft<std::vector<std::complex<float>>::iterator>:
|
||||
public fft<std::complex<float> *>
|
||||
{
|
||||
typedef fft<std::complex<float> *> base;
|
||||
typedef std::vector<std::complex<float> >::iterator I;
|
||||
typedef std::vector<std::complex<float>>::iterator I;
|
||||
fft(unsigned int n_): fft<std::complex<float> *>(n_) { }
|
||||
void
|
||||
transform(I a) {
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// Fast Fourier transform using the Apple vDSP framework
|
||||
//
|
||||
// Copyright (C) 2013-2020 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2013-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -177,14 +177,14 @@ private:
|
|||
FFTSetup setup;
|
||||
};
|
||||
|
||||
// Support transforming std::vector<std::complex<float> >::iterator
|
||||
// Support transforming std::vector<std::complex<float>>::iterator
|
||||
|
||||
template <>
|
||||
struct fft<typename std::vector<std::complex<float> >::iterator>:
|
||||
struct fft<typename std::vector<std::complex<float>>::iterator>:
|
||||
public fft<std::complex<float> *>
|
||||
{
|
||||
typedef fft<std::complex<float> *> base;
|
||||
typedef typename std::vector<std::complex<float> >::iterator I;
|
||||
typedef typename std::vector<std::complex<float>>::iterator I;
|
||||
fft(unsigned int n_): fft<std::complex<float> *>(n_) { }
|
||||
void
|
||||
transform(I a) {
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// Constant Q spectrum analysis and resynthesis
|
||||
//
|
||||
// Copyright (C) 2015-2020 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2015-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -27,7 +27,7 @@
|
|||
|
||||
#include "gaborator/fft.h"
|
||||
#include "gaborator/gaussian.h"
|
||||
#include "gaborator/linear_transform.h"
|
||||
#include "gaborator/affine_transform.h"
|
||||
#include "gaborator/pod_vector.h"
|
||||
#include "gaborator/pool.h"
|
||||
#include "gaborator/ref.h"
|
||||
|
|
@ -187,19 +187,6 @@ private:
|
|||
}
|
||||
|
||||
public:
|
||||
// Note: Pointer returned becomes invalid when range_vector
|
||||
// is changed
|
||||
T *
|
||||
get(I i, bool create) {
|
||||
if (! has_index(i)) {
|
||||
if (create)
|
||||
extend(i);
|
||||
else
|
||||
return 0;
|
||||
}
|
||||
return unchecked_get(i);
|
||||
}
|
||||
|
||||
// Get a pointer to an existing element, or null if out of range
|
||||
const T *
|
||||
get(I i) const {
|
||||
|
|
@ -212,7 +199,10 @@ public:
|
|||
// is changed
|
||||
T &
|
||||
get_or_create(I i) {
|
||||
return *get(i, true);
|
||||
if (! has_index(i))
|
||||
extend(i);
|
||||
return *unchecked_get(i);
|
||||
|
||||
}
|
||||
|
||||
// Get a reference to the element at index i, which must be valid
|
||||
|
|
@ -299,27 +289,37 @@ static inline unsigned fat_part(unsigned int fftsize) {
|
|||
return fftsize >> 2;
|
||||
}
|
||||
|
||||
// Frequency band parameters shared between octaves
|
||||
// Per-band, per-plan data
|
||||
|
||||
template<class T>
|
||||
struct band_params: public refcounted {
|
||||
struct band_plan {
|
||||
typedef complex<T> C;
|
||||
bool dc; // True iff this is the DC band
|
||||
unsigned int sftsize; // Size of "short FFT" spanning the band
|
||||
unsigned int sftsize_log2; // log2(sftsize)
|
||||
unsigned int step; // Signal samples per coefficient sample
|
||||
unsigned int step_log2; // log2(step)
|
||||
fft<C *> *sft; // Fourier transform for windows, of size sftsize
|
||||
std::vector<T> kernel; // Frequency-domain filter kernel
|
||||
std::vector<T> dual_kernel; // Dual of the above
|
||||
pod_vector<C> shift_kernel; // Complex exponential for fractional frequency compensation
|
||||
pod_vector<C> shift_kernel_conj; // Conjugate of the above
|
||||
int fq_offset_int; // Frequency offset in bins (big-fft bin of left window edge)
|
||||
double ff; // Center (bp) or corner (lp) frequency in units of the sampling frequency
|
||||
int fq_offset_int; // Frequency offset in bins (big-FFT bin of left window edge)
|
||||
double center; // Center frequency in units of FFT bins
|
||||
int icenter; // Center frequency rounded to nearest integer FFT bin
|
||||
};
|
||||
|
||||
// Frequency band parameters shared between octaves
|
||||
|
||||
template<class T>
|
||||
struct band_params: public refcounted {
|
||||
typedef complex<T> C;
|
||||
bool dc; // True iff this is the lowpass (DC) band
|
||||
double ff; // Center (bp) or corner (lp) frequency in units of the sampling frequency
|
||||
double ffsd; // Standard deviation of the bandpass Gaussian, as fractional frequency
|
||||
unsigned int time_support; // Filter support in time domain, in octave subsamples
|
||||
unsigned int step; // Signal samples per coefficient sample
|
||||
unsigned int step_log2; // log2(step)
|
||||
double ff_support; // Filter support in frequency domain
|
||||
double time_support; // Filter support in time domain, in octave subsamples
|
||||
std::vector<band_plan<T>> anl_plans;
|
||||
std::vector<band_plan<T>> syn_plans;
|
||||
};
|
||||
|
||||
// Downsampling parameters. These have some similarity to band
|
||||
|
|
@ -342,14 +342,14 @@ struct downsampling_params {
|
|||
// Forward declarations
|
||||
template<class T> struct analyzer;
|
||||
template<class T> struct zone;
|
||||
template<class T, class C = complex<T> > struct coefs;
|
||||
template<class T, class C = complex<T>> struct coefs;
|
||||
template<class C> struct sliced_coefs;
|
||||
template <class T, class OI = complex<T> *, class C = complex<T> >
|
||||
struct row_source;
|
||||
template <class T, class II = complex<T> *, class C = complex<T> >
|
||||
struct row_dest;
|
||||
template <class T, class II = complex<T> *, class C = complex<T> >
|
||||
struct row_add_dest;
|
||||
template <class T, class OI = complex<T> *, class C = complex<T>>
|
||||
struct row_source;
|
||||
template <class T, class II = complex<T> *, class C = complex<T>>
|
||||
struct row_dest;
|
||||
template <class T, class II = complex<T> *, class C = complex<T>>
|
||||
struct row_add_dest;
|
||||
|
||||
// Abstract class for tracking changes to a coefficient set.
|
||||
// This may be used for updating a resolution pyramid of
|
||||
|
|
@ -358,10 +358,10 @@ template <class T, class II = complex<T> *, class C = complex<T> >
|
|||
template<class T>
|
||||
struct shadow {
|
||||
virtual ~shadow() { }
|
||||
virtual void update(const coefs<complex<T> > &msc,
|
||||
virtual void update(const coefs<complex<T>> &msc,
|
||||
sample_index_t i0, sample_index_t i1) = 0;
|
||||
// Deprecated
|
||||
virtual void update(int oct, int ze, const sliced_coefs<complex<T> > &sc,
|
||||
virtual void update(int oct, int ze, const sliced_coefs<complex<T>> &sc,
|
||||
slice_index_t sli0, slice_index_t sli1) = 0;
|
||||
};
|
||||
|
||||
|
|
@ -388,9 +388,9 @@ struct zone_coefs_meta {
|
|||
void init(const band_vector &bands_) {
|
||||
bands = bands_;
|
||||
unsigned int offset = 0;
|
||||
for (unsigned int i = 0; i < bands.size(); i++) {
|
||||
bands[i].band_offset = offset;
|
||||
offset += bands[i].slice_len;
|
||||
for (band_coefs_meta &b: bands) {
|
||||
b.band_offset = offset;
|
||||
offset += b.slice_len;
|
||||
}
|
||||
total_size = offset;
|
||||
}
|
||||
|
|
@ -452,24 +452,20 @@ bno_split(const coefs_meta &meta, int gbno, int &oct, unsigned int &obno, bool d
|
|||
} else {
|
||||
// Within bandpass region
|
||||
// Start by determining the octave
|
||||
int oct_tmp = gbno / meta.bands_per_octave;
|
||||
int obno_tmp = gbno % meta.bands_per_octave;
|
||||
// Octave 0 is twice the size
|
||||
if (oct_tmp < 2) {
|
||||
if (oct_tmp == 1) {
|
||||
// Octave 1 gets appended to octave 0
|
||||
obno_tmp += meta.bands_per_octave;
|
||||
oct_tmp--;
|
||||
}
|
||||
} else {
|
||||
// Because octave 0 is twice the size, the
|
||||
// other octave numbers are one less
|
||||
oct_tmp--;
|
||||
int n_bands_top_octave = (int)meta.octaves[0].z->bands.size();
|
||||
if (gbno < n_bands_top_octave) {
|
||||
// Top octave
|
||||
oct = 0;
|
||||
obno = n_bands_top_octave - 1 - gbno;
|
||||
return true;
|
||||
}
|
||||
gbno -= n_bands_top_octave;
|
||||
int oct_tmp = 1 + gbno / meta.bands_per_octave;
|
||||
int obno_tmp = gbno % meta.bands_per_octave;
|
||||
oct = oct_tmp;
|
||||
// Now determine the band within the octave.
|
||||
// obno_tmp counts down, but obno counts up.
|
||||
obno = (unsigned int)meta.octaves[oct].z->bands.size() - 1 - obno_tmp;
|
||||
obno = (unsigned int)meta.octaves[oct_tmp].z->bands.size() - 1 - obno_tmp;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
|
@ -553,12 +549,12 @@ void add(oct_coefs<C> &a, const oct_coefs<C> &b) {
|
|||
|
||||
template<class C>
|
||||
struct sliced_coefs {
|
||||
typedef range_vector<ref<oct_coefs<C> >, slice_index_t> slices_t;
|
||||
typedef range_vector<ref<oct_coefs<C>>, slice_index_t> slices_t;
|
||||
uint64_t estimate_memory_usage() const {
|
||||
unsigned int n = 0;
|
||||
size_t size_each = 0;
|
||||
for (slice_index_t sl = slices.begin_index(); sl < slices.end_index(); sl++) {
|
||||
const ref<oct_coefs<C> > &t = slices.get_existing(sl);
|
||||
const ref<oct_coefs<C>> &t = slices.get_existing(sl);
|
||||
if (t) {
|
||||
if (! size_each)
|
||||
size_each = (size_t)t->estimate_memory_usage();
|
||||
|
|
@ -597,7 +593,7 @@ oct_coefs<CT> *get_existing_coefs(const sliced_coefs<CT> &sc,
|
|||
|
||||
template <class CT>
|
||||
oct_coefs<CT> &get_or_create_coefs(sliced_coefs<CT> &sc, slice_index_t i) {
|
||||
ref<oct_coefs<CT> > &p(sc.slices.get_or_create(i));
|
||||
ref<oct_coefs<CT>> &p(sc.slices.get_or_create(i));
|
||||
if (! p)
|
||||
p.reset(new oct_coefs<CT>(*sc.meta));
|
||||
return *p;
|
||||
|
|
@ -652,8 +648,7 @@ struct coefs {
|
|||
coefs(const analyzer<T> &anl_, shadow<T> *shadow_ = 0):
|
||||
octaves(anl_.n_octaves), shadow0(shadow_)
|
||||
{
|
||||
// Use the analysis structure of the largest plan
|
||||
meta = anl_.anl_plans[anl_.anl_plans.size() - 1]->cmeta;
|
||||
meta = anl_.cmeta_any.get();
|
||||
// Set up shortcut pointer to zone metadata in each octave
|
||||
for (unsigned int oct = 0; oct < octaves.size(); oct++)
|
||||
octaves[oct].meta = meta->octaves[oct].z;
|
||||
|
|
@ -668,8 +663,8 @@ struct coefs {
|
|||
for (unsigned int oct = 0; oct < octaves.size(); oct++)
|
||||
octaves[oct].clear();
|
||||
}
|
||||
ref<coefs_meta> meta;
|
||||
std::vector<sliced_coefs<C> > octaves;
|
||||
coefs_meta *meta;
|
||||
std::vector<sliced_coefs<C>> octaves;
|
||||
shadow<T> *shadow0;
|
||||
};
|
||||
|
||||
|
|
@ -775,12 +770,30 @@ static inline int signed_index(unsigned int i, unsigned int size) {
|
|||
return (i ^ t) - t;
|
||||
}
|
||||
|
||||
// Construct a frequency-domain lowpass filter whose response is the
|
||||
// convolution of a rectangle and a gaussian. The cutoff freqency is
|
||||
// ff_cutoff (a fractional frequency), and the standard deviation of
|
||||
// the gaussian is ff_sd. The returned filter covers the full
|
||||
// frequency range from 0 to fs (with negative frequencies at the end,
|
||||
// the usual convention for FFT spectra).
|
||||
// Evaluate a Gaussian windowed lowpass filter frequency response.
|
||||
// This is the convolution of a rectangle centered at f=0 and a Gaussian,
|
||||
// and corresponds to a Gaussian windowed sinc in the time domain.
|
||||
// The -6 dB cutoff freqency is ff_cutoff (a fractional frequency),
|
||||
// the standard deviation of the Gaussian is ff_sd, and the frequency
|
||||
// response is evaluated at ff. The frequency response is smooth at
|
||||
// f=0 even if the transition bands overlap.
|
||||
|
||||
inline double
|
||||
gaussian_windowed_lowpass_1(double ff_cutoff, double ff_sd, double ff) {
|
||||
return
|
||||
// A rectangle is the sum of a rising step and a later falling
|
||||
// step, or the difference between a rising step and a later
|
||||
// rising step. By linearity, a Gaussian filtered rectangle
|
||||
// is the difference between two Gaussian filtered rising
|
||||
// steps.
|
||||
gaussian_edge(ff_sd, -ff + ff_cutoff) -
|
||||
gaussian_edge(ff_sd, -ff - ff_cutoff);
|
||||
}
|
||||
|
||||
// Fill a sequence with a frequency-ddomain lowpass filter as above.
|
||||
// The returned filter covers the full frequency range from 0 to fs
|
||||
// (with negative frequencies at the end, the usual convention for FFT
|
||||
// spectra).
|
||||
//
|
||||
// When center=true, construct a time-domain window instead,
|
||||
// passing the center of the time-domain signal.
|
||||
|
|
@ -803,9 +816,7 @@ inline void gaussian_windowed_lowpass(double ff_cutoff, double ff_sd,
|
|||
else
|
||||
// Symmetric around zero
|
||||
thisff = (i > len / 2 ? len - i : i) * inv_len;
|
||||
double x = thisff - ff_cutoff;
|
||||
double v = gaussian_edge(ff_sd, -x);
|
||||
*it = v;
|
||||
*it = gaussian_windowed_lowpass_1(ff_cutoff, ff_sd, thisff);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -820,15 +831,11 @@ struct zone: public refcounted {
|
|||
unsigned int zno;
|
||||
// Total number of bands, including DC band if lowest octave
|
||||
unsigned int n_bands;
|
||||
unsigned int max_step_log2;
|
||||
// Band parameters by increasing frequency; DC band is index 0 if
|
||||
// present
|
||||
std::vector<ref<band_params<T> > > bandparams;
|
||||
// Pseudo-bands mimicing the response of bands in the
|
||||
// neighboring octaves, used for calculating the duals only
|
||||
std::vector<ref<band_params<T> > > mock_bandparams;
|
||||
pod_vector<T> power;
|
||||
pod_vector<T> power_nodc;
|
||||
unsigned int max_step_log2;
|
||||
std::vector<ref<band_params<T>>> bandparams;
|
||||
std::vector<ref<band_params<T>>> mock_bandparams;
|
||||
};
|
||||
|
||||
template <class T>
|
||||
|
|
@ -839,38 +846,68 @@ struct octave {
|
|||
// Helper function for pushing parameters onto the vectors in struct zone
|
||||
|
||||
template <class T>
|
||||
void push(std::vector<ref<band_params<T> > > &v, band_params<T> *p) {
|
||||
v.push_back(ref<band_params<T> >(p));
|
||||
void push(std::vector<ref<band_params<T>>> &v, band_params<T> *p) {
|
||||
v.push_back(ref<band_params<T>>(p));
|
||||
}
|
||||
|
||||
// Phase conventions: phase_convention::absolute means the phase of a
|
||||
// coefficient at time tc is relative to e(i tau f t), and
|
||||
// phase_convention::relative means it is relative to
|
||||
// e(i tau f (t - tc))
|
||||
// Phase conventions: coef_phase::absolute means the phase of a
|
||||
// coefficient at time tc is relative to e^(i tau f t), and
|
||||
// coef_phase::local means it is relative to
|
||||
// e^(i tau f (t - tc))
|
||||
|
||||
enum class phase_convention { absolute, relative };
|
||||
|
||||
enum class bandwidth_formula { v1, v2 };
|
||||
enum class coef_phase { global, local };
|
||||
|
||||
// A set of spectrum analysis parameters
|
||||
|
||||
struct parameters {
|
||||
parameters(unsigned int bands_per_octave_, double ff_min_,
|
||||
parameters(unsigned int bands_per_octave_,
|
||||
double ff_min_,
|
||||
double ff_ref_ = 1.0,
|
||||
double overlap_ = 0.7,
|
||||
double max_error_ = 1e-5,
|
||||
phase_convention phase_ = phase_convention::absolute,
|
||||
bandwidth_formula bandwidth_ = bandwidth_formula::v1):
|
||||
double max_error_ = 1e-5):
|
||||
bands_per_octave(bands_per_octave_),
|
||||
ff_min(ff_min_),
|
||||
ff_ref(ff_ref_),
|
||||
overlap(overlap_),
|
||||
max_error(max_error_),
|
||||
phase(phase_),
|
||||
bandwidth(bandwidth_),
|
||||
coef_scale(1.0),
|
||||
synthesis(true)
|
||||
{ }
|
||||
synthesis(true),
|
||||
multirate(true)
|
||||
{
|
||||
init_v1();
|
||||
}
|
||||
// Pseudo-constructor with version 1 defaults
|
||||
static parameters v1(unsigned int bands_per_octave_,
|
||||
double ff_min_,
|
||||
double ff_ref_ = 1.0,
|
||||
double overlap_ = 0.7,
|
||||
double max_error_ = 1e-5)
|
||||
{
|
||||
parameters p(bands_per_octave_, ff_min_, ff_ref_, overlap_, max_error_);
|
||||
p.init_v1();
|
||||
return p;
|
||||
}
|
||||
// Pseudo-constructor with version 2 defaults
|
||||
static parameters v2(unsigned int bands_per_octave_,
|
||||
double ff_min_,
|
||||
double ff_ref_ = 1.0,
|
||||
double overlap_ = 0.7,
|
||||
double max_error_ = 1e-5)
|
||||
{
|
||||
parameters p(bands_per_octave_, ff_min_, ff_ref_, overlap_, max_error_);
|
||||
p.init_v2();
|
||||
return p;
|
||||
}
|
||||
void init_v1() {
|
||||
phase = coef_phase::global;
|
||||
bandwidth_version = 1;
|
||||
lowpass_version = 1;
|
||||
}
|
||||
void init_v2() {
|
||||
phase = coef_phase::local;
|
||||
bandwidth_version = 2;
|
||||
lowpass_version = 2;
|
||||
}
|
||||
// Provide an operator< so that we can create a set or map of parameters
|
||||
bool operator<(const parameters &b) const {
|
||||
#define GABORATOR_COMPARE_LESS(member) do { \
|
||||
|
|
@ -885,6 +922,11 @@ struct parameters {
|
|||
GABORATOR_COMPARE_LESS(overlap);
|
||||
GABORATOR_COMPARE_LESS(max_error);
|
||||
GABORATOR_COMPARE_LESS(phase);
|
||||
GABORATOR_COMPARE_LESS(bandwidth_version);
|
||||
GABORATOR_COMPARE_LESS(lowpass_version);
|
||||
GABORATOR_COMPARE_LESS(coef_scale);
|
||||
GABORATOR_COMPARE_LESS(synthesis);
|
||||
GABORATOR_COMPARE_LESS(multirate);
|
||||
#undef GABORATOR_COMPARE_LESS
|
||||
// Equal
|
||||
return false;
|
||||
|
|
@ -903,7 +945,7 @@ struct parameters {
|
|||
// The standard deviation of the Gaussian in units of the mean
|
||||
double sd() const {
|
||||
return overlap *
|
||||
(bandwidth == bandwidth_formula::v1 ?
|
||||
(bandwidth_version == 1 ?
|
||||
band_spacing() - 1 :
|
||||
log(2) / bands_per_octave);
|
||||
}
|
||||
|
|
@ -926,10 +968,12 @@ struct parameters {
|
|||
double ff_ref;
|
||||
double overlap;
|
||||
double max_error;
|
||||
phase_convention phase;
|
||||
bandwidth_formula bandwidth;
|
||||
coef_phase phase;
|
||||
int bandwidth_version;
|
||||
int lowpass_version;
|
||||
double coef_scale;
|
||||
bool synthesis; // Synthesis is supported
|
||||
bool multirate;
|
||||
};
|
||||
|
||||
// Like std::fill, but returns the end iterator
|
||||
|
|
@ -945,19 +989,20 @@ I fill(I b, I e, T v) {
|
|||
|
||||
template <class V, class S>
|
||||
void scale_vector(V &v, S s) {
|
||||
for (size_t i = 0; i < v.size(); i++)
|
||||
v[i] *= s;
|
||||
for (auto &e: v)
|
||||
e *= s;
|
||||
}
|
||||
|
||||
// Zero-padding source wrapper. This returns data from the underlying
|
||||
// source within the interval src_i0 to src_i1, and zero elsewhere.
|
||||
|
||||
template <class S, class T>
|
||||
template <class S, class OI>
|
||||
struct zeropad_source {
|
||||
typedef typename std::iterator_traits<OI>::value_type T;
|
||||
zeropad_source(const S &source_, int64_t src_i0_, int64_t src_i1_):
|
||||
source(source_), src_i0(src_i0_), src_i1(src_i1_)
|
||||
{ }
|
||||
T *operator()(int64_t i0, int64_t i1, T *output) const {
|
||||
OI operator()(int64_t i0, int64_t i1, OI output) const {
|
||||
int64_t overlap_begin = std::max(i0, src_i0);
|
||||
int64_t overlap_end = std::min(i1, src_i1);
|
||||
if (overlap_end <= overlap_begin) {
|
||||
|
|
@ -1002,16 +1047,15 @@ void copy_overlapping_zerofill(T *dst, size_t dstlen, const T *src,
|
|||
int64_t i0, int64_t i1)
|
||||
{
|
||||
pointer_source<T> ps(src, i0, i1);
|
||||
zeropad_source<pointer_source<T>, T> zs(ps, i0, i1);
|
||||
zeropad_source<pointer_source<T>, T *> zs(ps, i0, i1);
|
||||
zs(0, dstlen, dst);
|
||||
}
|
||||
|
||||
// Given a set of FFT coefficients "coefs" of a real
|
||||
// sequence, where only positive-frequency coefficients
|
||||
// (including DC and Nyquist) are valid, return the
|
||||
// coefficient for an arbitrary frequency index "i"
|
||||
// which may correspond to a negative frequency, or
|
||||
// even an alias outside the range (0..fftsize-1).
|
||||
// Given a set of FFT coefficients "coefs" of a real sequence, where
|
||||
// only positive-frequency coefficients (including DC and Nyquist) are
|
||||
// valid, return the coefficient for an arbitrary frequency index "i"
|
||||
// which may correspond to a negative frequency, or even an alias
|
||||
// outside the range (0..fftsize-1).
|
||||
|
||||
template<class T>
|
||||
complex<T> get_real_spectrum_coef(complex<T> *coefs, int i, unsigned int fftsize) {
|
||||
|
|
@ -1100,6 +1144,12 @@ void get_band_coef_bounds(const coefs<T, C> &msc, int oct, unsigned int obno,
|
|||
{
|
||||
const sliced_coefs<C> &sc = msc.octaves[oct];
|
||||
const typename sliced_coefs<C>::slices_t &slices = sc.slices;
|
||||
if (slices.empty()) {
|
||||
// Don't try to convert int64t_min/max slices to coef time
|
||||
ci0_ret = 0;
|
||||
ci1_ret = 0;
|
||||
return;
|
||||
}
|
||||
// Convert from slices to coefficient samples
|
||||
ci0_ret = coef_time(*sc.meta, slices.begin_index(), oct, obno);
|
||||
ci1_ret = coef_time(*sc.meta, slices.end_index(), oct, obno);
|
||||
|
|
@ -1116,6 +1166,38 @@ void get_band_coef_bounds(const coefs<T, C> &msc, int gbno,
|
|||
get_band_coef_bounds(msc, oct, obno, ci0_ret, ci1_ret);
|
||||
}
|
||||
|
||||
// Evaluate the frequency-domain analysis filter kernel of band "bp"
|
||||
// at frequency "ff"
|
||||
|
||||
template <class T>
|
||||
double eval_kernel(parameters *, band_params<T> *bp, double ff) {
|
||||
if (bp->dc) {
|
||||
return gaussian_windowed_lowpass_1(bp->ff, bp->ffsd, ff);
|
||||
} else {
|
||||
return norm_gaussian(bp->ffsd, ff - bp->ff);
|
||||
}
|
||||
}
|
||||
|
||||
// Evaluate the frequency-domain synthesis filter kernel of band "bp"
|
||||
// at frequency "ff"
|
||||
|
||||
template <class T>
|
||||
double eval_dual_kernel(parameters *params, band_params<T> *bp, double ff) {
|
||||
double gain = 1.0;
|
||||
if (params->lowpass_version == 2) {
|
||||
if (bp->dc) {
|
||||
// Adjust the gain of the reconstruction lowpass
|
||||
// filter to make the overall gain similar to to
|
||||
// the bandpass region.
|
||||
double adjusted_overlap = params->sd() /
|
||||
(log(2) / params->bands_per_octave);
|
||||
double avg_bandpass_gain = adjusted_overlap * sqrt(M_PI);
|
||||
gain = avg_bandpass_gain * 0.5;
|
||||
}
|
||||
}
|
||||
return eval_kernel(params, bp, ff) * gain;
|
||||
}
|
||||
|
||||
template<class T>
|
||||
struct analyzer: public refcounted {
|
||||
typedef complex<T> C;
|
||||
|
|
@ -1123,6 +1205,7 @@ struct analyzer: public refcounted {
|
|||
|
||||
analyzer(const parameters ¶ms_):
|
||||
params(params_),
|
||||
max_step_log2(0),
|
||||
fftsize_max(0),
|
||||
sftsize_max(0)
|
||||
{
|
||||
|
|
@ -1141,12 +1224,6 @@ struct analyzer: public refcounted {
|
|||
// where one band would fall on fs exactly.
|
||||
tuning_log2ff = sane_fmod(log2(params.ff_ref), band_spacing_log2);
|
||||
|
||||
// Calculate the frequency of band0, the "fs/8 band", the
|
||||
// lowest band in each octave except possibly the lowest octave.
|
||||
// Its frequency will fs/8, or a fraction of the band spacing higher
|
||||
// due to tuning. The magic -3 derives from fs/8 as log2(1/8).
|
||||
band0_log2ff = -3 + tuning_log2ff;
|
||||
|
||||
// Calculate the total number of bands needed so that
|
||||
// the lowest band has a frequency <= params.ff_min.
|
||||
// end_log2ff = the log2ff of the band after the last (just past fs/2);
|
||||
|
|
@ -1162,7 +1239,7 @@ struct analyzer: public refcounted {
|
|||
ffref_gbno = (int)rint((top_band_log2ff - log2(params.ff_ref)) /
|
||||
band_spacing_log2);
|
||||
|
||||
// Establish linar transforms for converting between
|
||||
// Establish affine transforms for converting between
|
||||
// log-frequencies (log2(ff)) and bandpass band numbers.
|
||||
// Derivation:
|
||||
//ff = exp2(tuning_log2ff - 1 - (gbno + 1) * band_spacing_log2)
|
||||
|
|
@ -1181,41 +1258,40 @@ struct analyzer: public refcounted {
|
|||
// negative of an unsigned int.
|
||||
double a = -(double)params.bands_per_octave;
|
||||
double b = -a * tuning_log2ff + a - 1;
|
||||
log2ff_bandpass_band = linear_transform(a, b);
|
||||
log2ff_bandpass_band = affine_transform(a, b);
|
||||
bandpass_band_log2ff = log2ff_bandpass_band.inverse();
|
||||
|
||||
// Calculate the kernel support needed for band0,
|
||||
// and size the FFT accordingly. This duplicates some
|
||||
// code at the beginning of make_band().
|
||||
double band0_ff = exp2(band0_log2ff);
|
||||
double band0_time_sd = time_sd(band0_ff);
|
||||
|
||||
// Figure out the number of octaves, keeping in mind that
|
||||
// the top octave is twice the size, and DC band is added
|
||||
// to the bottom octave even if that makes it larger than
|
||||
// the others.
|
||||
n_octaves = (n_bandpass_bands_total + params.bands_per_octave - 1) /
|
||||
params.bands_per_octave;
|
||||
if (n_octaves > 1)
|
||||
n_octaves--;
|
||||
|
||||
{
|
||||
// Precalculate the parameters of the downsampling filter.
|
||||
// These are the same for all plans, and need to be
|
||||
// calculated before creating the plans; in particular, we
|
||||
// need to know the support before we can create the
|
||||
// plans, because in low-bpo cases, it can determine the
|
||||
// minimum amount of fat needed. The actual filter
|
||||
// coefficients are plan-dependent and will be calculated
|
||||
// later.
|
||||
// minimum amount of fat needed. The filter kernel is
|
||||
// specific to the plan as it depends on the FFT size,
|
||||
// and will be calculated later.
|
||||
|
||||
// In terms of the post-downsampling sampling frequency,
|
||||
// the downsampling lowpass filter transition band starts
|
||||
// at the top edge of the highest band, and ranges to
|
||||
// fs/2. We use the worst-case center frequency for the
|
||||
// highest band, fs/4, and add the support.
|
||||
double f0 = 0.25 + gaussian_support(ff_sd(0.25), params.max_error);
|
||||
// When operating at a high Q, the we will need to use
|
||||
// large FFTs in any case, and it makes sense to use a
|
||||
// narrow transition band because we can get that
|
||||
// essentially for free, and the passband will be
|
||||
// correspondingly wider, which will allow processing more
|
||||
// bands at the lower sample rate. Conversely, at low
|
||||
// Q, we should use a wide transition band so that the
|
||||
// FFTs can be kept short.
|
||||
|
||||
// The filter is defined in terms of the lower
|
||||
// (downsampled) sample rate.
|
||||
|
||||
// Make the transition band the same width as the width
|
||||
// (two-sided support) of a band at ff=0.25, but don't let
|
||||
// the low edge go below 0.25 to make sure we have a
|
||||
// reasonable amount of passband left.
|
||||
double f1 = 0.5;
|
||||
double f0 =
|
||||
std::max(f1 - 2 * gaussian_support(ff_sd(0.25), params.max_error),
|
||||
0.25);
|
||||
assert(f0 < f1);
|
||||
|
||||
// The cutoff frequency is in the center of the transition band
|
||||
double ff = (f0 + f1) * 0.5;
|
||||
|
|
@ -1227,39 +1303,77 @@ struct analyzer: public refcounted {
|
|||
double time_sd = sd_f2t(ff_sd);
|
||||
|
||||
// Set members
|
||||
ds_passband = f0;
|
||||
ds_ff = ff;
|
||||
ds_ff_sd = ff_sd;
|
||||
ds_time_support = gaussian_support(time_sd, params.max_error);
|
||||
// Since the filter is designed at the lower sample rate,
|
||||
// ds_time_support is in the unit of lower octave samples
|
||||
ds_time_support = gaussian_support(time_sd, params.max_error * 0.5);
|
||||
}
|
||||
|
||||
band0_time_analysis_support =
|
||||
gaussian_support(band0_time_sd, params.max_error);
|
||||
band0_time_synthesis_support =
|
||||
band0_time_analysis_support * synthesis_support_multiplier();
|
||||
// Determine the octave structure, packing each band into the
|
||||
// lowest octave possible. For now, while bpo is restricted
|
||||
// to integer values, this just means determining how many
|
||||
// bands need to go in the top octave, and the remaining ones
|
||||
// will be divided into groups of bpo bands (except possibly
|
||||
// the last).
|
||||
int gbno;
|
||||
for (gbno = bandpass_bands_begin(); gbno < bandpass_bands_end(); gbno++) {
|
||||
double ff = bandpass_band_ff(gbno);
|
||||
double ffsd = ff_sd(ff);
|
||||
double ff_support = gaussian_support(ffsd, params.max_error * 0.5);
|
||||
// If the bandpass support falls within the downsampling filter
|
||||
// passband of the next octave, we can switch octaves.
|
||||
if (params.multirate && ff + ff_support <= ds_passband / 2)
|
||||
break;
|
||||
}
|
||||
n_bands_top_octave = gbno;
|
||||
|
||||
// Figure out the number of octaves, keeping in mind that
|
||||
// the top octave is of variable size, and DC band is added
|
||||
// to the bottom octave even if that makes it larger than
|
||||
// the others.
|
||||
n_octaves = 1 // The top octave
|
||||
+ (n_bandpass_bands_total - n_bands_top_octave +
|
||||
(params.bands_per_octave - 1)) / params.bands_per_octave;
|
||||
|
||||
// Calculate the kernel support needed for the lowest-frequency
|
||||
// bandpass band to use as a basis for an initial estimate of
|
||||
// the FFT size needed. This duplicates some code at the
|
||||
// beginning of make_band().
|
||||
assert(n_bands_top_octave >= 1);
|
||||
int low_bp_band = n_bands_top_octave - 1;
|
||||
double low_bp_band_time_sd = time_sd(band_ff(low_bp_band));
|
||||
|
||||
double low_bp_band_time_analysis_support =
|
||||
gaussian_support(low_bp_band_time_sd, params.max_error);
|
||||
double low_bp_band_time_synthesis_support =
|
||||
low_bp_band_time_analysis_support * synthesis_support_multiplier();
|
||||
|
||||
make_zones();
|
||||
|
||||
// Make analysis plans
|
||||
// Since ds_time_support is in the unit of lower octave samples,
|
||||
// we need to multiply it by two to get upper octave samples.
|
||||
unsigned int max_support =
|
||||
std::max(ceil(band0_time_analysis_support),
|
||||
std::max(ceil(low_bp_band_time_analysis_support),
|
||||
ds_time_support * 2);
|
||||
|
||||
unsigned int size = next_power_of_two(max_support * 2);
|
||||
plan *p;
|
||||
ref<plan> p;
|
||||
for (;;) {
|
||||
p = new plan(this, false, size, max_support);
|
||||
if (p->ok)
|
||||
break;
|
||||
delete p;
|
||||
size *= 2;
|
||||
}
|
||||
anl_plans.push_back(p); // Smallest possible plan
|
||||
size *= 2;
|
||||
p = new plan(this, false, size, max_support);
|
||||
p = new plan(this, false, size * 2, max_support);
|
||||
assert(p->ok);
|
||||
anl_plans.push_back(p); // Next larger plan
|
||||
|
||||
if (params.synthesis) {
|
||||
// Make synthesis plan (only one for now)
|
||||
max_support = std::max(ceil(band0_time_synthesis_support),
|
||||
max_support = std::max(ceil(low_bp_band_time_synthesis_support),
|
||||
ds_time_support * 2);
|
||||
// Room for at at least the two fats + as much filet
|
||||
size = next_power_of_two(max_support * 2) * 2;
|
||||
|
|
@ -1268,6 +1382,11 @@ struct analyzer: public refcounted {
|
|||
syn_plans.push_back(p);
|
||||
}
|
||||
|
||||
for (int i = 0; i < (int)anl_plans.size(); i++)
|
||||
make_band_plans(i, false);
|
||||
for (int i = 0; i < (int)syn_plans.size(); i++)
|
||||
make_band_plans(i, true);
|
||||
|
||||
// Find the largest fftsize and sftsize of any plan
|
||||
for (size_t i = 0; i < anl_plans.size(); i++) {
|
||||
fftsize_max = std::max(fftsize_max, anl_plans[i]->fftsize);
|
||||
|
|
@ -1278,7 +1397,186 @@ struct analyzer: public refcounted {
|
|||
sftsize_max = std::max(sftsize_max, syn_plans[i]->sftsize_max);
|
||||
}
|
||||
|
||||
cmeta_any = anl_plans[0]->cmeta;
|
||||
// Lay out the coefficient structures according to the
|
||||
// synthesis plan if we have one, or the largest analysis
|
||||
// plan if not.
|
||||
std::vector<ref<plan>> *cmeta_source =
|
||||
params.synthesis ? &syn_plans : &anl_plans;
|
||||
ref<plan> &largest_plan(((*cmeta_source)[cmeta_source->size() - 1]));
|
||||
cmeta_any = make_meta(filet_part(largest_plan->fftsize));
|
||||
}
|
||||
|
||||
void make_zones() {
|
||||
// Band number starting at 0 close to fs/2 and increasing
|
||||
// with decreasing frequency
|
||||
int tbno = 0;
|
||||
int oct = 0;
|
||||
int zno = 0;
|
||||
// Loop over the octaves, from high to low frequencies,
|
||||
// creating new zones where needed
|
||||
for (;;) {
|
||||
int max_bands_this_octave = (zno == 0) ?
|
||||
n_bands_top_octave : params.bands_per_octave;
|
||||
int bands_remaining = n_bandpass_bands_total - tbno;
|
||||
int bands_this_octave = std::min(max_bands_this_octave, bands_remaining);
|
||||
int bands_below = bands_remaining - bands_this_octave;
|
||||
bool dc_zone = (bands_below == 0);
|
||||
bool dc_adjacent_zone = (bands_below < (int)params.bands_per_octave);
|
||||
if (zno < 2 || dc_zone || dc_adjacent_zone ||
|
||||
params.bands_per_octave < 6)
|
||||
{
|
||||
make_zone(oct, zno, tbno, tbno + bands_this_octave,
|
||||
dc_zone, bands_below);
|
||||
zno++;
|
||||
}
|
||||
octaves.push_back(octave<T>());
|
||||
octaves.back().z = zones[zno - 1].get();
|
||||
oct++;
|
||||
tbno += bands_this_octave;
|
||||
if (dc_zone)
|
||||
break;
|
||||
}
|
||||
assert(octaves.size() == n_octaves);
|
||||
}
|
||||
|
||||
// Create a zone consisting of the bandpass bands band0
|
||||
// (inclusive) to band1 (exclusive), using the usual gbno
|
||||
// numbering going from high to low frequencies, and
|
||||
// possible a lowpass band band1.
|
||||
|
||||
void make_zone(int oct, unsigned int zno,
|
||||
int band0, int band1,
|
||||
bool dc_zone, int bands_below)
|
||||
{
|
||||
assert(zones.size() == zno);
|
||||
zone<T> *z = new zone<T>();
|
||||
z->zno = zno;
|
||||
zones.push_back(ref<zone<T>>(z));
|
||||
|
||||
pod_vector<T> power;
|
||||
|
||||
// Create the real (non-mock) bands, from low to high
|
||||
// frequency.
|
||||
if (dc_zone)
|
||||
// This zone has a lowpass band
|
||||
push(z->bandparams, make_band(oct, band1, true, false));
|
||||
// The actual (non-mock) bandpass bands of this zone
|
||||
for (int i = band1 - 1; i >= band0; i--)
|
||||
push(z->bandparams, make_band(oct, i, false, false));
|
||||
|
||||
if (! dc_zone) {
|
||||
// There are other zones below this; add mock bands
|
||||
// to simulate them for purposes of calculating the dual.
|
||||
|
||||
// Identify the lowest frequency of interest in the zone
|
||||
assert(z->bandparams.size() >= 1);
|
||||
band_params<T> *low_band = z->bandparams[0].get();
|
||||
double zone_bottom_ff = low_band->ff - low_band->ff_support;
|
||||
|
||||
int i = band1;
|
||||
for (; i < band1 + bands_below; i++) {
|
||||
band_params<T> *mock_band = make_band(oct, i, false, true);
|
||||
push(z->mock_bandparams, mock_band);
|
||||
// There's no point in creating further mock bands
|
||||
// once they no longer overlap with the current zone.
|
||||
// The condition used here may cause the creation of
|
||||
// one more mock band than is actually needed, as it
|
||||
// is easier to create the band first and check for
|
||||
// overlap later than the other way round.
|
||||
if (mock_band->ff + mock_band->ff_support < zone_bottom_ff) {
|
||||
i++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
// Create a mock lowpass band. This may correspond to the
|
||||
// actual lowpass band, or if the loop above exited early,
|
||||
// it may merely be be a placeholder that serves no real
|
||||
// purpose other than making the power vector look better.
|
||||
push(z->mock_bandparams, make_band(oct, i, true, true));
|
||||
}
|
||||
|
||||
// If there are other zones above this, add mock bands
|
||||
// to simulate them for purposes of calculating the dual,
|
||||
// but only up to the Nyquist frequency of the current
|
||||
// octave.
|
||||
int nyquist_band = oct * (int)params.bands_per_octave;
|
||||
for (int i = band0 - 1; i >= nyquist_band; i--)
|
||||
push(z->mock_bandparams, make_band(oct, i, false, true));
|
||||
|
||||
z->n_bands = (unsigned int)z->bandparams.size();
|
||||
|
||||
// Find the largest coefficient step in the zone, as this will
|
||||
// determine the necessary alignment of signal slices in time,
|
||||
// but make it at least two (corresponding max_step_log2 = 1)
|
||||
// because the downsampling code requires alignement to even
|
||||
// indices.
|
||||
unsigned int m = 1;
|
||||
for (unsigned int obno = 0; obno < z->bandparams.size(); obno++) {
|
||||
m = std::max(m, z->bandparams[obno]->step_log2);
|
||||
}
|
||||
z->max_step_log2 = m;
|
||||
|
||||
max_step_log2 = std::max(max_step_log2, m);
|
||||
}
|
||||
|
||||
// Calculate band parameters for a single band.
|
||||
//
|
||||
// If dc is true, this is the DC band, and gbno indicates
|
||||
// the cutoff frequency; it is one more than the gbno of
|
||||
// the lowest-frequency bandpass band.
|
||||
|
||||
band_params<T> *
|
||||
make_band(int oct, double gbno, bool dc, bool mock) {
|
||||
band_params<T> *bp = new band_params<T>;
|
||||
if (dc)
|
||||
// Make the actual DC band cutoff frequency a bit higher,
|
||||
// by an empirically chosen fraction of a band, to reduce
|
||||
// power fluctuations.
|
||||
gbno -= 0.8750526596806952;
|
||||
|
||||
// For bandpass bands, the center frequency, or for the
|
||||
// lowpass band, the lowpass cutoff frequency, as a
|
||||
// fractional frequency, in terms of the octave's sample
|
||||
// rate.
|
||||
double ff = ldexp(bandpass_band_ff(gbno), oct);
|
||||
|
||||
// Standard deviation of the bandpass Gaussian, as a
|
||||
// fractional frequency
|
||||
double ffsd = ff_sd(ff);
|
||||
double time_sd = sd_f2t(ffsd);
|
||||
double time_support =
|
||||
gaussian_support(time_sd, params.max_error);
|
||||
|
||||
// The support of the Gaussian, i.e., the smallest standard
|
||||
// deviation at which it can be truncated on each side
|
||||
// without the error exceeding our part of the error budget,
|
||||
// which is some fraction of params.max_error. Note
|
||||
// that this is one-sided; the full width of the support
|
||||
// is 2 * ff_support.
|
||||
double bp_ff_support = gaussian_support(ffsd, params.max_error * 0.5);
|
||||
// Additional support for the flat portion of the DC band lowpass
|
||||
double dc_support = dc ? ff : 0;
|
||||
// Total frequency-domain support for this band, one-sided
|
||||
double band_support = bp_ff_support + dc_support;
|
||||
// Total support needed for this band, two-sided
|
||||
double band_2support = band_support * 2;
|
||||
|
||||
// Determine the downsampling factor for this band.
|
||||
int exp = 0;
|
||||
while (band_2support <= 0.5) {
|
||||
band_2support *= 2;
|
||||
exp++;
|
||||
}
|
||||
bp->step_log2 = exp;
|
||||
bp->step = 1U << bp->step_log2;
|
||||
|
||||
bp->dc = dc;
|
||||
bp->ff = ff;
|
||||
bp->ff_support = band_support;
|
||||
bp->time_support = time_support;
|
||||
bp->ffsd = ffsd;
|
||||
|
||||
return bp;
|
||||
}
|
||||
|
||||
// Given a fractional frequency, return the standard deviation
|
||||
|
|
@ -1289,7 +1587,8 @@ struct analyzer: public refcounted {
|
|||
// of the time-domain window in samples.
|
||||
//
|
||||
// ff_sd = 1.0 / (tau * t_sd)
|
||||
// per http://users.ece.gatech.edu/mrichard/Gaussian%20FT%20and%20random%20process.pdf
|
||||
// per http://users.ece.gatech.edu/mrichard/
|
||||
// Gaussian%20FT%20and%20random%20process.pdf
|
||||
// and python test program gaussian-overlap.py
|
||||
// => (tau * t_sd) * ff_sd = 1.0
|
||||
// => t_sd = 1.0 / (tau * f_sd)
|
||||
|
|
@ -1301,13 +1600,15 @@ struct analyzer: public refcounted {
|
|||
// the largest distance in time between a signal sample and a
|
||||
// coefficient affected by that sample.
|
||||
double analysis_support() const {
|
||||
return analysis_support(n_bands_total - 2); // Last before DC
|
||||
// The lowpass filter is the steepest one
|
||||
return analysis_support(band_lowpass());
|
||||
}
|
||||
|
||||
// Find the time support of the analysis filter for bandpass band
|
||||
// number gbno
|
||||
double analysis_support(double gbno) const {
|
||||
return gaussian_support(time_sd(bandpass_band_ff(gbno)), params.max_error);
|
||||
return gaussian_support(time_sd(bandpass_band_ff(gbno)),
|
||||
params.max_error);
|
||||
}
|
||||
|
||||
// Ditto for the resynthesis filters.
|
||||
|
|
@ -1323,15 +1624,14 @@ struct analyzer: public refcounted {
|
|||
// synthesis filters are wider than the analysis filters in
|
||||
// the time domain.
|
||||
double synthesis_support_multiplier() const {
|
||||
// Empirical formula
|
||||
return params.synthesis ? 1.8 / params.overlap : 1;
|
||||
}
|
||||
|
||||
// Return the fractional frequency of relative band number
|
||||
// "rbno" (relative to the sampling frequency of its octave).
|
||||
double rbno_ff(double rbno) const {
|
||||
double log2ff = band0_log2ff + rbno * band_spacing_log2;
|
||||
return exp2(log2ff);
|
||||
if (! params.synthesis)
|
||||
return 1.0;
|
||||
// At high ff_min, the top band may be the one with the widest support.
|
||||
if (params.ff_min > 0.25)
|
||||
return 4;
|
||||
if (params.bands_per_octave <= 4)
|
||||
return 2.5;
|
||||
return 2.3;
|
||||
}
|
||||
|
||||
// Get the center frequency of bandpass band "gbno", which
|
||||
|
|
@ -1349,23 +1649,16 @@ struct analyzer: public refcounted {
|
|||
return log2ff_bandpass_band(log2(ff));
|
||||
}
|
||||
|
||||
// Return the fractional frequency of bandpass band number "obno"
|
||||
// in octave "oct", scaling according to the octave
|
||||
double bandpass_band_ff(int oct, int obno) const {
|
||||
return bandpass_band_ff(bno_merge(oct, obno));
|
||||
}
|
||||
|
||||
const plan &choose_plan(const std::vector<ref<plan>> &plans,
|
||||
int64_t size) const
|
||||
int choose_plan(const std::vector<ref<plan>> &plans, int64_t size) const
|
||||
{
|
||||
unsigned int i = 0;
|
||||
while (i < plans.size() - 1 && plans[i]->filet_size < size)
|
||||
i++;
|
||||
return *plans[i];
|
||||
return i;
|
||||
}
|
||||
|
||||
void
|
||||
synthesize_one_slice(int oct, const plan &plan, const coefs<T> &msc,
|
||||
synthesize_one_slice(int oct, int pno, const coefs<T> &msc,
|
||||
const pod_vector<T> &downsampled,
|
||||
sample_index_t t0,
|
||||
T *signal_out,
|
||||
|
|
@ -1374,7 +1667,8 @@ struct analyzer: public refcounted {
|
|||
pod_vector<C> &buf3 // largest sftsize
|
||||
) const
|
||||
{
|
||||
zone<T> &z = *plan.octaves[oct].z;
|
||||
const plan &plan(*syn_plans[pno]);
|
||||
zone<T> &z = *octaves[oct].z;
|
||||
pod_vector<C> &signal(buf0);
|
||||
std::fill(signal.begin(), signal.end(), (T)0);
|
||||
|
||||
|
|
@ -1382,44 +1676,48 @@ struct analyzer: public refcounted {
|
|||
|
||||
for (unsigned int obno = 0; obno < z.bandparams.size(); obno++) {
|
||||
band_params<T> *bp = z.bandparams[obno].get();
|
||||
band_plan<T> *bpl = &bp->syn_plans[pno];
|
||||
|
||||
// log2 of the coefficient downsampling factor
|
||||
int coef_shift = bp->step_log2;
|
||||
coef_index_t ii = t0 >> coef_shift;
|
||||
|
||||
read(msc, bno_merge(oct, obno), ii, ii + bp->sftsize, coefbuf.data());
|
||||
read(msc, bno_merge(oct, obno), ii, ii + bpl->sftsize, coefbuf.data());
|
||||
C *indata = coefbuf.data();
|
||||
|
||||
pod_vector<C> &sdata(buf2);
|
||||
|
||||
T scale_factor = (T)1 / ((T)params.coef_scale * bp->sftsize);
|
||||
T scale_factor = (T)1 / ((T)params.coef_scale * bpl->sftsize);
|
||||
|
||||
// Apply phase correction, adjust for non-integer center
|
||||
// frequency, and apply scale factor. Note that phase
|
||||
// must be calculated in double precision.
|
||||
double ff = bp->center * plan.inv_fftsize_double;
|
||||
double arg = (params.phase == phase_convention::absolute) ?
|
||||
|
||||
// We can't use bp->ff here because in the case of the
|
||||
// lowpass band, it's the cutoff rather than the center.
|
||||
double ff = bpl->center * plan.inv_fftsize_double;
|
||||
double arg = (params.phase == coef_phase::global) ?
|
||||
tau * t0 * ff : 0;
|
||||
C phase_times_scale = C(cos(arg), sin(arg)) * scale_factor;
|
||||
elementwise_product_times_scalar(sdata.data(), indata,
|
||||
bp->shift_kernel_conj.data(),
|
||||
phase_times_scale, bp->sftsize);
|
||||
bpl->shift_kernel_conj.data(),
|
||||
phase_times_scale, bpl->sftsize);
|
||||
|
||||
// Switch to frequency domain
|
||||
bp->sft->transform(sdata.data());
|
||||
bpl->sft->transform(sdata.data());
|
||||
|
||||
// Multiply signal spectrum by frequency-domain dual window,
|
||||
// accumulating result in signal.
|
||||
|
||||
for (unsigned int i = 0; i < bp->sftsize; i++) {
|
||||
int iii = (bp->fq_offset_int + i) & (plan.fftsize - 1);
|
||||
for (unsigned int i = 0; i < bpl->sftsize; i++) {
|
||||
int iii = (bpl->fq_offset_int + i) & (plan.fftsize - 1);
|
||||
// Note the ifftshift of the input index, as f=0
|
||||
// appears in the middle of the window
|
||||
C v = sdata[i ^ (bp->sftsize >> 1)] * bp->dual_kernel[i];
|
||||
C v = sdata[i ^ (bpl->sftsize >> 1)] * bpl->dual_kernel[i];
|
||||
// Frequency symmetry
|
||||
signal[iii] += v;
|
||||
if (! bp->dc)
|
||||
signal[(plan.fftsize - iii) & (plan.fftsize - 1)] += conj(v);
|
||||
if (params.lowpass_version == 2 || ! bp->dc)
|
||||
signal[-iii & (plan.fftsize - 1)] += conj(v);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -1440,13 +1738,14 @@ struct analyzer: public refcounted {
|
|||
[i ^ (plan.dsparams.sftsize >> 1)];
|
||||
}
|
||||
|
||||
// This implicitly zero pads the spectrum, by
|
||||
// not adding anything to the middle part
|
||||
// The splitting of the nyquist band is per
|
||||
// http://dsp.stackexchange.com/questions/14919/upsample-data-using-ffts-how-is-this-exactly-done
|
||||
// but should not really matter because
|
||||
// there should be no energy there to speak of
|
||||
// thanks to the windowing above.
|
||||
// This implicitly zero pads the spectrum, by not adding
|
||||
// anything to the middle part. The splitting of the
|
||||
// Nyquist band is per http://dsp.stackexchange.com/
|
||||
// questions/14919/upsample-data-using-ffts-how-is-this-
|
||||
// exactly-done but should not really matter because there
|
||||
// should be no energy there to speak of thanks to the
|
||||
// windowing above.
|
||||
|
||||
assert(plan.dsparams.sftsize == plan.fftsize / 2);
|
||||
unsigned int i;
|
||||
for (i = 0; i < plan.dsparams.sftsize / 2; i++)
|
||||
|
|
@ -1482,9 +1781,12 @@ private:
|
|||
void
|
||||
analyze_sliced(buffers<T> &buf, int oct, const T *real_signal,
|
||||
sample_index_t t0, sample_index_t t1,
|
||||
double included_ds_support,
|
||||
coefs<T> &msc) const
|
||||
{
|
||||
const plan &plan(choose_plan(anl_plans, t1 - t0));
|
||||
assert(t1 >= t0);
|
||||
int pno = choose_plan(anl_plans, t1 - t0);
|
||||
const plan &plan(*anl_plans[pno]);
|
||||
|
||||
// Even though we don't align the FFTs to full filet-size
|
||||
// slices in this code path, we still need to align them to
|
||||
|
|
@ -1493,7 +1795,7 @@ private:
|
|||
// them to the largest coefficient time step of the octave.
|
||||
// slice_t0 is the sample time of the first sample in the
|
||||
// filet (not the FFT as a whole).
|
||||
zone<T> &z = *plan.octaves[oct].z;
|
||||
zone<T> &z = *octaves[oct].z;
|
||||
sample_index_t slice_t0 = t0 & ~((1 << z.max_step_log2) - 1);
|
||||
//printf("slice t0 = %d\n", (int)slice_t0);
|
||||
|
||||
|
|
@ -1517,7 +1819,7 @@ private:
|
|||
// nonzero data. Calculate adjusted bounds to use in the
|
||||
// recursive analysis so that we don't needlessly analyze
|
||||
// zeroes.
|
||||
int ds_support = ds_time_support;
|
||||
int ds_support = (int)ceil(ds_time_support - included_ds_support);
|
||||
sample_index_t dst0a = std::max(dst0, (t0 >> 1) - ds_support);
|
||||
sample_index_t dst1a = std::min(dst1, (t1 >> 1) + 1 + ds_support);
|
||||
|
||||
|
|
@ -1566,6 +1868,7 @@ private:
|
|||
|
||||
for (unsigned int obno = 0; obno < z.bandparams.size(); obno++) {
|
||||
band_params<T> *bp = z.bandparams[obno].get();
|
||||
band_plan<T> *bpl = &bp->anl_plans[pno];
|
||||
C *sdata = tmp.data();
|
||||
|
||||
// Multiply a slice of the spectrum by the frequency-
|
||||
|
|
@ -1579,21 +1882,21 @@ private:
|
|||
// checking in the inner loop, use a separate slow path
|
||||
// for the rare cases where wrapping happens.
|
||||
|
||||
int start_index = bp->fq_offset_int;
|
||||
int end_index = bp->fq_offset_int + bp->sftsize;
|
||||
int start_index = bpl->fq_offset_int;
|
||||
int end_index = bpl->fq_offset_int + bpl->sftsize;
|
||||
if (start_index >= 0 && end_index < (int)((plan.fftsize >> 1) + 1)) {
|
||||
// Fast path: the slice lies entirely within the
|
||||
// positive-frequency half of the spectrum (including
|
||||
// DC and Nyquist).
|
||||
elementwise_product(sdata,
|
||||
spectrum.data() + start_index,
|
||||
bp->kernel.data(),
|
||||
bp->sftsize);
|
||||
bpl->kernel.data(),
|
||||
bpl->sftsize);
|
||||
} else {
|
||||
// Slow path
|
||||
for (size_t i = 0; i < bp->sftsize; i++)
|
||||
for (size_t i = 0; i < bpl->sftsize; i++)
|
||||
sdata[i] = get_real_spectrum_coef(spectrum.data(),
|
||||
(int)(start_index + i), plan.fftsize) * bp->kernel[i];
|
||||
(int)(start_index + i), plan.fftsize) * bpl->kernel[i];
|
||||
}
|
||||
// The band center frequency is at the center of the
|
||||
// spectrum slice and at the center of the window, so
|
||||
|
|
@ -1609,22 +1912,22 @@ private:
|
|||
|
||||
// Switch to time domain
|
||||
auto band(buf.template get<C>(3, plan.sftsize_max));
|
||||
bp->sft->itransform(sdata, band.data());
|
||||
bpl->sft->itransform(sdata, band.data());
|
||||
|
||||
// Apply ifftshift, adjust for non-integer center
|
||||
// frequency, correct phase, scale amplitude, and add
|
||||
// to the output coefficients.
|
||||
double ff = bp->center * plan.inv_fftsize_double;
|
||||
double ff = bpl->center * plan.inv_fftsize_double;
|
||||
double arg;
|
||||
if (params.phase == phase_convention::absolute)
|
||||
if (params.phase == coef_phase::global)
|
||||
arg = -tau * (slice_t0 - plan.fat_size) * ff;
|
||||
else
|
||||
arg = 0;
|
||||
C phase_times_scale = C(cos(arg), sin(arg)) * scale_factor;
|
||||
elementwise_product_times_scalar(coefbuf.data(), band.data(),
|
||||
bp->shift_kernel.data(),
|
||||
bpl->shift_kernel.data(),
|
||||
phase_times_scale,
|
||||
bp->sftsize);
|
||||
bpl->sftsize);
|
||||
// log2 of the coefficient downsampling factor
|
||||
int coef_shift = bp->step_log2;
|
||||
assert(((slice_t0 - (int)plan.fat_size) & ((1 << coef_shift) - 1)) == 0);
|
||||
|
|
@ -1635,9 +1938,10 @@ private:
|
|||
// t0..t1 + the actual support of the filter for this band.
|
||||
// There's no point adding the zeros to the coefficients,
|
||||
// so trim.
|
||||
coef_index_t ii0 = std::max(ii, (t0 - (int)bp->time_support) >> coef_shift);
|
||||
coef_index_t ii1 = std::min(ii + bp->sftsize,
|
||||
((t1 + (int)bp->time_support) >> coef_shift) + 1);
|
||||
int support = (int)ceil(bp->time_support);
|
||||
coef_index_t ii0 = std::max(ii, (t0 - support) >> coef_shift);
|
||||
coef_index_t ii1 = std::min(ii + bpl->sftsize,
|
||||
((t1 + support) >> coef_shift) + 1);
|
||||
add(msc, bno_merge(oct, obno), ii0, ii1, coefbuf.data() + (ii0 - ii));
|
||||
}
|
||||
|
||||
|
|
@ -1692,7 +1996,7 @@ private:
|
|||
// Recurse
|
||||
if (oct + 1 < (int)n_octaves)
|
||||
analyze_sliced(buf, oct + 1, downsampled.data() + (dst0a - dst0),
|
||||
dst0a, dst1a, msc);
|
||||
dst0a, dst1a, ds_time_support / 2, msc);
|
||||
}
|
||||
|
||||
// Resynthesize audio from the coefficients in "msc". The audio will
|
||||
|
|
@ -1705,7 +2009,8 @@ private:
|
|||
sample_index_t t0, sample_index_t t1,
|
||||
T *real_signal) const
|
||||
{
|
||||
const plan &plan(choose_plan(syn_plans, t1 - t0));
|
||||
int pno = choose_plan(syn_plans, t1 - t0);
|
||||
const plan &plan(*syn_plans[pno]);
|
||||
|
||||
// XXX clean up - no need to pass support arg
|
||||
slice_index_t si0 = plan.affected_slice_b(t0, plan.oct_support);
|
||||
|
|
@ -1750,7 +2055,7 @@ private:
|
|||
downsampled.begin());
|
||||
}
|
||||
|
||||
synthesize_one_slice(oct, plan, msc, downsampled, slice_t0,
|
||||
synthesize_one_slice(oct, pno, msc, downsampled, slice_t0,
|
||||
signal_slice.data(), buf0, buf2, buf3);
|
||||
|
||||
// Copy overlapping part
|
||||
|
|
@ -1778,7 +2083,7 @@ public:
|
|||
assert(msc.octaves.size() == n_octaves);
|
||||
(void)n_threads;
|
||||
buffers<T> buf(fftsize_max, sftsize_max);
|
||||
analyze_sliced(buf, 0, real_signal, t0, t1, msc);
|
||||
analyze_sliced(buf, 0, real_signal, t0, t1, 0, msc);
|
||||
}
|
||||
|
||||
// The main synthesis entry point
|
||||
|
|
@ -1898,10 +2203,10 @@ public:
|
|||
double band_spacing_log2;
|
||||
double band_spacing;
|
||||
double tuning_log2ff;
|
||||
double band0_log2ff;
|
||||
linear_transform log2ff_bandpass_band;
|
||||
linear_transform bandpass_band_log2ff;
|
||||
affine_transform log2ff_bandpass_band;
|
||||
affine_transform bandpass_band_log2ff;
|
||||
unsigned int n_bandpass_bands_total;
|
||||
unsigned int n_bands_top_octave;
|
||||
|
||||
struct plan: public refcounted {
|
||||
plan(const plan &) = delete;
|
||||
|
|
@ -1909,7 +2214,6 @@ public:
|
|||
ok(false),
|
||||
synthesis(synthesis_),
|
||||
fftsize(fftsize_),
|
||||
max_step_log2(0),
|
||||
oct_support(support_),
|
||||
sftsize_max(0)
|
||||
{
|
||||
|
|
@ -1923,37 +2227,6 @@ public:
|
|||
#else
|
||||
ft = pool<fft<C *>, int>::shared.get(fftsize);
|
||||
#endif
|
||||
|
||||
// Band number starting at 0 close to fs/2 and increasing
|
||||
// with decreasing frequency
|
||||
int tbno = 0;
|
||||
// Band number of the fs/8 band in the current octave
|
||||
int base = 0;
|
||||
int zno = 0;
|
||||
// Loop over the octaves, from high to low frequencies,
|
||||
// creating new zones where needed
|
||||
for (;;) {
|
||||
int max_bands_this_octave = (zno == 0) ?
|
||||
anl->params.bands_per_octave * 2 : anl->params.bands_per_octave;
|
||||
base += max_bands_this_octave;
|
||||
int bands_remaining = anl->n_bandpass_bands_total - tbno;
|
||||
int bands_this_octave = std::min(max_bands_this_octave, bands_remaining);
|
||||
int bands_below = bands_remaining - bands_this_octave;
|
||||
bool dc_zone = (bands_below == 0);
|
||||
bool dc_adjacent_zone = (bands_below < (int)anl->params.bands_per_octave);
|
||||
if (zno < 2 || dc_zone || dc_adjacent_zone) {
|
||||
make_zone(anl, zno, base - (tbno + bands_this_octave),
|
||||
base - tbno, dc_zone, bands_below);
|
||||
zno++;
|
||||
}
|
||||
octaves.push_back(octave<T>());
|
||||
octaves.back().z = zones[zno - 1].get();
|
||||
tbno += bands_this_octave;
|
||||
if (dc_zone)
|
||||
break;
|
||||
}
|
||||
assert(octaves.size() == anl->n_octaves);
|
||||
|
||||
// Set up the downsampling parameters in dsparams.
|
||||
|
||||
// Downsampling is always by a factor of two.
|
||||
|
|
@ -1965,13 +2238,15 @@ public:
|
|||
if (synthesis)
|
||||
dsparams.dual_kernel.resize(dsparams.sftsize);
|
||||
|
||||
// Use the convolution of a rectangle and a gaussian.
|
||||
// Use the convolution of a rectangle and a Gaussian.
|
||||
// A piecewise function composed from two half-gaussians
|
||||
// joined by a horizontal y=1 segment is not quite smooth
|
||||
// enough. Put the passband in the middle.
|
||||
gaussian_windowed_lowpass(anl->ds_ff, anl->ds_ff_sd,
|
||||
dsparams.kernel.begin(),
|
||||
dsparams.kernel.end(), true);
|
||||
for (int i = 0; i < (int)dsparams.sftsize; i++)
|
||||
dsparams.kernel[i] =
|
||||
gaussian_windowed_lowpass_1(anl->ds_ff, anl->ds_ff_sd,
|
||||
((double)i / dsparams.sftsize) - 0.5);
|
||||
|
||||
if (synthesis) {
|
||||
// The dual_kernel field of the downsampling pseudo-band holds
|
||||
// the upsampling filter, identical to the downsampling filter
|
||||
|
|
@ -1992,9 +2267,9 @@ public:
|
|||
#endif
|
||||
// It may be possible to reduce the size of the fat from 1/4
|
||||
// of the fftsize, but we need to keep things aligned with the
|
||||
// coefficients
|
||||
// coefficients, and it needs to be even for downsampling.
|
||||
if (! synthesis) {
|
||||
unsigned int align = 1 << max_step_log2;
|
||||
unsigned int align = 1 << std::max(anl->max_step_log2, 2U);
|
||||
fat_size = (oct_support + (align - 1)) & ~(align - 1);
|
||||
// There must be room for at least one signal sample in each
|
||||
// half of the FFT; it can't be all fat
|
||||
|
|
@ -2005,268 +2280,10 @@ public:
|
|||
}
|
||||
filet_size = fftsize - 2 * fat_size;
|
||||
|
||||
// Initialize coefficient metadata
|
||||
cmeta = new coefs_meta;
|
||||
cmeta->n_octaves = anl->n_octaves;
|
||||
cmeta->n_bands_total = anl->n_bands_total;
|
||||
cmeta->bands_per_octave = anl->params.bands_per_octave;
|
||||
cmeta->slice_len = fftsize >> 1;
|
||||
cmeta->zones.resize(zones.size());
|
||||
for (unsigned int zi = 0; zi < zones.size(); zi++) {
|
||||
zone<T> *z = zones[zi].get();
|
||||
typename zone_coefs_meta::band_vector bv(z->bandparams.size());
|
||||
for (unsigned int i = 0; i < z->bandparams.size(); i++) {
|
||||
bv[i].slice_len = z->bandparams[i]->sftsize >> 1;
|
||||
bv[i].slice_len_log2 = z->bandparams[i]->sftsize_log2 - 1;
|
||||
bv[i].step_log2 = fftsize_log2 - z->bandparams[i]->sftsize_log2;
|
||||
}
|
||||
cmeta->zones[zi].init(bv);
|
||||
}
|
||||
cmeta->octaves.resize(octaves.size());
|
||||
tbno = 0;
|
||||
for (unsigned int i = 0; i < anl->n_octaves; i++) {
|
||||
cmeta->octaves[i].z = &cmeta->zones[this->octaves[i].z->zno];
|
||||
cmeta->octaves[i].n_bands_above = tbno;
|
||||
tbno += cmeta->octaves[i].z->bands.size();
|
||||
}
|
||||
// Constructor was successful
|
||||
ok = true;
|
||||
}
|
||||
|
||||
void make_zone(analyzer<T> *anl, unsigned int zno, int bb, int be, bool dc_zone, int bands_below) {
|
||||
assert(zones.size() == zno);
|
||||
zone<T> *z = new zone<T>();
|
||||
z->zno = zno;
|
||||
zones.push_back(ref<zone<T> >(z));
|
||||
|
||||
// The maximum number of mock bands to insert between the DC
|
||||
// band and the real bands.
|
||||
int n_mock_bands = 6;
|
||||
|
||||
if (dc_zone) {
|
||||
// This zone goes all the way to DC
|
||||
band_params<T> *dc_band = make_band(anl, bb - 1, true);
|
||||
push(z->bandparams, dc_band);
|
||||
} else {
|
||||
// There are other zones below this; add mock bands
|
||||
// to simulate them for purposes of calculating the dual
|
||||
if (n_mock_bands > bands_below)
|
||||
n_mock_bands = bands_below;
|
||||
// Mock DC band
|
||||
push(z->mock_bandparams,
|
||||
make_band(anl, bb - 1 - n_mock_bands, true));
|
||||
// Mock bandpass bands
|
||||
for (int i = bb - 1; i > bb - 1 - n_mock_bands; i--)
|
||||
push(z->mock_bandparams, make_band(anl, i, false));
|
||||
}
|
||||
|
||||
// The actual bandpass bands of this zone
|
||||
for (int i = bb; i < be; i++)
|
||||
push(z->bandparams, make_band(anl, i, false));
|
||||
|
||||
// If there are other zones above this, add mock bands
|
||||
// to simulate them for purposes of calculating the dual
|
||||
for (int i = be; i < (int)anl->params.bands_per_octave * 2; i++)
|
||||
push(z->mock_bandparams, make_band(anl, i, false));
|
||||
|
||||
z->n_bands = (unsigned int)z->bandparams.size();
|
||||
|
||||
// Find the largest coefficient step
|
||||
unsigned int m = 0;
|
||||
for (unsigned int obno = 0; obno < z->bandparams.size(); obno++) {
|
||||
m = std::max(m, z->bandparams[obno]->step_log2);
|
||||
}
|
||||
z->max_step_log2 = m;
|
||||
max_step_log2 = std::max(max_step_log2, m);
|
||||
|
||||
// Calculate complex exponentials for non-integer center
|
||||
// frequency adjustment and phase convention adjustment
|
||||
for (unsigned int obno = 0; obno < z->bandparams.size(); obno++) {
|
||||
band_params<T> *bp = z->bandparams[obno].get();
|
||||
for (unsigned int i = 0; i < bp->sftsize; i++) {
|
||||
double center =
|
||||
(anl->params.phase == phase_convention::absolute) ? bp->center : 0;
|
||||
double arg = tau * ((double)i / bp->sftsize) * -(center - bp->icenter);
|
||||
C t(cos(arg), sin(arg));
|
||||
// Apply ifftshift of spectrum in time domain
|
||||
if (i & 1)
|
||||
t = -t;
|
||||
bp->shift_kernel[i] = t;
|
||||
}
|
||||
}
|
||||
|
||||
if (synthesis) {
|
||||
// Accumulate window power for calculating dual
|
||||
z->power.resize(fftsize);
|
||||
std::fill(z->power.data(), z->power.data() + fftsize, (T)0);
|
||||
for (unsigned int obno = 0; obno < z->bandparams.size(); obno++)
|
||||
accumulate_power(z->bandparams[obno].get(), z->power.data());
|
||||
for (unsigned int obno = 0; obno < z->mock_bandparams.size(); obno++)
|
||||
accumulate_power(z->mock_bandparams[obno].get(), z->power.data());
|
||||
|
||||
// Calculate duals
|
||||
for (unsigned int obno = 0; obno < z->bandparams.size(); obno++) {
|
||||
band_params<T> *bp = z->bandparams[obno].get();
|
||||
for (unsigned int i = 0; i < bp->sftsize; i++) {
|
||||
// ii = large-FFT bin number
|
||||
int ii = i + bp->fq_offset_int;
|
||||
bp->dual_kernel[i] = bp->kernel[i] / z->power[ii & (fftsize - 1)];
|
||||
}
|
||||
for (unsigned int i = 0; i < bp->sftsize; i++) {
|
||||
C t = bp->shift_kernel[i];
|
||||
// Undo time-domain ifftshift
|
||||
if (i & 1)
|
||||
t = -t;
|
||||
bp->shift_kernel_conj[i] = conj(t);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Free the mock band parameters
|
||||
z->mock_bandparams.clear();
|
||||
z->power.clear();
|
||||
}
|
||||
|
||||
// Add the power of the kernel in "*bp" to "power"
|
||||
void
|
||||
accumulate_power(band_params<T> *bp, T *power) {
|
||||
for (unsigned int i = 0; i < bp->sftsize; i++) {
|
||||
// ii = large-FFT bin number
|
||||
unsigned int ii = (i + bp->fq_offset_int) & (fftsize - 1);
|
||||
assert(ii >= 0 && ii < fftsize);
|
||||
T y = bp->kernel[i];
|
||||
T p = y * y;
|
||||
power[ii] += p;
|
||||
if (! bp->dc) {
|
||||
unsigned int ni = fftsize - ii;
|
||||
ni &= fftsize - 1; // XXX is this correct?
|
||||
assert(ni < fftsize);
|
||||
power[ni] += p;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate band parameters for a single band.
|
||||
//
|
||||
// rbno is a "relative band number" indicating a frequency
|
||||
// within its octave: it is 0 for the fs/8 band and increases
|
||||
// with frequency. It is not always the same as the index
|
||||
// into its octaves' bandparams[] or the coefficient bands
|
||||
// (those would be called obno, not rbno).
|
||||
//
|
||||
// If dc is true, this is the DC band, and rbno indicates
|
||||
// the cutoff frequency; it is one less than the rbno of
|
||||
// the lowest-frequency bandpass band.
|
||||
|
||||
band_params<T> *
|
||||
make_band(analyzer *anl, double rbno, bool dc) {
|
||||
if (dc)
|
||||
// Make the actual DC band cutoff frequency a bit higher,
|
||||
// by an empirically chosen fraction of a band, to reduce
|
||||
// power fluctuations.
|
||||
rbno += 0.8750526596806952;
|
||||
|
||||
// For bandpass bands, the center frequency, or
|
||||
// for the DC band, the lowpass cutoff frequency,
|
||||
// as a fractional frequency.
|
||||
double ff = anl->rbno_ff(rbno);
|
||||
|
||||
// Standard deviation of the bandpass Gaussian,
|
||||
// as a fractional frequency
|
||||
double ffsd = anl->ff_sd(ff);
|
||||
// The support of the Gaussian, i.e., the smallest standard
|
||||
// deviation at which it can be truncated on each side
|
||||
// without the error exceeding our part of the error budget,
|
||||
// which is some fraction of params.max. Note
|
||||
// that this is one-sided; the full width of the support
|
||||
// is 2 * ff_support.
|
||||
double ff_support = gaussian_support(ffsd, anl->params.max_error * 0.5);
|
||||
// Additional support for the flat portion of the DC band lowpass
|
||||
double dc_support = dc ? ff : 0;
|
||||
// The support as the number of FFT frequency bands needed,
|
||||
// allowing for both sides of the Gaussian.
|
||||
int fq_2support = int(ceil((ff_support + dc_support) * 2 * fftsize));
|
||||
|
||||
unsigned int sftsize = next_power_of_two(fq_2support);
|
||||
|
||||
// If there is a very small number of bandpass bands, the
|
||||
// taper of the lowpass band may extend past the Nyquist
|
||||
// frequency, causing sftsize to become greater than fftsize.
|
||||
// This would mean sampling the lowpass band at a higher sample
|
||||
// rate than the original signal, which would be a bit silly,
|
||||
// and having sftsize > fftsize will break things by causing
|
||||
// shift amounts to go negative, so instead we force sftsize
|
||||
// = fftsize but remember the original sftsize. The lowpass
|
||||
// kernel is calculated using the original sftsize and allowed
|
||||
// to alias over the truncated sftsize.
|
||||
unsigned int orig_sftsize = sftsize;
|
||||
if (sftsize > fftsize) {
|
||||
assert(dc);
|
||||
sftsize = fftsize;
|
||||
}
|
||||
|
||||
band_params<T> *bp = new band_params<T>;
|
||||
bp->dc = dc;
|
||||
bp->sftsize = sftsize;
|
||||
bp->sftsize_log2 = whichp2(bp->sftsize);
|
||||
bp->step_log2 = fftsize_log2 - bp->sftsize_log2;
|
||||
bp->step = 1U << bp->step_log2;
|
||||
sftsize_max = std::max(sftsize_max, bp->sftsize);
|
||||
bp->sft = pool<fft<C *>, int>::shared.get(bp->sftsize);
|
||||
|
||||
bp->kernel.resize(bp->sftsize);
|
||||
bp->shift_kernel.resize(bp->sftsize);
|
||||
if (synthesis) {
|
||||
bp->dual_kernel.resize(bp->sftsize);
|
||||
bp->shift_kernel_conj.resize(bp->sftsize);
|
||||
}
|
||||
|
||||
if (dc)
|
||||
bp->center = 0;
|
||||
else
|
||||
bp->center = ff * fftsize;
|
||||
bp->icenter = (int)rint(bp->center);
|
||||
bp->ff = ff;
|
||||
bp->fq_offset_int = bp->icenter - (bp->sftsize >> 1);
|
||||
bp->ffsd = ffsd;
|
||||
|
||||
// Calculate frequency-domain window kernel
|
||||
|
||||
if (dc) {
|
||||
// The cutoff frequency is a fraction of the fftsize, but
|
||||
// gaussian_windowed_lowpass() designs the filter in terms
|
||||
// of the the sftsize, so we need to scale the frequencies
|
||||
// accordingly. The cast to double is necessary to get
|
||||
// a float division in pathological cases where orig_sftsize
|
||||
// is larger than fftsize.
|
||||
double scale = (double)fftsize / orig_sftsize;
|
||||
pod_vector<T> win(orig_sftsize);
|
||||
gaussian_windowed_lowpass(
|
||||
ff * scale,
|
||||
ffsd * scale,
|
||||
win.begin(), win.end());
|
||||
for (unsigned int i = 0; i < bp->kernel.size(); i++)
|
||||
bp->kernel[i] = 0;
|
||||
for (unsigned int i = 0; i < win.size(); i++)
|
||||
bp->kernel[i & (bp->kernel.size() - 1)] += win[i];
|
||||
fftshift(bp->kernel.begin(), bp->kernel.end());
|
||||
} else {
|
||||
for (unsigned int i = 0; i < bp->sftsize; i++) {
|
||||
// ii = large-FFT band number
|
||||
unsigned int ii = i + bp->fq_offset_int;
|
||||
// this_ff = fractional frequency of this kernel sample
|
||||
double this_ff = ii * inv_fftsize_double;
|
||||
T y = norm_gaussian(ffsd, this_ff - ff);
|
||||
bp->kernel[i] = y;
|
||||
}
|
||||
}
|
||||
|
||||
double time_sd = sd_f2t(ffsd);
|
||||
double time_support =
|
||||
gaussian_support(time_sd, anl->params.max_error);
|
||||
bp->time_support = ceil(time_support);
|
||||
return bp;
|
||||
}
|
||||
|
||||
// Index of first slice affected by sample at t0
|
||||
|
||||
// fft number i covers the sample range
|
||||
|
|
@ -2299,7 +2316,6 @@ public:
|
|||
unsigned int fftsize; // The size of the main FFT, a power of two.
|
||||
unsigned int fat_size;
|
||||
unsigned int filet_size;
|
||||
unsigned int max_step_log2;
|
||||
|
||||
// The width of the widest filter in the time domain, in
|
||||
// octave subsamples
|
||||
|
|
@ -2309,7 +2325,6 @@ public:
|
|||
T inv_fftsize_t; // 1.0f / fftsize (if using floats)
|
||||
|
||||
unsigned int sftsize_max; // The size of the largest band FFT, a power of two
|
||||
std::vector<ref<zone<T> > > zones;
|
||||
downsampling_params<T> dsparams;
|
||||
|
||||
// Fourier transform object for transforming a full slice
|
||||
|
|
@ -2318,11 +2333,207 @@ public:
|
|||
#else
|
||||
fft<C *> *ft;
|
||||
#endif
|
||||
std::vector<octave<T> > octaves; // Per-octave parameters
|
||||
|
||||
ref<coefs_meta> cmeta; // Coefficient metadata
|
||||
};
|
||||
|
||||
// Calculate per-plan, per-band coefficients for plan "pno",
|
||||
// a synthesis plan if "syn" is true, otherwise an analysis plan.
|
||||
|
||||
void make_band_plans(int pno, bool syn) {
|
||||
std::vector<ref<typename analyzer<T>::plan>> &plans
|
||||
(syn ? syn_plans : anl_plans);
|
||||
plan &plan(*plans[pno].get());
|
||||
|
||||
for (int zno = 0; zno < (int)zones.size(); zno++) {
|
||||
zone<T> *z = zones[zno].get();
|
||||
|
||||
make_band_plans_2(z->bandparams, pno, syn, false);
|
||||
make_band_plans_2(z->mock_bandparams, pno, syn, true);
|
||||
|
||||
if (plan.synthesis) {
|
||||
// Accumulate window power for calculating dual
|
||||
std::vector<T> power(plan.fftsize);
|
||||
// Real bands
|
||||
for (unsigned int i = 0; i < z->bandparams.size(); i++) {
|
||||
band_params<T> *bp = z->bandparams[i].get();
|
||||
band_plan<T> *bpl = syn ? &bp->syn_plans[pno] : &bp->anl_plans[pno];
|
||||
accumulate_power(plan, bp, bpl, power.data());
|
||||
}
|
||||
// Mock bands
|
||||
for (unsigned int i = 0; i < z->mock_bandparams.size(); i++) {
|
||||
band_params<T> *bp = z->mock_bandparams[i].get();
|
||||
band_plan<T> *bpl = syn ? &bp->syn_plans[pno] : &bp->anl_plans[pno];
|
||||
accumulate_power(plan, bp, bpl, power.data());
|
||||
}
|
||||
|
||||
// Calculate duals
|
||||
for (unsigned int obno = 0; obno < z->bandparams.size(); obno++) {
|
||||
band_params<T> *bp = z->bandparams[obno].get();
|
||||
band_plan<T> *bpl = syn ? &bp->syn_plans[pno] : &bp->anl_plans[pno];
|
||||
for (unsigned int i = 0; i < bpl->sftsize; i++) {
|
||||
// ii = large-FFT bin number
|
||||
int ii = i + bpl->fq_offset_int;
|
||||
bpl->dual_kernel[i] /= power[ii & (plan.fftsize - 1)];
|
||||
}
|
||||
// The analysis kernels are no longer needed
|
||||
bpl->kernel = std::vector<T>();
|
||||
bpl->shift_kernel = pod_vector<C>();
|
||||
}
|
||||
z->mock_bandparams.clear();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void make_band_plans_2(std::vector<ref<band_params<T>>> &bv, int pno,
|
||||
bool syn, bool mock)
|
||||
{
|
||||
std::vector<ref<typename analyzer<T>::plan>> &plans
|
||||
(syn ? syn_plans : anl_plans);
|
||||
plan &plan(*plans[pno].get());
|
||||
|
||||
for (unsigned int obno = 0; obno < bv.size(); obno++) {
|
||||
band_params<T> *bp = bv[obno].get();
|
||||
std::vector<band_plan<T>> *bplv = syn ? &bp->syn_plans : &bp->anl_plans;
|
||||
// XXX redundant resizes
|
||||
bplv->resize(plans.size());
|
||||
band_plan<T> *bpl = &(*bplv)[pno];
|
||||
|
||||
// Note that bp->step_log2 cannot be negative, meaning
|
||||
// that the bands can only be subsampled, not oversampled.
|
||||
unsigned int sftsize = plan.fftsize >> bp->step_log2;
|
||||
|
||||
// PFFFT has a minimum size
|
||||
sftsize = std::max(sftsize, (unsigned int)GABORATOR_MIN_FFT_SIZE);
|
||||
|
||||
bpl->sftsize = sftsize;
|
||||
bpl->sftsize_log2 = whichp2(bpl->sftsize);
|
||||
|
||||
if (! mock) {
|
||||
plan.sftsize_max = std::max(plan.sftsize_max, bpl->sftsize);
|
||||
bpl->sft = pool<fft<C *>, int>::shared.get(bpl->sftsize);
|
||||
}
|
||||
|
||||
bpl->kernel.resize(bpl->sftsize);
|
||||
bpl->shift_kernel.resize(bpl->sftsize);
|
||||
if (plan.synthesis) {
|
||||
bpl->dual_kernel.resize(bpl->sftsize);
|
||||
bpl->shift_kernel_conj.resize(bpl->sftsize);
|
||||
}
|
||||
|
||||
if (bp->dc)
|
||||
bpl->center = 0;
|
||||
else
|
||||
bpl->center = bp->ff * plan.fftsize;
|
||||
bpl->icenter = (int)rint(bpl->center);
|
||||
bpl->fq_offset_int = bpl->icenter - (bpl->sftsize >> 1);
|
||||
|
||||
// Calculate frequency-domain window kernel, possibly with
|
||||
// wrap-around
|
||||
for (unsigned int i = 0; i < bpl->sftsize; i++)
|
||||
bpl->kernel[i] = 0;
|
||||
// i loops over the kernel, with i=0 at the center.
|
||||
// The range is twice the support on each side so that
|
||||
// any excess space in the kernel due to rounding up
|
||||
// the size to a power of two is filled in with actual
|
||||
// Gaussian values rather than zeros.
|
||||
int fq_support = (int)ceil(bp->ff_support * plan.fftsize);
|
||||
for (int i = - 2 * fq_support; i < 2 * fq_support; i++) {
|
||||
// ii = large-FFT band number of this kernel sample
|
||||
int ii = i + bpl->fq_offset_int + (int)bpl->sftsize / 2;
|
||||
// this_ff = fractional frequency of this kernel sample
|
||||
double this_ff = ii * plan.inv_fftsize_double;
|
||||
// ki = kernel index
|
||||
int ki = ii - bpl->fq_offset_int;
|
||||
// When sftsize == fftsize, the support of the kernel can
|
||||
// exceed sftsize, and in this case, it should be allowed
|
||||
// to wrap so that it remains smooth. When sftsize < fftsize,
|
||||
// sftsize is large enough for the support and no wrapping
|
||||
// is needed or wanted.
|
||||
if (bpl->kernel.size() == plan.fftsize && !mock) {
|
||||
bpl->kernel[ki & (plan.fftsize - 1)] +=
|
||||
eval_kernel(¶ms, bp, this_ff);
|
||||
if (plan.synthesis)
|
||||
bpl->dual_kernel[ki & (plan.fftsize - 1)] +=
|
||||
eval_dual_kernel(¶ms, bp, this_ff);
|
||||
} else {
|
||||
if (ki >= 0 && ki < (int)bpl->kernel.size()) {
|
||||
bpl->kernel[ki] += eval_kernel(¶ms, bp, this_ff);
|
||||
if (plan.synthesis)
|
||||
bpl->dual_kernel[ki] = eval_dual_kernel(¶ms, bp, this_ff);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate complex exponentials for non-integer center
|
||||
// frequency adjustment and phase convention adjustment
|
||||
for (unsigned int obno = 0; obno < bv.size(); obno++) {
|
||||
band_params<T> *bp = bv[obno].get();
|
||||
band_plan<T> *bpl = syn ? &bp->syn_plans[pno] : &bp->anl_plans[pno];
|
||||
for (unsigned int i = 0; i < bpl->sftsize; i++) {
|
||||
double center =
|
||||
(params.phase == coef_phase::global) ? bpl->center : 0;
|
||||
double arg = tau * ((double)i / bpl->sftsize) * -(center - bpl->icenter);
|
||||
C t(cos(arg), sin(arg));
|
||||
// Apply ifftshift of spectrum in time domain
|
||||
bpl->shift_kernel[i] = (i & 1) ? -t : t;
|
||||
if (plan.synthesis)
|
||||
// Conjugate kernel does not have ifftshift
|
||||
bpl->shift_kernel_conj[i] = conj(t);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Add the power of the kernel in "*bp" to "power"
|
||||
void
|
||||
accumulate_power(plan &plan, band_params<T> *bp, band_plan<T> *bpl, T *power) {
|
||||
for (unsigned int i = 0; i < bpl->sftsize; i++) {
|
||||
// ii = large-FFT bin number
|
||||
unsigned int ii = i + bpl->fq_offset_int;
|
||||
ii &= plan.fftsize - 1;
|
||||
assert(ii >= 0 && ii < plan.fftsize);
|
||||
T p = bpl->kernel[i] * bpl->dual_kernel[i];
|
||||
power[ii] += p;
|
||||
if (params.lowpass_version == 2 || ! bp->dc) {
|
||||
unsigned int ni = -ii;
|
||||
ni &= plan.fftsize - 1;
|
||||
power[ni] += p;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Create coefficient metadata based on a slice length
|
||||
coefs_meta *make_meta(int slice_len) const {
|
||||
coefs_meta *cmeta = new coefs_meta;
|
||||
cmeta->n_octaves = n_octaves;
|
||||
cmeta->n_bands_total = n_bands_total;
|
||||
cmeta->bands_per_octave = params.bands_per_octave;
|
||||
cmeta->slice_len = slice_len;
|
||||
cmeta->zones.resize(zones.size());
|
||||
for (unsigned int zi = 0; zi < zones.size(); zi++) {
|
||||
zone<T> *z = zones[zi].get();
|
||||
typename zone_coefs_meta::band_vector bv(z->bandparams.size());
|
||||
for (unsigned int i = 0; i < z->bandparams.size(); i++) {
|
||||
unsigned int step_log2 = z->bandparams[i]->step_log2;
|
||||
bv[i].step_log2 = step_log2;
|
||||
bv[i].slice_len_log2 = whichp2(slice_len) - step_log2;
|
||||
bv[i].slice_len = 1 << bv[i].slice_len_log2;
|
||||
}
|
||||
cmeta->zones[zi].init(bv);
|
||||
}
|
||||
cmeta->octaves.resize(octaves.size());
|
||||
int tbno = 0;
|
||||
for (unsigned int i = 0; i < n_octaves; i++) {
|
||||
cmeta->octaves[i].z = &cmeta->zones[this->octaves[i].z->zno];
|
||||
cmeta->octaves[i].n_bands_above = tbno;
|
||||
tbno += cmeta->octaves[i].z->bands.size();
|
||||
}
|
||||
return cmeta;
|
||||
}
|
||||
|
||||
std::vector<octave<T>> octaves; // Per-octave parameters
|
||||
std::vector<ref<zone<T>>> zones;
|
||||
unsigned int max_step_log2;
|
||||
|
||||
std::vector<ref<plan>> anl_plans;
|
||||
std::vector<ref<plan>> syn_plans;
|
||||
|
||||
|
|
@ -2332,17 +2543,13 @@ public:
|
|||
double top_band_log2ff; // log2 of fractional frequency of the highest-frequency band
|
||||
int ffref_gbno; // Band number of the reference frequency
|
||||
|
||||
// Width of the downsampling filter passband in terms of the
|
||||
// downsampled sample rate (between 0.25 and 0.5)
|
||||
double ds_passband;
|
||||
double ds_ff; // Downsampling filter -6 dB transition frequency
|
||||
double ds_ff_sd; // Downsampling filter standard deviation
|
||||
double ds_time_support; // Downsampling filter time-domain kernel support, each side
|
||||
|
||||
// The largest time-domain support of any bandpass filter,
|
||||
// in octave subsamples. This is determined by the "band0"
|
||||
// filter (near fs/8).
|
||||
|
||||
double band0_time_analysis_support;
|
||||
double band0_time_synthesis_support;
|
||||
|
||||
unsigned int fftsize_max; // Largest FFT size of any plan
|
||||
unsigned int sftsize_max; // Largest SFT size of any plan
|
||||
|
||||
|
|
@ -2394,7 +2601,7 @@ void foreach_slice(unsigned int sh, coef_index_t i0, coef_index_t i1, F f) {
|
|||
// D is the dest object type
|
||||
// C is the coefficient type
|
||||
|
||||
template <class T, class D, class C = complex<T> >
|
||||
template <class T, class D, class C = complex<T>>
|
||||
struct row_foreach_slice {
|
||||
typedef C value_type;
|
||||
row_foreach_slice(const coefs<T, C> &msc,
|
||||
|
|
@ -2729,7 +2936,7 @@ void forget_before(const analyzer<T> &, coefs<T> &msc,
|
|||
sc.slices.erase_before(sli);
|
||||
if (clean_cut) {
|
||||
// Partially erase slice at boundary, if any
|
||||
const ref<oct_coefs<C> > *t = sc.slices.get(sli);
|
||||
const ref<oct_coefs<C>> *t = sc.slices.get(sli);
|
||||
if (! t)
|
||||
continue;
|
||||
if (! *t)
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// The Gaussian and related functions
|
||||
//
|
||||
// Copyright (C) 2015-2018 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2015-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -48,7 +48,7 @@ double gaussian_edge(double sd, double x) {
|
|||
return (erf(erf_arg) + 1) * 0.5;
|
||||
}
|
||||
|
||||
// Translate the time-domain standard deviation of a gaussian
|
||||
// Translate the time-domain standard deviation of a Gaussian
|
||||
// (in samples) into the corresponding frequency-domain standard
|
||||
// deviation (as a fractional frequency), or vice versa.
|
||||
|
||||
|
|
@ -60,21 +60,64 @@ static inline double sd_f2t(double ff_sd) {
|
|||
return sd_t2f(ff_sd);
|
||||
}
|
||||
|
||||
// Given a gaussian with standard deviation "sd" and a maximum error
|
||||
// Given a Gaussian with standard deviation "sd" and a maximum error
|
||||
// "max_error", calculate the support needed on each side to keep the
|
||||
// area below the curve within max_error of the exact value.
|
||||
|
||||
static inline double gaussian_support(double sd, double max_error) {
|
||||
static inline
|
||||
double gaussian_area_support(double sd, double max_error) {
|
||||
return sd * M_SQRT2 * erfc_inv(max_error);
|
||||
}
|
||||
|
||||
// Inverse of the above: given a support and maximum error, calculate
|
||||
// the standard deviation.
|
||||
|
||||
static inline double gaussian_support_inv(double support, double max_error) {
|
||||
static inline
|
||||
double gaussian_area_support_inv(double support, double max_error) {
|
||||
return support / (M_SQRT2 * erfc_inv(max_error));
|
||||
}
|
||||
|
||||
// Given a gaussian with standard deviation "sd" and a maximum error
|
||||
// "max_error", calculate the support needed on each side for the
|
||||
// value to fall to a factor of "max_error" of the peak.
|
||||
|
||||
static inline
|
||||
double gaussian_value_support(double sd, double max_error) {
|
||||
return sd * M_SQRT2 * sqrt(-log(max_error));
|
||||
}
|
||||
|
||||
// Inverse of the above: given a support and maximum error, calculate
|
||||
// the standard deviation.
|
||||
|
||||
static inline
|
||||
double gaussian_value_support_inv(double support, double max_error) {
|
||||
return support / (M_SQRT2 * sqrt(-log(max_error)));
|
||||
}
|
||||
|
||||
// Choose which criterion to use
|
||||
|
||||
#if 1
|
||||
static inline
|
||||
double gaussian_support(double support, double max_error) {
|
||||
return gaussian_area_support(support, max_error);
|
||||
};
|
||||
|
||||
static inline
|
||||
double gaussian_support_inv(double support, double max_error) {
|
||||
return gaussian_area_support_inv(support, max_error);
|
||||
};
|
||||
#else
|
||||
static inline
|
||||
double gaussian_support(double support, double max_error) {
|
||||
return gaussian_value_support(support, max_error);
|
||||
};
|
||||
|
||||
static inline
|
||||
double gaussian_support_inv(double support, double max_error) {
|
||||
return gaussian_value_support_inv(support, max_error);
|
||||
};
|
||||
#endif
|
||||
|
||||
} // namespace
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -1,35 +0,0 @@
|
|||
//
|
||||
// A class for linear transforms of the form ax + b
|
||||
//
|
||||
// Copyright (C) 2020 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
||||
#ifndef _GABORATOR_LINEAR_TRANSFORM_H
|
||||
#define _GABORATOR_LINEAR_TRANSFORM_H
|
||||
|
||||
namespace gaborator {
|
||||
|
||||
struct linear_transform {
|
||||
linear_transform(): a(0), b(0) { }
|
||||
linear_transform(double a_, double b_): a(a_), b(b_) { }
|
||||
linear_transform(const linear_transform &rhs): a(rhs.a), b(rhs.b) { }
|
||||
double operator()(double x) const { return a * x + b; }
|
||||
linear_transform inverse() const {
|
||||
return linear_transform(1.0 / a, -b / a);
|
||||
}
|
||||
static linear_transform identity() { return linear_transform(1, 0); }
|
||||
double a, b;
|
||||
};
|
||||
|
||||
// Composition
|
||||
|
||||
static inline linear_transform
|
||||
operator *(const linear_transform &a, const linear_transform &b) {
|
||||
return linear_transform(a.a * b.a, a.a * b.b + a.b);
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
#endif
|
||||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// A vector class without default-initialization, for "plain old data"
|
||||
//
|
||||
// Copyright (C) 2016-2019 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2016-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -84,7 +84,7 @@ struct pod_vector {
|
|||
if (&a == this)
|
||||
return *this;
|
||||
_free();
|
||||
b = new T[a.size()];
|
||||
b = static_cast<T *>(::operator new(a.size() * sizeof(T)));
|
||||
e = b + a.size();
|
||||
std::copy(a.b, a.e, b);
|
||||
//if (size()) fprintf(stderr, "pod_vector = %d\n", (int)size());
|
||||
|
|
@ -94,6 +94,7 @@ struct pod_vector {
|
|||
pod_vector &operator=(pod_vector &&x) noexcept {
|
||||
if (&x == this)
|
||||
return *this;
|
||||
_free();
|
||||
b = x.b;
|
||||
e = x.e;
|
||||
x.b = x.e = 0;
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// Rendering of spectrogram images
|
||||
//
|
||||
// Copyright (C) 2015-2020 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2015-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -36,7 +36,6 @@ unsigned int float2pixel_8bit(T val) {
|
|||
return (unsigned int)(sqrtf(val) * 255.0f);
|
||||
}
|
||||
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
// Power-of-two rendering
|
||||
|
|
@ -116,7 +115,7 @@ struct abs_writer_dest {
|
|||
// for vectorization.
|
||||
|
||||
template <class T, class C, class OI>
|
||||
struct abs_row_source<T, C, OI, struct complex_abs_fob<T> > {
|
||||
struct abs_row_source<T, C, OI, struct complex_abs_fob<T>> {
|
||||
// Note unused last arg
|
||||
abs_row_source(const coefs<T, C> &coefs_,
|
||||
int oct_, unsigned int obno_,
|
||||
|
|
@ -139,7 +138,7 @@ template <class OI, class T, class C, class NORMF, class RESAMPLER>
|
|||
OI render_line(const analyzer<T> &anl,
|
||||
const coefs<T, C> &msc,
|
||||
int gbno,
|
||||
linear_transform xf,
|
||||
affine_transform xf,
|
||||
sample_index_t i0, sample_index_t i1,
|
||||
OI output,
|
||||
NORMF normf)
|
||||
|
|
@ -160,8 +159,7 @@ OI render_line(const analyzer<T> &anl,
|
|||
// Scale by the downsampling factor of the band
|
||||
int scale_exp = band_scale_exp(*msc.octaves[oct].meta, oct, obno);
|
||||
RESAMPLER x_resampler(zoom_p2(xf, -scale_exp));
|
||||
output = x_resampler.resample(abs_rowsource,
|
||||
i0, i1, output);
|
||||
output = x_resampler(abs_rowsource, i0, i1, output);
|
||||
return output;
|
||||
}
|
||||
|
||||
|
|
@ -173,10 +171,11 @@ OI render_line(const analyzer<T> &anl,
|
|||
template <class OI, class T, class C, class NORMF, class RESAMPLER>
|
||||
OI render_noyscale(const analyzer<T> &anl,
|
||||
const coefs<T, C> &msc,
|
||||
linear_transform x_xf,
|
||||
affine_transform x_xf,
|
||||
int64_t xi0, int64_t xi1,
|
||||
int64_t yi0, int64_t yi1,
|
||||
OI output,
|
||||
int64_t output_stride,
|
||||
NORMF normf)
|
||||
{
|
||||
assert(xi1 >= xi0);
|
||||
|
|
@ -185,18 +184,11 @@ OI render_noyscale(const analyzer<T> &anl,
|
|||
int gbno0 = (int)yi0;
|
||||
int gbno1 = (int)yi1;
|
||||
for (int gbno = gbno0; gbno < gbno1; gbno++) {
|
||||
int oct;
|
||||
unsigned int obno; // Band number within octave
|
||||
bool clip = ! bno_split(*msc.meta, gbno, oct, obno, false);
|
||||
if (clip) {
|
||||
for (int x = 0; x < w; x++)
|
||||
*output++ = (T)0;
|
||||
} else {
|
||||
output = render_line<OI, T, C, NORMF, RESAMPLER>
|
||||
(anl, msc, gbno, x_xf,
|
||||
xi0, xi1,
|
||||
output, normf);
|
||||
}
|
||||
render_line<OI, T, C, NORMF, RESAMPLER>
|
||||
(anl, msc, gbno, x_xf,
|
||||
xi0, xi1,
|
||||
output, normf);
|
||||
output += output_stride;
|
||||
}
|
||||
return output;
|
||||
}
|
||||
|
|
@ -294,7 +286,7 @@ template <class OI, class T, class C, class NORMF = complex_abs_fob<T>,
|
|||
void render_incremental(
|
||||
const analyzer<T> &anl,
|
||||
const coefs<T, C> &msc,
|
||||
linear_transform x_xf, linear_transform y_xf,
|
||||
affine_transform x_xf, affine_transform y_xf,
|
||||
int64_t xi0, int64_t xi1,
|
||||
int64_t yi0, int64_t yi1,
|
||||
int64_t inc_i0, int64_t inc_i1,
|
||||
|
|
@ -457,25 +449,33 @@ void render_incremental(
|
|||
return;
|
||||
}
|
||||
|
||||
// Construct a temporary image resampled in the
|
||||
// X direction, but not yet in the Y direction. Include
|
||||
// extra scanlines at the top and bottom for interpolation.
|
||||
size_t n_pixels = (size_t)((ysi1 - ysi0) * (xi1 - xi0));
|
||||
pod_vector<RST> render_data(n_pixels);
|
||||
if (y_xf.a == 1 && y_xf.b == 0) {
|
||||
// No Y resampling needed, render data resampled in the X
|
||||
// direction only.
|
||||
render_noyscale<OI, T, C, NORMF, RESAMPLER>
|
||||
(anl, msc, x_xf, xi0, xi1, yi0, yi1,
|
||||
output, output_stride, normf);
|
||||
} else {
|
||||
// Construct a temporary image resampled in the
|
||||
// X direction, but not yet in the Y direction. Include
|
||||
// extra scanlines at the top and bottom for interpolation.
|
||||
size_t n_pixels = (size_t)((ysi1 - ysi0) * (xi1 - xi0));
|
||||
pod_vector<RST> render_data(n_pixels);
|
||||
|
||||
// Render data resampled in the X direction
|
||||
RST *p = render_data.data();
|
||||
render_noyscale<OI, T, C, NORMF, RESAMPLER>
|
||||
(anl, msc, x_xf, xi0, xi1,
|
||||
ysi0, ysi1, p, normf);
|
||||
// Render data resampled in the X direction
|
||||
RST *p = render_data.data();
|
||||
render_noyscale<OI, T, C, NORMF, RESAMPLER>
|
||||
(anl, msc, x_xf, xi0, xi1,
|
||||
ysi0, ysi1, p, xi1 - xi0, normf);
|
||||
|
||||
// Resample in the Y direction
|
||||
for (int64_t xi = xi0; xi < xi1; xi++) {
|
||||
transverse_source<RST, OI> src(render_data.data(),
|
||||
xi0, xi1, ysi0, ysi1,
|
||||
xi);
|
||||
stride_iterator<OI> dest(output + (xi - xi0), output_stride);
|
||||
y_resampler.resample(src, yi0, yi1, dest);
|
||||
// Resample in the Y direction
|
||||
for (int64_t xi = xi0; xi < xi1; xi++) {
|
||||
transverse_source<RST, OI> src(render_data.data(),
|
||||
xi0, xi1, ysi0, ysi1,
|
||||
xi);
|
||||
stride_iterator<OI> dest(output + (xi - xi0), output_stride);
|
||||
y_resampler(src, yi0, yi1, dest);
|
||||
}
|
||||
}
|
||||
updated(xi0, xi1, yi0, yi1);
|
||||
}
|
||||
|
|
@ -493,8 +493,8 @@ void render_p2scale(const analyzer<T> &anl,
|
|||
// Provide default inc_i0, inc_i1, and output_stride
|
||||
render_incremental<OI, T, C, NORMF, RESAMPLER>
|
||||
(anl, msc,
|
||||
linear_transform(ldexp(1, xe), xorigin),
|
||||
linear_transform(ldexp(1, ye), yorigin),
|
||||
affine_transform(ldexp(1, xe), xorigin),
|
||||
affine_transform(ldexp(1, ye), yorigin),
|
||||
xi0, xi1, yi0, yi1,
|
||||
INT64_MIN, INT64_MAX,
|
||||
output, xi1 - xi0, normf);
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// Fast resampling by powers of two
|
||||
//
|
||||
// Copyright (C) 2016-2020 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2016-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -17,7 +17,7 @@
|
|||
#include <math.h>
|
||||
#include <algorithm> // std::copy
|
||||
|
||||
#include "gaborator/linear_transform.h"
|
||||
#include "gaborator/affine_transform.h"
|
||||
#include "gaborator/pod_vector.h"
|
||||
|
||||
namespace gaborator {
|
||||
|
|
@ -56,7 +56,7 @@ namespace gaborator {
|
|||
struct p2_transform {
|
||||
p2_transform(int e_, int64_t origin_): e(e_), origin(origin_) { }
|
||||
// Convert a linear transform into a p2_transform
|
||||
p2_transform(linear_transform xf) {
|
||||
p2_transform(affine_transform xf) {
|
||||
int exp;
|
||||
double m = frexp(xf.a, &exp);
|
||||
assert(m == 0.5);
|
||||
|
|
@ -75,9 +75,9 @@ zoom_p2(p2_transform xf, int e) {
|
|||
return p2_transform(xf.e + e, xf.origin);
|
||||
}
|
||||
|
||||
static inline linear_transform
|
||||
zoom_p2(linear_transform xf, int e) {
|
||||
return linear_transform(ldexp(xf.a, e), xf.b);
|
||||
static inline affine_transform
|
||||
zoom_p2(affine_transform xf, int e) {
|
||||
return affine_transform(ldexp(xf.a, e), xf.b);
|
||||
}
|
||||
|
||||
// Resample data from "source", generating a view between indices
|
||||
|
|
@ -89,13 +89,14 @@ zoom_p2(linear_transform xf, int e) {
|
|||
// resolution (e=0).
|
||||
//
|
||||
// S is the type of the data source
|
||||
// T is the sample type
|
||||
// OI is the output iterator type
|
||||
|
||||
template <class S, class T>
|
||||
T *resample2_ptr(const S &source, p2_transform xf,
|
||||
int64_t i0, int64_t i1,
|
||||
bool interpolate, T *out)
|
||||
template <class S, class OI>
|
||||
OI resample2_p2xf(const S &source, p2_transform xf,
|
||||
int64_t i0, int64_t i1,
|
||||
bool interpolate, OI out)
|
||||
{
|
||||
typedef typename std::iterator_traits<OI>::value_type T;
|
||||
assert(i1 >= i0);
|
||||
if (xf.e > 0) {
|
||||
// Downsample
|
||||
|
|
@ -110,8 +111,8 @@ T *resample2_ptr(const S &source, p2_transform xf,
|
|||
// Get super-octave data
|
||||
gaborator::pod_vector<T> super_data(si1 - si0);
|
||||
T *p = super_data.data();
|
||||
p = resample2_ptr(source, p2_transform(xf.e - 1, xf.origin),
|
||||
si0, si1, interpolate, p);
|
||||
p = resample2_p2xf(source, p2_transform(xf.e - 1, xf.origin),
|
||||
si0, si1, interpolate, p);
|
||||
assert(p == super_data.data() + si1 - si0);
|
||||
for (int64_t i = i0; i < i1; i++) {
|
||||
int64_t si = 2 * i - si0;
|
||||
|
|
@ -180,8 +181,8 @@ T *resample2_ptr(const S &source, p2_transform xf,
|
|||
// Get sub-octave data
|
||||
gaborator::pod_vector<T> sub_data(si1 - si0);
|
||||
T *p = sub_data.data();
|
||||
p = resample2_ptr(source, p2_transform(xf.e + 1, xf.origin),
|
||||
si0, si1, interpolate, p);
|
||||
p = resample2_p2xf(source, p2_transform(xf.e + 1, xf.origin),
|
||||
si0, si1, interpolate, p);
|
||||
assert(p == sub_data.data() + si1 - si0);
|
||||
for (int64_t i = i0; i < i1; i++) {
|
||||
int64_t si = (i >> 1) - si0;
|
||||
|
|
@ -207,14 +208,10 @@ T *resample2_ptr(const S &source, p2_transform xf,
|
|||
return out;
|
||||
}
|
||||
|
||||
// As above, but generating the output through an iterator
|
||||
// rather than a pointer.
|
||||
//
|
||||
// S is the type of the data source
|
||||
// OI is the output iterator type
|
||||
// As above, but taking an affine_transform.
|
||||
|
||||
template <class S, class OI>
|
||||
OI resample2(const S &source, linear_transform lxf,
|
||||
OI resample2(const S &source, affine_transform lxf,
|
||||
int64_t i0, int64_t i1,
|
||||
bool interpolate, OI out)
|
||||
{
|
||||
|
|
@ -222,7 +219,7 @@ OI resample2(const S &source, linear_transform lxf,
|
|||
typedef typename std::iterator_traits<OI>::value_type T;
|
||||
gaborator::pod_vector<T> data(i1 - i0);
|
||||
T *p = data.data();
|
||||
p = resample2_ptr(source, xf, i0, i1, interpolate, p);
|
||||
p = resample2_p2xf(source, xf, i0, i1, interpolate, p);
|
||||
return std::copy(data.data(), data.data() + (i1 - i0), out);
|
||||
}
|
||||
|
||||
|
|
@ -234,13 +231,13 @@ OI resample2(const S &source, linear_transform lxf,
|
|||
// return an unnecessarily large support when interpolation is off
|
||||
|
||||
inline void
|
||||
resample2_support(linear_transform lxf, int64_t i0, int64_t i1,
|
||||
resample2_support(affine_transform lxf, int64_t i0, int64_t i1,
|
||||
int64_t &si0_ret, int64_t &si1_ret)
|
||||
{
|
||||
p2_transform xf(lxf);
|
||||
int margin = 2;
|
||||
if (xf.e > 0) {
|
||||
// Note code duplication wrt resample2_ptr().
|
||||
// Note code duplication wrt resample2_p2xf().
|
||||
// Also note tail recursion.
|
||||
int64_t si0 = i0 * 2 - margin + 1;
|
||||
int64_t si1 = i1 * 2 + margin;
|
||||
|
|
@ -261,7 +258,7 @@ resample2_support(linear_transform lxf, int64_t i0, int64_t i1,
|
|||
// destination indices that depend on a given range of source indices.
|
||||
|
||||
inline void
|
||||
resample2_inv_support(linear_transform lxf, int64_t si0, int64_t si1,
|
||||
resample2_inv_support(affine_transform lxf, int64_t si0, int64_t si1,
|
||||
int64_t &i0_ret, int64_t &i1_ret)
|
||||
{
|
||||
p2_transform xf(lxf);
|
||||
|
|
@ -290,9 +287,9 @@ resample2_inv_support(linear_transform lxf, int64_t si0, int64_t si1,
|
|||
// Lanczos2 power-of-two resampler
|
||||
|
||||
struct lanczos2_pow2_resampler {
|
||||
lanczos2_pow2_resampler(linear_transform xf_): xf(xf_) { }
|
||||
lanczos2_pow2_resampler(affine_transform xf_): xf(xf_) { }
|
||||
template <class S, class OI>
|
||||
OI resample(const S &source, int64_t i0, int64_t i1, OI out) const {
|
||||
OI operator()(const S &source, int64_t i0, int64_t i1, OI out) const {
|
||||
return resample2(source, xf, i0, i1, true, out);
|
||||
}
|
||||
void support(int64_t i0, int64_t i1, int64_t &si0_ret, int64_t &si1_ret) const {
|
||||
|
|
@ -301,16 +298,16 @@ struct lanczos2_pow2_resampler {
|
|||
void inv_support(int64_t si0, int64_t si1, int64_t &i0_ret, int64_t &i1_ret) {
|
||||
return resample2_inv_support(xf, si0, si1, i0_ret, i1_ret);
|
||||
}
|
||||
linear_transform xf;
|
||||
affine_transform xf;
|
||||
};
|
||||
|
||||
// Nearest-neighbor power-of-two resampler
|
||||
// XXX simplify
|
||||
|
||||
struct nn_pow2_resampler {
|
||||
nn_pow2_resampler(linear_transform xf_): xf(xf_){ }
|
||||
nn_pow2_resampler(affine_transform xf_): xf(xf_){ }
|
||||
template <class S, class OI>
|
||||
OI resample(const S &source, int64_t i0, int64_t i1, OI out) const {
|
||||
OI operator()(const S &source, int64_t i0, int64_t i1, OI out) const {
|
||||
return resample2(source, xf, i0, i1, false, out);
|
||||
}
|
||||
void support(int64_t i0, int64_t i1, int64_t &si0_ret, int64_t &si1_ret) const {
|
||||
|
|
@ -319,7 +316,7 @@ struct nn_pow2_resampler {
|
|||
void inv_support(int64_t si0, int64_t si1, int64_t &i0_ret, int64_t &i1_ret) {
|
||||
return resample2_inv_support(xf, si0, si1, i0_ret, i1_ret);
|
||||
}
|
||||
linear_transform xf;
|
||||
affine_transform xf;
|
||||
};
|
||||
|
||||
} // namespace
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
//
|
||||
// Version number
|
||||
//
|
||||
// Copyright (C) 2015-2020 Andreas Gustafsson. This file is part of
|
||||
// Copyright (C) 2015-2021 Andreas Gustafsson. This file is part of
|
||||
// the Gaborator library source distribution. See the file LICENSE at
|
||||
// the top level of the distribution for license information.
|
||||
//
|
||||
|
|
@ -10,6 +10,6 @@
|
|||
#define _GABORATOR_VERSION_H
|
||||
|
||||
#define GABORATOR_VERSION_MAJOR 1
|
||||
#define GABORATOR_VERSION_MINOR 6
|
||||
#define GABORATOR_VERSION_MINOR 7
|
||||
|
||||
#endif
|
||||
|
|
|
|||
Loading…
Reference in a new issue