Use standard intrinsic types, adjust __builtin with intrinsics

Release 2.0
2023-10-15 16:43:27 +02:00 · 2023-10-15 16:42:31 +02:00
33 changed files with 2082 additions and 1183 deletions
--- a/30
+++ b/30
@ -1,4 +1,34 @@

+2.0
+
+Add new analyzer methods bandpass_band_ff() and band_q().
+
+Add new analyzer methods band_analysis_support() and
+band_synthesis_support().
+
+Add support for adjusting the amount of overlap between bands.
+
+Add support for selecting a local or global phase convention, and make
+the local phase convention the default in the 2.0 API.
+
+Extend the supported ff_min range to 0.0001...0.5.
+
+Add support for linear and mel frequency scales.
+
+Analyzers no longer share data via static variables as this could
+cause race conditions in multithreaded applications.
+
+Disable the copy constructor and assignment operator of the analyzer
+class.
+
+The number of bands per octave is no longer restricted to integer
+values.
+
+Further improve the performance of analyzing short signal blocks.
+
+Document bands_begin() and bands_end(), and expand the documentation
+of band numbering in general.
+
 1.7

 Miscellaneous bug fixes.
--- a/2
+++ b/2
@ -1,5 +1,5 @@

-The Gaborator library is Copyright (C) 1992-2019 Andreas Gustafsson.
+The Gaborator library is Copyright (C) 1992-2023 Andreas Gustafsson.

 License to distribute and modify the code is hereby granted under the
 terms of the GNU Affero General Public License, version 3 (henceforth,
--- a/doc/._filter-response.png
+++ b/doc/._filter-response.png
--- a/doc/filter.html
+++ b/doc/filter.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2017-2020 Andreas Gustafsson.  This file is part of
+Copyright (C) 2017-2023 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.
 -->
@ -77,14 +77,17 @@ that in <a href="render.html">Example 1</a>:</p>

 <h2>Spectrum Analysis Parameters</h2>

-<p>The spectrum analysis works much the same as in Example 1,
-but uses slightly different parameters.
-We use a larger number of frequency bands per octave (100)
-to minimize ripple in the frequency response, and the
-reference frequency argument is omitted as we don't care about the
-exact alignment of the bands with respect to a musical scale.</p>
+<p>The spectrum analysis works much the same as in Example 1, but uses
+slightly different parameters.  The number of frequency bands per
+octave determines how steep the transitions between filter bands can
+be made, and since this particular filter has no steep transitions, a
+low value like 10 is sufficient.  The reference frequency argument is
+omitted as we don't care about the exact alignment of the bands with
+respect to a musical scale.</p>
+
 <pre>
-    gaborator::parameters params(100, 20.0 / fs);
+    gaborator::log_fq_scale scale(10, 20.0 / fs);
+    gaborator::parameters params(scale);
    gaborator::analyzer&lt;float&gt; analyzer(params);
 </pre>

@ -107,7 +110,7 @@ proportional to the square root of the inverse of the frequency.
 To get the frequency of each band, we call the
 <code>analyzer</code> method <code>band_ff()</code>, which
 returns the center frequency of the band in units of the
-sampling frequency.  The gain is normalized to unity at 20 Hz.
+sample rate.  The gain is normalized to unity at 20 Hz.
 </p>
 <pre>
    for (int band = analyzer.bandpass_bands_begin(); band &lt; analyzer.bandpass_bands_end(); band++) {
@ -201,7 +204,7 @@ Note that we use <code>SFC_SET_CLIPPING</code>
 to make sure that any samples too loud for the file format
 will saturate; by default, <i>libsndfile</i> makes them
 wrap around, which sounds really bad.</p>
-<a name="writing_audio_code">
+<a id="writing_audio_code">
 <pre>
    SNDFILE *sf_out = sf_open(argv[2], SFM_WRITE, &amp;sfinfo);
    if (! sf_out) {
@ -233,16 +236,16 @@ 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 FreeBSD">
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` filter.cc `pkg-config --libs sndfile` -o filter
+<pre class="build Darwin Linux NetBSD FreeBSD OpenBSD">
+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>
 <pre class="build Darwin">
-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
+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">
-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
+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>
--- a/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png
+++ b/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_anl_wob.png
--- a/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png
+++ b/doc/gen/allkernels_v1_bpo12_ffmin0.03_ffref0.5_syn_wob.png
--- a/doc/gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png
+++ b/doc/gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png
--- a/doc/index.html
+++ b/doc/index.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2018-2021 Andreas Gustafsson.  This file is part of
+Copyright (C) 2018-2023 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.
 -->
@ -12,23 +12,16 @@ the top level of the distribution for license information.
 </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 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&ouml;rfler, and Grill, described in the papers
-<i><a href="http://www.univie.ac.at/nonstatgab/pdf_files/dohogrve11_amsart.pdf">
-Constructing an invertible constant-Q transform with nonstationary Gabor frames, 2011</a></i>
-and <i><a href="http://www.univie.ac.at/nonstatgab/pdf_files/dogrhove12_amsart.pdf">
-A Framework for invertible, real-time constant-Q transforms, 2012</a></i>,
-using Gaussian bandpass filters and an efficient multi-rate architecture.
-</p>
+<p>The Gaborator is a C++ library that generates spectrograms for the
+visualization and analysis of audio signals.  It supports constant-Q,
+constant-bandwidth and mel scale spectrograms, as well as a fast and
+accurate inverse transformation of the spectrogram coefficients back
+into audio for spectral effects and editing.</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, FreeBSD, and iOS, on Intel
-x86_64 and ARM processors.</p>
+<p>The Gaborator is compatible with C++ versions from C++11 through
+C++20 and beyond.  It has been tested on macOS, iOS, Linux, NetBSD,
+FreeBSD, and OpenBSD, 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.
@ -38,7 +31,7 @@ See the file <a href="../LICENSE">LICENSE</a> for details.</p>

 <p>The following examples demonstrate the use of the library in
 various scenarios.  They are presented in a "literate
-programming" style, with the code embedded in the commentary
+programming" style, where the code is embedded in the commentary
 rather than the other way around.
 Concatenating the code fragments in each example yields a complete C++
 program, which can also be found as a <code>.cc</code> file in
--- a/doc/overview.html
+++ b/doc/overview.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2018-2021 Andreas Gustafsson.  This file is part of
+Copyright (C) 2018-2023 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.
 -->
@ -33,20 +33,21 @@ of these functions.</p>
 an <i>analysis filter bank</i>, to split it into a number of
 overlapping frequency <i>bands</i>.</p>

-<p>The filter bank consists of a number of logarithmically spaced
-Gaussian bandpass filters and a single lowpass filter.  Each bandpass
-filter has a bandwidth proportional to its center frequency, which
-means they all have the same quality factor Q and form
-a <i>constant-Q</i> filter bank.  The highest-frequency bandpass
-filter will have a center frequency close to half the sample rate; in
-the graphs below, this is simple labeled 0.5 because all frequencies
-in the Gaborator are in units of the sample rate.  The
-lowest-frequency bandpass filter should be centered at, or slightly
-below, the lowest frequency of interest to the application at hand.
-For example, when analyzing audio, this is often the lower limit of
-human hearing; at a sample rate of 44100 Hz, this means 20 Hz / 44100
-Hz &asymp; 0.00045.  This lower frequency limit is referred to as
-the <i>minimum frequency</i> or f<sub>min</sub>.
+<p>When using a logarithmic frequency scale,
+the filter bank consists of a number of
+logarithmically spaced Gaussian bandpass filters and a single lowpass
+filter.  Each bandpass filter has a bandwidth proportional to its
+center frequency, which means they all have the same quality factor Q
+and form a <i>constant-Q</i> filter bank.  The highest-frequency
+bandpass filter will have a center frequency close to half the sample
+rate. In the graphs below, this is labeled 0.5 because
+frequencies in the Gaborator are generally given in units of the
+sample rate.  The lowest-frequency bandpass filter should be centered
+at, or slightly below, the lowest frequency of interest to the
+application at hand.  For example, when analyzing audio, this is often
+the lower limit of human hearing; at a sample rate of 44100 Hz, this
+means 20 Hz / 44100 Hz &asymp; 0.00045.  This lower frequency limit is
+referred to as the <i>minimum frequency</i> or f<sub>min</sub>.
 </p>

 <p>Although frequencies below f<sub>min</sub> are assumed to not be of
@ -55,12 +56,27 @@ reconstruction, and that is what the lowpass filter is for.  Together,
 the lowpass filter and the bandpass filters overlap to cover the full
 frequency range from 0 to 0.5.</P>

-<p>The spacing of the bandpass filters is specified by the user as an
-integer number of filters (or, equivalently, bands) per octave.  For
+<p>The spacing of the bandpass filters is specified by the user as
+a number of filters (or, equivalently, bands) per octave.  For
 example, when analyzing music, this is often 12 bands per octave (one
 band per semitone in the equal-tempered scale), or if a finer
 frequency resolution is needed, some multiple of 12.</p>

+<p id="overlap">The bandwidth of each individual bandpass filter is
+chosen to achieve a reasonable amount of overlap with the adjacent
+filters.  If the bandwidth is too narrow, there will be too little
+overlap, causing deep gaps between the bands.  If it is too wide,
+there will be a great deal of overlap, resulting in a blurred
+spectrogram with poor frequency selectivity and highly redundant
+coefficients.</p>
+
+<p>Since the Gaborator uses Gaussian bandpass filters, it defines the
+width of each filter in terms of its standard deviation.  The overlap
+is defined as the ratio of this standard deviation to the spacing between
+adjacent bands.  The default value for the overlap is 0.7, meaning the
+standard deviation of each Gaussian filter is 0.7 times the local
+spacing between adjacent filters.</p>
+
 <p>The following plot shows the frequency responses of the analysis
 filters at 12 bands per octave and f<sub>min</sub> = 0.03.  A more
 typical f<sub>min</sub> for audio work would be 0.00045, but
@ -69,21 +85,25 @@ and the lowest-frequency bandpass filters would be extremely narrow.</p>

 <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
-at a reduced sample rate, lower than that of the orignal signal but
-high enough that there is negligible aliasing given the bandwidth of
-the filter in case.  The Gaborator uses sample rates related to the
-original signal sample rate by powers of two.  This means some of
-frequency bands are sampled a bit more often than strictly
-necessary, but has the advantage that the sampling can be synchronized
-to make the samples of many frequency bands coincide in time, which
-can be convenient in later analysis or spectrogram rendering.  The
-complex samples resulting from this process are the spectrogram
-coefficients.</p>
+<p>The bandpass filters produce a complex-valued output representing
+the amplitude and the phase of the signal within each band, and sampling
+this output produces the final spectrogram coefficients.
+Since the bandwidth of an individual band is smaller than that of
+the input signal as a whole, it can be sampled at a reduced sample rate.
+</p>

-<p>The center frequencies of the analysis filters and the points in
-time at which they are sampled form a two-dimensional,
+<p>To minimize the amount of coefficient data, each band should in
+principle be sampled at a different sample rate, but dealing with a
+large number of different sample rates would be cumbersome.  Instead,
+all bands are sampled at rates that are the input sample rate divided by
+some power two, oversampling the coefficients at the next higher such
+rate as needed.  This also has the advantage that the sampling can be
+synchronized to make the samples of many frequency bands coincide in
+time, which can be convenient in later analysis or spectrogram
+rendering.</p>
+
+<p>The center frequencies of the bands and the sample points in
+time together form a two-dimensional,
 multi-resolution <i>time-frequency grid</i>, where high frequencies
 are sampled sparsely in frequency but densely in time, and low
 frequencies are sampled densely in frequency but sparsely in time.</p>
@ -97,21 +117,35 @@ positive and the negative direction.</p>

 <img src="gen/grid_v1_bpo12_ffmin0.03_ffref0.5_wob.png" alt="Sampling grid">

+<p>When using a linear or mel frequency scale, no special lowpass band
+is needed because frequency scale extends to zero.
+In the case of a linear frequency scale, no multirate processing is
+needed, either, and the sampling grid is uniformly spaced in both the
+time and frequency dimensions.</p>
+
 <h2>Resynthesis</h2>

 <p>Resynthesizing a signal from the coefficients is more or less the
-reverse of the analysis process.  The coefficients are frequency
-shifted from the complex baseband back to their original center
-frequencies and run through a <i>reconstruction filter bank</i>
-that is a <i>dual</i> of the analysis filter bank.  The following
+reverse of the analysis process.  The coefficients are upsampled
+to the original signal sample rate and run through a <i>reconstruction filter bank</i>
+that is a <i>dual</i> of the analysis filter bank.  The construction
+of the dual filters is based on the methods described by
+Velasco, Holighaus, D&ouml;rfler, and Grill in the papers
+<i><a href="http://www.univie.ac.at/nonstatgab/pdf_files/dohogrve11_amsart.pdf">
+Constructing an invertible constant-Q transform with nonstationary Gabor frames, 2011</a></i>
+and <i><a href="http://www.univie.ac.at/nonstatgab/pdf_files/dogrhove12_amsart.pdf">
+A Framework for invertible, real-time constant-Q transforms, 2012</a></i>.
+</p>
+
+<p>The following
 plot shows the frequency responses of the reconstruction filters
 corresponding to the analysis filters shown earlier.</p>

 <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
-different.</p>
+<p>Although the bandpass filters look superficially similar to the
+Gaussian filters of the analysis filter bank, their shapes are
+actually subtly different.</p>

 <h2>Spectrogram Rendering</h2>

--- a/doc/realtime.html
+++ b/doc/realtime.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2018-2020 Andreas Gustafsson.  This file is part of
+Copyright (C) 2018-2021, 2023 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.
 -->
@ -98,7 +98,7 @@ For example, as of version 1.4, analyzing a signal in blocks of
 2<sup>10</sup> samples takes roughly five times as much CPU as
 analyzing it in blocks of 2<sup>20</sup> samples.</p>

-<p>For sufficiently small blocks, the processing time will exceed the
+<p>For very small blocks, the processing time will exceed the
 duration of the signal, at which point the system can no longer be
 considered real-time.  For example, analyzing a 48&nbsp;kHz audio
 stream on a 2.5 GHz Intel Core i5 CPU, this happens at block sizes
@ -111,9 +111,15 @@ blocks.</p>

 <h2>Can it process a signal stream of any length?</h2>

-<p>Not in practice &mdash; the length is limited by floating point
-precision.  At typical audio sample rates, roundoff errors start to
-become significant after some hours.</p>
+<p>As of version 2, the length of a stream is only limited by the
+64-bit signed integers used to represent sample times, allowing for
+more than a million years of 192 kHz audio both before and after t=0.
+This assumes the local phase convention is used (which is the default
+in version 2).  When using the global phase convention (as in version
+1), the floating point precision of phase values used in internal
+calculations will decrease as the distance from t=0 increases, causing
+a noticeable decrease in the signal to noise ratio after some
+hours.</p>

 <h2>Does it avoid dynamic memory allocation in the audio processing path?</h2>

--- a/doc/ref/gaborator_h.html
+++ b/doc/ref/gaborator_h.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2018-2021 Andreas Gustafsson.  This file is part of
+Copyright (C) 2018-2023 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.
 -->
@ -12,34 +12,76 @@ the top level of the distribution for license information.
 <body>
 <h1>Gaborator reference: <code>gaborator.h</code></h1>

-<h2>Spectrum Analysis Parameters</h2>
+<pre class="forward_decl">
+template&lt;class T&gt; class analyzer;
+</pre>

-<p>A <code>parameters</code> object holds a set of parameters that
-determine the frequency range and resolution of the spectrum
-analysis.</p>
+<h2>Frequency scales</h2>
+
+<p>The Gaborator supports three types of frequency scales:
+logarithmic, linear, and mel scales.  Each is represented by a class
+deriving from the abstract base class <code>fq_scale</code>.</p>

 <pre>
-class parameters {
+class fq_scale {
 </pre>

 <div class="class_def">
-<h3>Constructor</h3>
+
+<h3>Comparison</h3>
+<p>
+Comparison operators are provided for compatibility with
+standard container classes.
+</p>
 <pre>
-parameters(unsigned int bands_per_octave,
-           double ff_min,
-           double ff_ref = 1.0);
+bool operator<(const fq_scale &amp;rhs) const;
+bool operator==(const fq_scale &amp;rhs) const;
 </pre>
+
+</div>
+
+<pre>
+};
+</pre>
+
+<h3>Logarithmic frequency scale</h3>
+
+<p>The <code>log_fq_scale</code> class represents a logarithmic
+frequency scale, corresponding to a constant-Q transform or
+constant-Q spectrogram.</p>
+
+<pre>
+class log_fq_scale: public fq_scale {
+</pre>
+
+<div class="class_def">
+
+<h4>Constructor</h4>
+
+<pre>
+log_fq_scale(double bands_per_octave,
+             double ff_min,
+             double ff_ref = 0.5);
+</pre>
+
 <dl>
  <dt><code>bands_per_octave</code></dt>
  <dd>The number of frequency bands per octave.
  Values from 4 to 384 (inclusive) are supported.
+  Integer values yield the best performance, but
+  non-integer values are also supported.
+  The minimum of 4 applies when the default overlap of 0.7 is used;
+  if the overlap is increased from the default, the minimum should be
+  increased proportionally.
+  Values greater than 384 probably work but are not regularly tested.
  </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
-  frequency that <code>ff_min</code> falls between the two lowest
-  frequency bandpass filters.
-  Values from 0.001 to 0.13 are supported.</dd>
+  frequency that the lowest bandpass filter center frequency falls at
+  or below <code>ff_min</code>.
+  Values from 0.0001 up to but not including 0.5 are supported.
+  Values below 0.0001 probably work but are not regularly tested.</dd>
  <dt><code>ff_ref</code></dt>
  <dd>The reference frequency, in units of the sample rate.
  This allows fine-tuning of the analysis and synthesis filter
@ -50,12 +92,204 @@ parameters(unsigned int bands_per_octave,
  <code>ff_ref</code>.  Must be positive.  A typical value
  when analyzing music is <code>440.0 / fs</code>, where
  <code>fs</code> is the sample rate in Hz.
+  </dd>
+</dl>
+
+</div>
+
+<pre>
+};
+</pre>
+
+<h3>Linear frequency scale</h3>
+
+<p>The <code>lin_fq_scale</code> class represents a linear
+frequency scale, corresponding to a short-term Fourier transform
+(STFT) or constant-bandwidth spectrogram.
+</p>
+
+<pre>
+class lin_fq_scale: public fq_scale {
+</pre>
+
+<div class="class_def">
+
+<h4>Constructor</h4>
+
+<pre>
+lin_fq_scale(double size);
+</pre>
+
+<dl>
+<dt><code>size</code></dt>
+
+<dd><p style="margin-top: 0pt">The size of the transform.  This
+corresponds to the "FFT size" parameter that is typically used with
+FFT based spectrum analyzers: the frequency spacing between bands is
+the sample rate divided by <i>size</i>, and the number of frequency
+bands produced is roughly half of <i>size</i> because negative
+frequencies are omitted.  However, in the Gaborator, the <i>size</i>
+parameter is not limited to powers of two, or even to integer values,
+so any arbitrary band spacing of <i>x</i> Hz can be achieved by
+specifying a <i>size</i> of <code>fs / x</code> where
+<code>fs</code> is the sample rate in Hz.</p>
+
+<p>Values from 2 to 10000 are supported.  Larger values probably work but
+are not regularly tested.</p>
+
 </dd>
 </dl>
+
+</div>
+
+<pre>
+};
+</pre>
+
+<h3>Mel frequency scale</h3>
+
+<p>The <code>mel_fq_scale</code> class represents a mel
+frequency scale, a scale designed to mimic human pitch perception
+such that one mel is close to the smallest percepible difference
+in frequency.
+It is defined by the formula <i>mel = 2595
+log<sub>10</sub>(1 + f / 700)</i>, where <i>f</i>&thinsp;
+is the frequency in Hz.</p>
+
+<p>The mel scale is approximately logarithmic
+at high frequencies, and approximately linear at low frequencies,
+with a smooth transition between the two.</p>
+
+<pre>
+class mel_fq_scale: public fq_scale {
+</pre>
+
+<div class="class_def">
+<h4>Constructor</h4>
+
+<pre>
+mel_fq_scale(double bands_per_mel,
+             double fs);
+</pre>
+
+<dl>
+<dt><code>bands_per_mel</code></dt>
+
+<dd>The number of frequency bands per mel.  Supported values range
+from 0.003 (yielding a total of 12 frequency bands between 0 and
+22.5&nbsp;kHz at a 44.1&nbsp;kHz sample rate) to 1.0 (yielding some
+3924 bands).  The minimum of 0.003 applies when the default overlap of
+0.7 is used; if the overlap is increased from the default, the minimum
+should be increased proportionally.  Values larger than 1.0 probably
+work but are not regularly tested.</dd>
+
+<dt><code>fs</code></dt>
+
+<dd>The sample rate in Hz.  Because the mel scale is defined in terms
+of absolute frequencies in Hz, the sample rate needs to be specified
+explicitly; this is an exception to the sample rate agnostic design
+of the rest of the Gaborator API.  Values from 8000 to 96000 are
+supported.  Values outside this range probably work but are not
+regularly tested.
+</dd>
+</dl>
+
+</div>
+
+<pre>
+};
+</pre>
+
+<h2 id="phase_conventions">Phase conventions</h2>
+
+<pre>
+enum class coef_phase { global, local };
+</pre>
+
+<p>The spectrogram cofficients produced by the Gaborator are complex
+numbers, and as such they have both a magnitude and a phase.  The
+phase can be defined in several different ways, using different
+reference points for what is considered a phase of zero.</p>
+
+<p>Under one possible definition, an input signal of a real and
+positive impulse at some time <i>t</i>&thinsp; will cause all
+spectrogram coefficients at that time <i>t</i>&thinsp; to also be real
+and positive, or in other words, to have a phase of zero, for any
+frequency.  We call this the <i>local</i> phase convention, and it is
+the default beginning with version 2 of the Gaborator.<p>
+
+<p>Under an alternative definition, an input signal of a cosine wave
+at some frequency <i>f</i>&thinsp; will cause all spectrogram coefficients at
+that frequency <i>f</i>&thinsp; to have a phase of zero, at any point in time.
+We call this the <i>global</i> phase convention, and it was the
+default in version 1 of the Gaborator.</p>
+
+<p>The local phase convention is preferred because it guarantees that
+the coefficients for a given point in time only depend on the local
+signal behavior around that point in time, and not on how far that
+point in time is from the reference point of t=0.  This avoids a
+loss of floating point precision for signals far from t=0.</p>
+
+<h2>Spectrum Analysis Parameters</h2>
+
+<p>A <code>parameters</code> object holds a set of parameters for
+the spectrum analysis and resynthesis.</p>
+
+<pre>
+class parameters {
+</pre>
+
+<div class="class_def">
+<h3>Constructor</h3>
+
+<p>The only required parameter is the frequency scale, which is passed
+as an argument to the <code>parameters</code> constructor.  Optional
+parameters may be specified by assigning to data members of
+the <code>parameters</code> object after construction.</p>
+
+<pre>
+parameters(const fq_scale &scale);
+</pre>
+<h3>Legacy Constructor</h3>
+<p>For backwards compatibilty with version 1 of the Gaborator,
+the following constrctor is also supported:</p>
+<pre>
+parameters(double bands_per_octave,
+           double ff_min,
+           double ff_ref = 1.0);
+</pre>
+<p>This constructs a set of parameters with a logarithmic frequency scale, like
+<code>parameters(log_fq_scale(bands_per_octave, ff_min, ff_ref))</code>
+except that the <code>phase</code> member is initialized to the version 1 default
+of <code>coef_phase::global</code> rather than the verison 2 default of
+<code>coef_phase::local</code>.
+</p>
+
+<h3>Phase convention</h3>
+
+<p>The <a href="#phase_conventions">phase convention</a>
+to be used for the analysis and resynthesis
+can be set or examined using the public data member <code>phase</code>.
+
+<pre>
+coef_phase phase;
+</pre>
+
+<h3>Overlap</h3>
+
+<p>The amount of <a href="../overview.html#overlap">overlap</a> between
+adjacent bandpass bands can be set or examined using the public data
+member <code>overlap</code>.  The default is 0.7, meaning each band
+will have a standard deviation of 0.7 times the band spacing.</p>
+
+<pre>
+double overlap;
+</pre>
+
 <h3>Comparison</h3>
 <p>
 Comparison operators are provided for compatibility with
-standard container classes.  The ordering is arbitrary but consistent.
+standard container classes.
 </p>
 <pre>
 bool operator<(const parameters &amp;rhs) const;
@ -69,10 +303,6 @@ bool operator==(const parameters &amp;rhs) const;

 <h2>Spectrogram Coefficients</h2>

-<pre class="forward_decl">
-template&lt;class T&gt; class analyzer;
-</pre>
-
 <p>
 A <code>coefs</code> object stores a set of spectrogram coefficients.
 It is a dynamic data structure and will be automatically grown to
@ -85,7 +315,7 @@ it will default to <code>std::complex&lt;T&gt;</code>.
 </p>

 <pre>
-template&lt;class T, class C = std::complex&lt;T&gt&gt;
+template &lt;class T, class C = std::complex&lt;T&gt&gt;
 class coefs {
 </pre>
 <div class="class_def">
@ -113,7 +343,7 @@ the floating-point type to use for the calculations.  This is typically <code>fl
 alternatively, <code>double</code> can be used for increased accuracy at the
 expense of speed and memory consumption.</p>

-<pre>template&lt;class T&gt;
+<pre>template &lt;class T&gt;
 class analyzer {</pre>
 <div class="class_def">

@ -175,12 +405,12 @@ etc.
 <pre>
 void
 synthesize(const coefs&lt;T&gt; &amp;coefs,
-           uint64_t t0,
-           uint64_t t1,
+           int64_t t0,
+           int64_t t1,
           T *signal) const;
 </pre>
-<p>Synthesize signal samples from the coefficients <code>coef</code> and store
-them at <code>*signal</code>.
+<p>Synthesize signal samples from the coefficients <code>coefs</code>
+and store them at <code>*signal</code>.
 </p>
 <dl>
  <dt><code>coefs</code></dt><dd>The coefficients to synthesize the signal from.</dd>
@ -207,73 +437,155 @@ be multiple of a large power of two.<p>
 <h3>Frequency Band Numbering</h3>

 <p>The frequency bands of the analysis filter bank are numbered by
-nonnegative integers that increase towards lower (sic) frequencies.
-There is a number of <i>bandpass bands</i> corresponding to the
+non-negative integers that increase towards lower (sic) frequencies.
+Although this numbering may seem backwards, it makes a certain
+amount of sense for constant-Q analysis where the
+frequency range has a hard upper limit (the Nyquist frequency) but no
+hard lower limit.  Numbering the bands from high to low frequency
+allows the frequency range to be extended into arbitrarily low
+frequencies without using negative band numbers or having to renumber
+existing bands.  It also maps logically to computer graphics
+coordinate systems where the Y coordinate increases downwards.</p>
+
+<p>When using a logarithmic frequency scale,
+there is a number of <i>bandpass bands</i> corresponding to the
 logarithmically spaced bandpass analysis filters, from near 0.5
-(half the sample rate) to
+(the Nyquist frequency, or half the sample rate) to
 near f<sub>min</sub>, and a single <i>lowpass band</i> containing the
 residual signal from frequencies below f<sub>min</sub>.
-The numbering can be examined using the following methods:
+Note that there is no special highpass band; frequencies at or
+close to the Nyquist frequency are included in the highest-frequency
+bandpass band(s).
 </p>

+<p>When using a linear or mel frequency scale, all the bands are
+considered <i>bandpass bands</i> and there is no special <i>lowpass band</i>.
+</p>
+
+<p>
+The band numbering can be examined using the following methods:
+</p>
+
+<pre>
+int bands_begin() const;
+</pre>
+<p>
+Return the lowest valid band number.  This is always 0.
+</p>
+<pre>
+int bands_end() const;
+</pre>
+<p>
+Return the highest valid band number plus one.  This is also the
+total number of bands.
+</p>
 <pre>
 int bandpass_bands_begin() const;
 </pre>
 <p>
-Return the smallest valid bandpass band number, corresponding to the
-highest-frequency bandpass filter.</p>
+Return the lowest valid bandpass band number, corresponding to the
+highest-frequency bandpass band.  This is currently 0, but should
+a future version of the Gaborator support a highpass band, that will
+become band 0 and <code>bandpass_bands_begin()</code> will return 1.
+</p>
 <pre>
 int bandpass_bands_end() const;
 </pre>
 <p>
-Return the bandpass band number one past the highest valid bandpass
-band number, corresponding to one past the lowest-frequency bandpass
-filter.
+Return the highest valid bandpass band number plus one, corresponding
+to one past the lowest-frequency bandpass band.
 </p>
 <pre>
 int band_lowpass() const;
 </pre>
 <p>
-Return the band number of the lowpass band.
+For an anzlyer using a logarithmic frequency scale only,
+return the band number of the lowpass band.  This returns the same
+number as <code>bandpass_bands_end()</code>, but is preferred for
+clarity when referring to the lowpass band rather than
+the excluded upper bound of the range of bandpass bands.
 </p>
 <pre>
 int band_ref() const;
 </pre>
 <p>
-Return the band number corresponding to the reference frequency
+For an anzlyer using a logarithmic frequency scale only,
+return the band number corresponding to the reference frequency
 <code>ff_ref</code>.   If <code>ff_ref</code> falls within
-the frequency range of the bandpass filter bank, this will
-be a valid bandpass band number, otherwise it will not.
+the frequency range of the bandpass filter bank, this is
+a valid bandpass band number, otherwise it is not.
 </p>
+
+<h3>Band properties</h3>
+
 <pre>
 double band_ff(int band) const;
 </pre>
 <p>
 Return the center frequency of band number <i>band</i>, in units of the
-sampling frequency.
+sample rate.  The center frequency of the lowpass band (if present) is 0.
+</p>
+
+<pre>
+double bandpass_band_ff(double band) const;
+</pre>
+<p>
+Return the center frequency of bandpass band number <i>band</i>, in
+units of the sample rate.  Unlike <code>band_ff</code>, this takes a
+floating point argument and supports interpolating between bands by
+giving a non-integer band number, and/or extrapolating outside the valid
+range of bandpass band numbers.
+</p>
+
+<pre>
+double band_q(double band) const;
+</pre>
+<p>
+Return the Q (quality factor) of the analysis filter of bandpass
+band <i>band</i>.  Q is defined as the -3&nbsp;dB bandwidth divided by
+the center frequency.  When using a logarithmic frequency scale, the
+bands form a constant-Q filter bank and Q is the same for all bands,
+and when using a linear or mel frequency scale, Q will vary from band
+to band.
 </p>

 <h3>Support</h3>
+
 <pre>
-double analysis_support() const;
+double band_analysis_support(int band) const;
 </pre>
-<p>Returns the one-sided worst-case time domain <i>support</i> of any of the
-analysis filters. When calling <code>analyze()</code> with a sample at time <i>t</i>,
-only spectrogram coefficients within the time range <i>t &plusmn; support</i>
-will be significantly changed.  Coefficients outside the range may change,
+<p>Return the one-sided time domain <i>support</i> of any of the analysis
+filter for band <i>band</i>.
+When calling <code>analyze()</code> with a sample at time <i>t</i>,
+the spectrogram coefficients of band <i>band</i> will change
+significantly only within the time range <i>t &plusmn; support</i>.
+Coefficients outside the range may change slightly,
 but the changes will sufficiently small that they may be ignored without
 significantly reducing accuracy.</p>

+<pre>
+double analysis_support() const;
+</pre>
+<p>Return the largest <code>band_analysis_support(band)</code> of any band.
+</p>
+
+<pre>
+double band_synthesis_support(int band) const;
+</pre>
+<p>Returns the one-sided time domain <i>support</i> of the
+reconstruction filter for band <i>band</i>.  When
+calling <code>synthesize()</code> to synthesize a sample at
+time <i>t</i>, the sample will only be significantly affected by
+spectrogram coefficients of band <i>band</i> in the time range <i>t
+&plusmn; support</i>.  Coefficients outside the range may be used in
+the synthesis, but substituting zeroes for the actual coefficient
+values will not significantly reduce accuracy.</p>
+
 <pre>
 double synthesis_support() const;
 </pre>
-<p>Returns the one-sided worst-case time domain <i>support</i> of any of the
-reconstruction filters.  When calling <code>synthesize()</code> to
-synthesize a sample at time <i>t</i>, the sample will only be
-significantly affected by spectrogram coefficients in the time
-range <i>t &plusmn; support</i>.  Coefficients outside the range may
-be used in the synthesis, but substituting zeroes for the actual
-coefficient values will not significantly reduce accuracy.</p>
+<p>Return the largest <code>band_synthesis_support(band)</code> of any band.
+</p>

 </div>
 <pre>
@ -312,7 +624,7 @@ and <code>t1</code>, respectively.
 <p>The function <code>f</code> should have the call signature</p>
 <dd>
 <pre>
-template&lt;class T&gt;
+template &lt;class T&gt;
 void f(int b, int64_t t, std::complex&lt;T&gt; &c0, std::complex&lt;T&gt; &ci...);
 </pre>
    <p>where</p>
@ -344,13 +656,19 @@ coefficient <i>is</i> iterated over and is missing from one of the additional
 coefficient sets <code>ci...</code>, it will be automatically created
 and initialized to zero in that additional coefficient set.</p>

+<p>For example, in the common case where the processing takes one input
+and produces one output, <code>c0</code> should be the input
+and <code>c1</code> should be the output.  This ensures that the
+entire input is iterated over, and that the output coefficients get
+created as needed.</p>
+
 <p><i>Note: The template parameters <code>C0</code>
 and <code>CI</code>...  exist to support the processing of coefficient
 sets containing data of types other
 than <code>std::complex&lt;T&gt;</code>, which is not currently part of the
-documented API.  In typical use, there is no need to specify them when
-calling <code>apply()</code> because the template parameter list
-can be deduced, but if they are expicitly specified, they should all
+documented API.  In typical use, there is no need to specify them
+because the template parameter list
+can be deduced, but if they are explicitly specified, they should all
 be <code>std::complex&lt;T&gt;</code>.
 </i></p>

@ -401,8 +719,8 @@ bounded.</p>

 <h3>Legacy API For Iterating Over Existing Coefficients</h3>

-<p>Prior to version 1.5, the only way to iterate over
-coefficients was the <code>apply()</code> function.
+<p>Prior to version 1.5, iteration over the coefficients was done
+using the <code>apply()</code> function.
 It is similar to <code>process()</code>, except that it
 </p>
 <ul>
@ -439,7 +757,7 @@ is applied to every coefficient.
  <dd>A function to apply to each coefficient in <code>c</code>,
  with the call signature
 <pre>
-template&lt;class T&gt;
+template &lt;class T&gt;
 void f(std::complex&lt;T&gt; &amp;coef, int band, int64_t t);
 </pre>
    <dl>
@ -450,10 +768,10 @@ void f(std::complex&lt;T&gt; &amp;coef, int band, int64_t t);
        This may be either a bandpass band or the lowpass band.</dd>
      <dt><code>t</code></dt>
      <dd>The point in time the coefficient <code>c0</code> pertains to, in samples</dd>
-      <dt><code>t0</code></dt><dd>When not <code>INT64_MIN</code>, only apply <code>f</code> to the coefficients for time &ge; <code>t0</code></dd>
-      <dt><code>t1</code></dt><dd>When not <code>INT64_MAX</code>, only apply <code>f</code> to the coefficients for time &lt; <code>t1</code></dd>
    </dl>
    </dd>
+    <dt><code>t0</code></dt><dd>When not <code>INT64_MIN</code>, only apply <code>f</code> to the coefficients for time &ge; <code>t0</code></dd>
+    <dt><code>t1</code></dt><dd>When not <code>INT64_MAX</code>, only apply <code>f</code> to the coefficients for time &lt; <code>t1</code></dd>
 </dl>

 <div class="nav"><span class="prev"><a href="intro.html">Previous: API Introduction</a></span><span class="next"><a href="render_h.html">Next: Spectrogram rendering: <code>render.h</code></a></span></div>
--- a/doc/ref/intro.html
+++ b/doc/ref/intro.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2018-2020 Andreas Gustafsson.  This file is part of
+Copyright (C) 2018-2023 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.
 -->
@ -15,17 +15,17 @@ the top level of the distribution for license information.
 <p>The public API of the Gaborator library is defined in the HTML
 documentation in the form of annotated C++ declarations.  These are
 similar to the actual declarations in the respective header files, but
-simplified for clarity and omitting implementation details.</p>
+simplified for clarity, omitting implementation details.</p>

 <p>The actual implementation in the header file may be different in a
-number of ways but nonetheless compatible with the documented API.
+number of ways, but nonetheless compatible with the documented API.
 For example, classes may be declared using the keyword <code>struct</code>
 rather than <code>class</code>, function parameter names may be
 different, types may be declared using different but equivalent
 typedefs, and functions or templates in the header file may have
 additional arguments with default values.  Any classes, functions, and
 other definitions not mentioned in the documentation should be
-considered private and are subject to change or deletion without
+considered private and subject to change or deletion without
 notice.
 </p>

@ -35,6 +35,22 @@ with <code>gaborator::</code>, or use <code>using namespace
 gaborator;</code>.
 </p>

+<h2>Frequency and Time Units</h2>
+
+<p>The Gaborator is generally <i>sample rate agnostic</i>: it does not
+know or care what the actual sample rate of the signal is, and will
+work equally well for signals measured in milliherz or gigaherz.
+Frequencies are specified as dimensionless numbers, in units of the
+sample rate.  Since the highest frequency that can be represented by a
+sampled signal is half the sample rate, frequencies in the Gaborator
+range from 0 to 0.5.  In the API, arguments that represent frequencies
+in this way have names beginning with <code>ff_</code> (short
+for <i>fractional frequency</i>) as a reminder that the values are
+fractions of the sample rate rather than having a unit such as Hz.</p>
+
+<p>Similarly, time is measured not in an absolute unit such as
+seconds, but the dimensionless unit of samples.</p>
+
 <div class="nav"><span class="prev"><a href="../synth.html">Previous: Example 5: Synthesis from Scratch</a></span><span class="next"><a href="gaborator_h.html">Next: Spectrum analysis and synthesis: <code>gaborator.h</code></a></span></div>

 </body>
--- a/doc/ref/render_h.html
+++ b/doc/ref/render_h.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2018-2020 Andreas Gustafsson.  This file is part of
+Copyright (C) 2018-2023 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,9 +32,15 @@ down by powers of two.
  <dt><code>c</code></dt>
  <dd>A set of spectrogram coefficients to render</dd>
  <dt><code>xorigin</code></dt>
-  <dd>The point in time corresponding to pixel X coordinate 0, in samples</dd>
+  <dd>The point in time corresponding to pixel X coordinate 0, in samples.
+  Due to limitations of the implementation, this must currently be 0.
+  Instead, horizontal translation must be done by means of converting the
+  origin to pixel units and adding it to both <code>xi0</code> and <code>yi0</code>.</dd>
  <dt><code>yorigin</code></dt>
-  <dd>The band number of the frequency band corresponding to pixel Y coordinate 0</dd>
+  <dd>The band number of the frequency band corresponding to pixel Y coordinate 0.
+  Typically 0, or the band number of a reference frequency such as 440&nbsp;Hz,
+  to make that frequency always coincide exactly with a pixel regardless of the
+  scale factor.</dd>
  <dt><code>xi0</code></dt>
  <dd>The X coordinate of the first pixel to render</dd>
  <dt><code>xi1</code></dt>
--- a/doc/render.html
+++ b/doc/render.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2017-2020 Andreas Gustafsson.  This file is part of
+Copyright (C) 2017-2023 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.
 -->
@ -91,11 +91,13 @@ mix down the channels to mono, into a new <code>std::vector&lt;float&gt;</code>:
    }
 </pre>

-<h2>The Spectrum Analysis Parameters</h2>
+<h2>The Frequency Scale</h2>

-<p>Next, we need to choose some parameters for the spectrum analysis:
-the frequency resolution, the frequency range, and optionally a
-reference frequency.</p>
+<p>To make a constant-Q spectrogram, we will use a logarithmic
+frequency scale, represented by a <code>gaborator::log_fq_scale</code> object.
+To create a <code>log_fq_scale</code>, we must choose some
+parameters: the frequency resolution, the frequency range, and
+optionally a reference frequency.</p>

 <p>The frequency resolution is specified as a number of frequency
 bands per octave.  A typical number for analyzing music signals is 48
@ -123,10 +125,22 @@ typically 440 Hz, the standard tuning of the note <i>A<sub>4</sub></i>.
 Like the minimum frequency, it is given in
 units of the sample rate, so we pass <code>440.0 / fs</code>.</p>

-<p>The parameters are held in an object of type
-<code>gaborator::parameters</code>:
+<p>Let's construct a <code>log_fq_scale</code>:</p>
+
 <pre>
-    gaborator::parameters params(48, 20.0 / fs, 440.0 / fs);
+    gaborator::log_fq_scale scale(48, 20.0 / fs, 440.0 / fs);
+</pre>
+
+<h2>The Spectrum Analysis Parameters</h2>
+
+<p>The frequency scale forms a part of the spectrum analysis parameters,
+which are stored in an object of type <code>gaborator::parameters</code>.
+In this example, there is no need to specify parameters other than the
+frequency scale, so the parameters can be constructed simply by calling
+the <code>gaborator::parameters</code> constructor with the frequency scale
+as an argument:</p>
+<pre>
+    gaborator::parameters params(scale);
 </pre>

 <h2>The Spectrum Analyzer</h2>
@ -170,7 +184,7 @@ array, in units of samples.  Since we are passing the whole file at
 once, the beginning of the range is sample number zero, and the end is
 sample number <code>mono.size()</code>.  The fourth argument is a
 reference to the set of coefficients that the results of the spectrum
-analysis will be stored in.
+analysis will be stored in (or to be precise, added to).
 </p>
 <pre>
    analyzer.analyze(mono.data(), 0, mono.size(), coefs);
@ -323,16 +337,16 @@ we just need to finish <code>main()</code>:
 }
 </pre>

-<a name="compiling"><h2>Compiling</h2></a>
+<a id="compiling"><h2>Compiling</h2></a>
 <p>
 If you are using macOS, Linux, NetBSD, or a similar system, you can build
 the example by running the following command in the <code>examples</code>
 subdirectory.
-You need to have <i>libsndfile</i> is installed and supported by
+You need to have <i>libsndfile</i> installed and supported by
 <code>pkg-config</code>.
 </p>
-<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 class="build Darwin Linux NetBSD FreeBSD OpenBSD">
+c++ -std=c++11 -I.. -O3 -ffast-math $(pkg-config --cflags sndfile) render.cc $(pkg-config --libs sndfile) -o render
 </pre>

 <h2>Compiling for Speed</h2>
@ -344,34 +358,34 @@ can use the FFT from Apple's vDSP library by defining
 framework:
 </p>
 <pre class="build Darwin">
-c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP `pkg-config --cflags sndfile` render.cc `pkg-config --libs sndfile` -framework Accelerate -o render
+c++ -std=c++11 -I.. -O3 -ffast-math -DGABORATOR_USE_VDSP $(pkg-config --cflags sndfile) render.cc $(pkg-config --libs sndfile) -framework Accelerate -o render
 </pre>

 <p>On Linux and NetBSD, you can use the PFFFT (Pretty Fast FFT) library.
 You can get the latest version from
 <a href="https://bitbucket.org/jpommier/pffft">https://bitbucket.org/jpommier/pffft</a>,
-or the exact version that was used for testing from gaborator.com:
+or the exact version the example code was tested with from gaborator.com:
 </p>
 <!-- ftp https://bitbucket.org/jpommier/pffft/get/29e4f76ac53b.zip -->
-<pre class="build Linux NetBSD FreeBSD">
+<pre class="build Linux NetBSD FreeBSD OpenBSD">
 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 FreeBSD">
+<pre class="build Linux NetBSD FreeBSD OpenBSD">
 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 FreeBSD">
+<pre class="build Linux NetBSD FreeBSD OpenBSD">
 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 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 class="build Linux NetBSD FreeBSD OpenBSD">
+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>

 <h2>Running</h2>
--- a/doc/snr.html
+++ b/doc/snr.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2019-2020 Andreas Gustafsson.  This file is part of
+Copyright (C) 2019-2023 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.
 -->
@ -20,7 +20,7 @@ and comparing the resynthesis result to the original.
 </p>

 <p>Since it does not involve any audio file I/O, this example
-does not require the sndfile library, making it the shortest
+does not require the libsndfile library, making it the shortest
 and simplest one by far.</p>

 <h2>Preamble</h2>
@ -34,8 +34,8 @@ and simplest one by far.</p>

 <h2>Amplitude Measurement</h2>
 <p>To calculate the signal-to-noise ratio, we need to measure the
-amplitude of the orignal signal and the error residue.  We will use
-the root-mean-square amplitude, which is calculcated by the
+amplitude of the original signal and the error residue.  We will use
+the root-mean-square amplitude, which is calculated by the
 function <code>rms()</code>.
 </p>
 <pre>
@ -64,10 +64,10 @@ int main(int argc, char **argv) {
 <p>Then we create a spectrum analyzer with 48 bands per octave
 and a frequency range of 3 decades (0.0005 to 0.5 times the sample rate):</p>
 <pre>
-    gaborator::parameters params(48, 5e-4);
+    gaborator::parameters params(gaborator::log_fq_scale(48, 5e-4));
    gaborator::analyzer&lt;float&gt; analyzer(params);
 </pre>
-<p>...and run the spectrum analyzis:</p>
+<p>...and run the spectrum analysis:</p>
 <pre>
    gaborator::coefs&lt;float&gt; coefs(analyzer);
    analyzer.analyze(signal_in.data(), 0, len, coefs);
@ -96,16 +96,16 @@ 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 FreeBSD">
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` snr.cc `pkg-config --libs sndfile` -o snr
+<pre class="build Darwin Linux NetBSD FreeBSD OpenBSD">
+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>
 <pre class="build Darwin">
-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
+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">
-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
+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>

 <h2>Running</h2>
@ -113,7 +113,7 @@ c++ -std=c++11 -I.. -Ipffft -O3 -ffast-math -DGABORATOR_USE_PFFFT `pkg-config --
 <pre class="run">
 ./snr
 </pre>
-<p>This will print the SNR which should be more than 100 dB if the library is working correctly.</p>
+<p>This will print the SNR which should be more than 100 dB when the library is working correctly.</p>

 <div class="nav"><span class="prev"><a href="stream.html">Previous: Example 3: Streaming</a></span><span class="next"><a href="synth.html">Next: Example 5: Synthesis from Scratch</a></span></div>

--- a/doc/stream.html
+++ b/doc/stream.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2018-2020 Andreas Gustafsson.  This file is part of
+Copyright (C) 2018-2023 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.
 -->
@ -18,7 +18,7 @@ the top level of the distribution for license information.
 rather than operating on a complete recording at once as in the previous
 examples.</p>

-<p>This program doesn't do anything particulary useful &mdash; it just
+<p>This program doesn't do anything particularly useful &mdash; it just
 inverts the phase of the signal, but not using the obvious method of
 changing the sign of each sample but by changing the sign of each
 spectrogram coefficient.  Consider the expression <code>coef =
@ -87,7 +87,8 @@ to reduce the latency, the number of frequency bands per octave is reduced
 from 48 to 12 (one per semitone), and the lower frequency limit of
 the bandpass filter bank is raised from 20&nbsp;Hz to 200&nbsp;Hz.</p>
 <pre>
-    gaborator::parameters params(12, 200.0 / fs, 440.0 / fs);
+    gaborator::log_fq_scale scale(12, 200.0 / fs, 440.0 / fs);
+    gaborator::parameters params(scale);
    gaborator::analyzer&lt;float&gt; analyzer(params);
 </pre>

@ -185,7 +186,7 @@ the last time and will not be affected by future input blocks.
 Therefore, it is now safe to examine and/or modify these
 coefficients as required by your application.  Here, by way
 of example, we simply change their signs to invert the phase of the signal.
-Note that unlike the earlier filter example where <code>prorcess()</code>
+Note that unlike the earlier filter example where <code>process()</code>
 applied a function to all the coefficients, here it is applied only to
 the coefficients within a limited time range.
 </p>
@ -225,7 +226,7 @@ blocksize</code>.</p>
 </pre>
 <p>Coefficients older than <code>t_out + blocksize - synthesis_support</code>
 will no longer be needed to synthesize the next block of output signal, so
-it's now OK to forget them and free the memory they used:
+we can now forget them to free up the memory they used:
 </p>
 <pre>
        forget_before(analyzer, coefs, t_out + blocksize - synthesis_support);
@ -247,16 +248,16 @@ 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 FreeBSD">
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` stream.cc `pkg-config --libs sndfile` -o stream
+<pre class="build Darwin Linux NetBSD FreeBSD OpenBSD">
+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>
 <pre class="build Darwin">
-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
+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">
-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
+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>
--- a/doc/synth.html
+++ b/doc/synth.html
@ -1,6 +1,6 @@
 <!DOCTYPE html>
 <!--
-Copyright (C) 2020 Andreas Gustafsson.  This file is part of
+Copyright (C) 2020-2023 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.
 -->
@ -43,7 +43,7 @@ int main(int argc, char **argv) {
 <p>Although this example does not perform any analysis, we nonetheless
 need to create an <code>analyzer</code> object, as it is used for both
 analysis and synthesis purposes.  To generate the frequencies of the
-12-note equal-tempered scale, we need 12 bands per octave; a multiple
+12-note equal-tempered scale, we need 12 bands per octave.  A multiple
 of 12 would also work, but here we don't need the added frequency
 resolution that would bring, and the time resolution would be
 worse.</p>
@ -55,7 +55,22 @@ bandpass filter bank, but that doesn't matter.</p>

 <pre>
    double fs = 44100;
-    gaborator::parameters params(12, 20.0 / fs, 8.18 / fs);
+    gaborator::log_fq_scale scale(12, 20.0 / fs, 8.18 / fs);
+    gaborator::parameters params(scale);
+</pre>
+
+<p>As we create new complex coefficients from scratch, we need to be
+careful about their phases.  For any given frequency, the portions of
+the output signal contributed by the coefficients at different points
+in time need to be in phase so that they combine constructively into a
+single (co)sine wave rather than interfering destructively.  The
+easiest way to achieve this is to select the global phase convention
+in the synthesis parameters.  With this setting, the coefficient phases
+will be interpreted relative to the common reference point of t=0
+rather than the time of each individual coefficient:</p>
+
+<pre>
+    params.phase = gaborator::coef_phase::global;
    gaborator::analyzer&lt;float&gt; analyzer(params);
 </pre>

@ -136,9 +151,9 @@ in front of <code>midi_note</code>.

 <p>We can now synthesize audio from the coefficients by
 calling <code>synthesize()</code>.  Audio will be generated
-starting half a second before the first note to allow for the pre-ringing
-of the synthesis filter, and ending a few seconds after the
-last note to allow for its decay.
+starting half a second before the first note to allow for pre-ringing
+of the synthesis filters, and ending a few seconds after the
+last note to give the note time to decay.
 </p>
 <pre>
    double audio_start_time = -0.5;
@ -188,16 +203,16 @@ 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 FreeBSD">
-c++ -std=c++11 -I.. -O3 -ffast-math `pkg-config --cflags sndfile` synth.cc `pkg-config --libs sndfile` -o synth
+<pre class="build Darwin Linux NetBSD FreeBSD OpenBSD">
+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>
 <pre class="build Darwin">
-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
+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">
-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
+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>

 <h2>Running</h2>
--- a/examples/filter.cc
+++ b/examples/filter.cc
@ -28,7 +28,8 @@ int main(int argc, char **argv) {
        exit(1);
    }
    sf_close(sf_in);
-    gaborator::parameters params(100, 20.0 / fs);
+    gaborator::log_fq_scale scale(10, 20.0 / fs);
+    gaborator::parameters params(scale);
    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++) {
--- a/examples/render.cc
+++ b/examples/render.cc
@ -36,7 +36,8 @@ int main(int argc, char **argv) {
            v += audio[i * sfinfo.channels + c];
        mono[i] = v;
    }
-    gaborator::parameters params(48, 20.0 / fs, 440.0 / fs);
+    gaborator::log_fq_scale scale(48, 20.0 / fs, 440.0 / fs);
+    gaborator::parameters params(scale);
    gaborator::analyzer<float> analyzer(params);
    gaborator::coefs<float> coefs(analyzer);
    analyzer.analyze(mono.data(), 0, mono.size(), coefs);
--- a/examples/snr.cc
+++ b/examples/snr.cc
@ -18,7 +18,7 @@ int main(int argc, char **argv) {
    std::uniform_real_distribution<> uniform(-1.0, 1.0);
    for (size_t i = 0; i < len; i++)
        signal_in[i] = uniform(rand);
-    gaborator::parameters params(48, 5e-4);
+    gaborator::parameters params(gaborator::log_fq_scale(48, 5e-4));
    gaborator::analyzer<float> analyzer(params);
    gaborator::coefs<float> coefs(analyzer);
    analyzer.analyze(signal_in.data(), 0, len, coefs);
--- a/examples/stream.cc
+++ b/examples/stream.cc
@ -32,7 +32,8 @@ int main(int argc, char **argv) {
    }
    sf_command(sf_in, SFC_SET_NORM_FLOAT, NULL, SF_FALSE);
    sf_command(sf_out, SFC_SET_NORM_FLOAT, NULL, SF_FALSE);
-    gaborator::parameters params(12, 200.0 / fs, 440.0 / fs);
+    gaborator::log_fq_scale scale(12, 200.0 / fs, 440.0 / fs);
+    gaborator::parameters params(scale);
    gaborator::analyzer<float> analyzer(params);
    size_t analysis_support = ceil(analyzer.analysis_support());
    size_t synthesis_support = ceil(analyzer.synthesis_support());
--- a/examples/synth.cc
+++ b/examples/synth.cc
@ -11,7 +11,9 @@ int main(int argc, char **argv) {
        exit(1);
    }
    double fs = 44100;
-    gaborator::parameters params(12, 20.0 / fs, 8.18 / fs);
+    gaborator::log_fq_scale scale(12, 20.0 / fs, 8.18 / fs);
+    gaborator::parameters params(scale);
+    params.phase = gaborator::coef_phase::global;
    gaborator::analyzer<float> analyzer(params);
    static int pentatonic[] = { 57, 60, 62, 64, 67 };
    int n_notes = 64;
--- a/gaborator/affine_transform.h
+++ b/gaborator/affine_transform.h
@ -23,7 +23,7 @@ struct affine_transform {
    double a, b;
 };

-// Composition
+// Composition: (a * b)(x) = a(b(x))

 static inline affine_transform
 operator *(const affine_transform &a, const affine_transform &b) {
--- a/gaborator/fft_naive.h
+++ b/gaborator/fft_naive.h
@ -1,7 +1,7 @@
 //
 // Fast Fourier transform, naive reference implementations
 //
-// Copyright (C) 1992-2021 Andreas Gustafsson.  This file is part of
+// Copyright (C) 1992-2022 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.
 //
@ -21,6 +21,14 @@

 namespace gaborator {

+// The template argument I is the iterator type used for accessing the
+// input and output data.  Typically this is "std::complex<float> *" or
+// "std::complex<double> *", but other types such as stride iterators
+// can also be used.
+//
+// This is not fully const correct as the same iterator type is used
+// for both read-only and read-write access.
+
 template <class I>
 struct fft {
    typedef typename std::iterator_traits<I>::value_type C; // complex
--- a/gaborator/gaborator.h
+++ b/gaborator/gaborator.h
--- a/gaborator/gaussian.h
+++ b/gaborator/gaussian.h
@ -1,7 +1,7 @@
 //
 // The Gaussian and related functions
 //
-// Copyright (C) 2015-2021 Andreas Gustafsson.  This file is part of
+// Copyright (C) 2015-2023 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,13 +17,13 @@ namespace gaborator {
 // error function.  This is good for arguments in the range 1e-8 to 1,
 // to within a few percent.

-inline float erfc_inv(float x) {
-    return sqrtf(-logf(x)) - 0.3f;
+static inline double erfc_inv(double x) {
+    return sqrt(-log(x)) - 0.25;
 }

 // Gaussian with peak = 1

-inline double norm_gaussian(double sd, double x) {
+static inline double norm_gaussian(double sd, double x) {
    return exp(-(x * x) / (2 * sd * sd));
 }

@ -100,22 +100,22 @@ double gaussian_value_support_inv(double support, double max_error) {
 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
--- a/gaborator/pod_vector.h
+++ b/gaborator/pod_vector.h
@ -1,7 +1,7 @@
 //
 // A vector class without default-initialization, for "plain old data"
 //
-// Copyright (C) 2016-2021 Andreas Gustafsson.  This file is part of
+// Copyright (C) 2016-2022 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.
 //
@ -28,9 +28,7 @@ struct pod_vector {
        b = e = 0;
    }
    explicit pod_vector(size_t size_) {
-        // Allocate raw uninitialized memory
-        b = static_cast<T *>(::operator new(size_ * sizeof(T)));
-        e = b + size_;
+        allocate(size_);
    }
    ~pod_vector()
 #if __cplusplus >= 201103L
@ -60,11 +58,8 @@ struct pod_vector {
        b = n;
        e = n + new_size;
    }
-    pod_vector(const pod_vector &a)
-
-    {
-        b = new T[a.size()];
-        e = b + a.size();
+    pod_vector(const pod_vector &a) {
+        allocate(a.size());
        std::copy(a.b, a.e, b);
        //if (size()) fprintf(stderr, "pod_vector cc %d\n", (int)size());
    }
@ -102,6 +97,11 @@ struct pod_vector {
    }
 #endif
 private:
+    void allocate(size_t n) {
+        // Allocate raw uninitialized memory
+        b = static_cast<T *>(::operator new(n * sizeof(T)));
+        e = b + n;
+    }
    void _free() {
        // Free as raw uninitialized memory
        ::operator delete(b);
--- a/gaborator/pool.h
+++ b/gaborator/pool.h
@ -1,7 +1,7 @@
 //
 // Pool of shared objects
 //
-// Copyright (C) 2015-2018 Andreas Gustafsson.  This file is part of
+// Copyright (C) 2015-2022 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.
 //
@ -15,8 +15,8 @@ namespace gaborator {

 // The "pool" class is for sharing FFT objects so that we don't
 // create multiple FFTs of the same size.  It could also be used
-// to share objects of some other class T where we don't want to
-// create multiple Ts with the same K.
+// to share objects of some other class T where we wish to avoid
+// creating multiple Ts with the same K.

 template <class T, class K>
 struct pool {
@ -27,7 +27,8 @@ struct pool {
        }
    }
    T *get(const K &k) {
-        std::pair<typename m_t::iterator, bool> r = m.insert(std::make_pair(k, (T *)0));
+        std::pair<typename m_t::iterator, bool> r =
+            m.insert(std::make_pair(k, (T *)0));
        if (r.second) {
            // New element was inserted
            assert((*(r.first)).second == 0);
@ -36,12 +37,8 @@ struct pool {
        return (*r.first).second;
    }
    m_t m;
-    static pool shared;
 };

-template <class T, class K>
-pool<T, K> pool<T, K>::shared;
-
 } // namespace

 #endif
--- a/gaborator/ref.h
+++ b/gaborator/ref.h
@ -1,7 +1,7 @@
 //
 // Intrusive reference counting smart pointer
 //
-// Copyright (C) 2016-2018 Andreas Gustafsson.  This file is part of
+// Copyright (C) 2016-2022 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.
 //
@ -24,19 +24,20 @@ struct refcounted {
 // object type and invoke operator delete on the base class.

 template <class T>
-void incref(T &r) {
-    r.refcount++;
+void incref(T *r) {
+    r->refcount++;
 }

 template <class T>
-void decref(T &r) {
-    r.refcount--;
-    if (r.refcount == 0)
-        delete &r;
+void decref(T *r) {
+    r->refcount--;
+    if (r->refcount == 0)
+        delete r;
 }

 template <class T>
 struct ref {
+    typedef T element_type;
    ref(): p(0) { }
    ref(T *p_): p(p_) {
        _incref();
@ -44,7 +45,15 @@ struct ref {
    ref(const ref &o): p(o.p) {
        _incref();
    }
-    ref &operator=(const ref &o) { reset(o.p); return *this; }
+    // Move constructor
+    ref(ref &&o) {
+        p = o.p;
+        o.p = 0;
+    }
+    ref &operator=(const ref &o) {
+        reset(o.p);
+        return *this;
+    }
    ~ref() { reset(); }
    void reset() {
        _decref();
@ -65,12 +74,12 @@ private:
    void _incref() {
        if (! p)
            return;
-        incref(*p);
+        incref(p);
    }
    void _decref() {
        if (! p)
            return;
-        decref(*p);
+        decref(p);
    }
    T *p;
 };
--- a/gaborator/render.h
+++ b/gaborator/render.h
@ -1,7 +1,7 @@
 //
 // Rendering of spectrogram images
 //
-// Copyright (C) 2015-2021 Andreas Gustafsson.  This file is part of
+// Copyright (C) 2015-2023 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 @@ unsigned int float2pixel_8bit(T val) {
    // of a negative number; those can arise from bicubic
    // interpolation.  While we're at it, let's also skip the gamma
    // correction for small numbers that will round to zero anyway,
-    // and especially denormals which could rigger GCC bug target/83240.
+    // and especially denormals which could trigger GCC bug target/83240.
    static const T almost_zero = 1.0 / 65536;
    if (val < almost_zero)
        val = 0;
@ -43,7 +43,7 @@ unsigned int float2pixel_8bit(T val) {

 // Magnitude

-template<class T>
+template <class T>
 struct complex_abs_fob {
    T operator()(const complex<T> &c) {
        return complex_abs(c);
@ -54,9 +54,16 @@ struct complex_abs_fob {

 // T -> f() -> OI::value_type

+
 template <class F, class OI, class T>
-struct transform_output_iterator: public std::iterator<std::output_iterator_tag, T> {
+struct transform_output_iterator {
+    // Iterator traits
+    typedef std::output_iterator_tag iterator_category;
    typedef T value_type;
+    typedef std::ptrdiff_t difference_type;
+    typedef T* pointer;
+    typedef T& reference;
+
    transform_output_iterator(F f_, OI output_): f(f_), output(output_) { }
    transform_output_iterator<F, OI, T>& operator=(T v) {
        *output++ = f(v);
@ -147,7 +154,7 @@ OI render_line(const analyzer<T> &anl,

    int oct;
    unsigned int obno; // Band number within octave
-    bool clip = ! bno_split(*msc.meta, gbno, oct, obno, false);
+    bool clip = ! bno_split(*msc.cmeta, gbno, oct, obno, false);
    if (clip) {
        for (sample_index_t i = i0; i < i1; i++)
            *output++ = (T)0;
@ -157,8 +164,8 @@ OI render_line(const analyzer<T> &anl,
        abs_rowsource(msc, oct, obno, normf);

    // 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));
+    int scale_exp = band_step_log2(*msc.cmeta, gbno);
+    RESAMPLER x_resampler(zoom_p2_src(xf, -scale_exp));
    output = x_resampler(abs_rowsource, i0, i1, output);
    return output;
 }
@ -228,10 +235,15 @@ struct transverse_source {
 };

 template <class I>
-struct stride_iterator: public std::iterator
-    <std::forward_iterator_tag, typename std::iterator_traits<I>::value_type>
-{
+struct stride_iterator {
    typedef typename std::iterator_traits<I>::value_type T;
+    // Iterator traits
+    typedef std::forward_iterator_tag iterator_category;
+    typedef T value_type;
+    typedef std::ptrdiff_t difference_type;
+    typedef T* pointer;
+    typedef T& reference;
+
    stride_iterator(I it_, size_t stride_): it(it_), stride(stride_) { }
    T& operator*() { return *it; }
    stride_iterator<I>& operator++() {
@ -271,7 +283,7 @@ struct updated_nop {
 // are updated, the time range being from inc_i0 (inclusive) to inc_i1
 // (exclusive).  The updated parts of the image consist of zero or
 // more non-overlapping rectangles; to find out what those are, pass a
-// function "update" which will be called will each rectangle in turn
+// function "update" which will be called with each rectangle in turn
 // after it has been updated.

 // For non-incremental rendering, pass inc_i0 = INT64_MIN and inc_i1 =
@ -299,7 +311,7 @@ void render_incremental(
    assert(yi1 >= yi0);
    assert(inc_i1 >= inc_i0);

-    // The data type to reasmple
+    // The data type to resample
    typedef typename NORMF::return_type RST;

    // Vertical resampler
@ -328,12 +340,12 @@ void render_incremental(
        y_resampler.support(y, y, ysi0, ysi1);
        int64_t band = ysi1;
        // Clamp the band to the valid range
-        band = std::max(band, (int64_t)anl.bandpass_bands_begin());
-        band = std::min(band, (int64_t)anl.bandpass_bands_end() - 1);
+        set_max(band, (int64_t)anl.bandpass_bands_begin());
+        set_min(band, (int64_t)anl.bandpass_bands_end() - 1);

        // Find the analysis support in the time (x) dimension,
        // in signal samples
-        double support = anl.analysis_support(band);
+        double support = anl.band_analysis_support(band);
        // Convert from signal samples to coefficient samples
        int scale_exp = anl.band_scale_exp((int)band);

@ -341,7 +353,7 @@ void render_incremental(
        // support, and map it back to pixel space to find the
        // affected pixel range, taking the resampler support
        // into account.
-        RESAMPLER x_resampler(zoom_p2(x_xf, -scale_exp));
+        RESAMPLER x_resampler(zoom_p2_src(x_xf, -scale_exp));
        int64_t ceil_support = ceil(support);

        // inv_support() calculates both sides of the support at once,
@ -360,17 +372,17 @@ void render_incremental(
        if (inc_i0 == INT64_MIN) {
            adj_x0 = xi0;
        } else {
-            adj_x0 = std::max(xi0, adj_x0);
+            set_max(adj_x0, xi0);
            // Don't let width go negative
-            adj_x0 = std::min(adj_x0, xi1);
+            set_min(adj_x0, xi1);
        }

        if (inc_i1 == INT64_MAX) {
            adj_x1 = xi1;
        } else {
-            adj_x1 = std::min(xi1, adj_x1);
+            set_min(adj_x1, xi1);
            // Don't let width go negative
-            adj_x1 = std::max(adj_x1, adj_x0);
+            set_max(adj_x1, adj_x0);
        }

        assert(adj_x0 <= adj_x1);
@ -423,7 +435,7 @@ void render_incremental(
    // the support also varies; use the largest resampling factor of
    // any band to get the worst-case support.
    int worstcase_band = anl.bandpass_bands_end() - 1;
-    RESAMPLER x_resampler(zoom_p2(x_xf, -anl.band_scale_exp(worstcase_band)));
+    RESAMPLER x_resampler(zoom_p2_src(x_xf, -anl.band_scale_exp(worstcase_band)));
    int64_t xsi0, xsi1;
    x_resampler.support(xi0, xi1, xsi0, xsi1);

--- a/gaborator/resample2.h
+++ b/gaborator/resample2.h
@ -55,7 +55,10 @@ namespace gaborator {

 struct p2_transform {
    p2_transform(int e_, int64_t origin_): e(e_), origin(origin_) { }
-    // Convert a linear transform into a p2_transform
+    // Convert an affine transform into a p2_transform.
+    // The affine transform must be representable as a
+    // p2_transform, meaning xf.a must be a power of two,
+    // and the origin must be an integer.
    p2_transform(affine_transform xf) {
        int exp;
        double m = frexp(xf.a, &exp);
@ -68,18 +71,22 @@ struct p2_transform {
    int64_t origin;
 };

-// Scale a transform by a power of two
-
-static inline p2_transform
-zoom_p2(p2_transform xf, int e) {
-    return p2_transform(xf.e + e, xf.origin);
-}
+// Scale a transform by a power of two, around the origin
+// of the destination space

 static inline affine_transform
-zoom_p2(affine_transform xf, int e) {
+zoom_p2_dst(affine_transform xf, int e) {
    return affine_transform(ldexp(xf.a, e), xf.b);
 }

+// Scale a transform by a power of two, around the origin
+// of the source space
+
+static inline affine_transform
+zoom_p2_src(affine_transform xf, int e) {
+    return affine_transform(ldexp(xf.a, e), ldexp(xf.b, e));
+}
+
 // Resample data from "source", generating a view between indices
 // i0 and i1 of the scale determined by exponent e, and storing
 // i1 - i0 samples starting at *out.
@ -211,11 +218,11 @@ OI resample2_p2xf(const S &source, p2_transform xf,
 // As above, but taking an affine_transform.

 template <class S, class OI>
-OI resample2(const S &source, affine_transform lxf,
+OI resample2(const S &source, affine_transform axf,
             int64_t i0, int64_t i1,
             bool interpolate, OI out)
 {
-    p2_transform xf(lxf);
+    p2_transform xf(axf);
    typedef typename std::iterator_traits<OI>::value_type T;
    gaborator::pod_vector<T> data(i1 - i0);
    T *p = data.data();
@ -231,22 +238,22 @@ OI resample2(const S &source, affine_transform lxf,
 // return an unnecessarily large support when interpolation is off

 inline void
-resample2_support(affine_transform lxf, int64_t i0, int64_t i1,
+resample2_support(affine_transform axf, int64_t i0, int64_t i1,
                  int64_t &si0_ret, int64_t &si1_ret)
 {
-    p2_transform xf(lxf);
+    p2_transform xf(axf);
    int margin = 2;
    if (xf.e > 0) {
        // Note code duplication wrt resample2_p2xf().
        // Also note tail recursion.
        int64_t si0 = i0 * 2 - margin + 1;
        int64_t si1 = i1 * 2 + margin;
-        resample2_support(zoom_p2(lxf, -1),
+        resample2_support(zoom_p2_dst(axf, -1),
                          si0, si1, si0_ret, si1_ret);
    } else if (xf.e < 0) {
        int64_t si0 = (i0 >> 1) - margin;
        int64_t si1 = ((i1 - 1) >> 1) + margin + 1;
-        resample2_support(zoom_p2(lxf, +1),
+        resample2_support(zoom_p2_dst(axf, +1),
                          si0, si1, si0_ret, si1_ret);
    } else {
        si0_ret = xf.origin + i0;
@ -258,21 +265,21 @@ resample2_support(affine_transform lxf, int64_t i0, int64_t i1,
 // destination indices that depend on a given range of source indices.

 inline void
-resample2_inv_support(affine_transform lxf, int64_t si0, int64_t si1,
+resample2_inv_support(affine_transform axf, int64_t si0, int64_t si1,
                      int64_t &i0_ret, int64_t &i1_ret)
 {
-    p2_transform xf(lxf);
+    p2_transform xf(axf);
    // Conservative
    int margin = 2;
    if (xf.e > 0) {
        int64_t di0, di1;
-        resample2_inv_support(zoom_p2(lxf, -1),
+        resample2_inv_support(zoom_p2_dst(axf, -1),
                              si0, si1, di0, di1);
        i0_ret = di0 >> 1;
        i1_ret = (di1 >> 1) + 1;
    } else if (xf.e < 0) {
        int64_t di0, di1;
-        resample2_inv_support(zoom_p2(lxf, +1),
+        resample2_inv_support(zoom_p2_dst(axf, +1),
                              si0, si1, di0, di1);
        i0_ret = di0 * 2 - margin + 1;
        i1_ret = di1 * 2 + margin;
--- a/gaborator/vector_math.h
+++ b/gaborator/vector_math.h
@ -1,7 +1,7 @@
 //
 // Vector math operations
 //
-// Copyright (C) 2016-2018 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,6 +17,13 @@

 #include <complex>

+// Use __restrict__ where supported
+#if defined(__GNUC__) || defined(__clang__)
+#define GABORATOR_RESTRICT __restrict__
+#else
+#define GABORATOR_RESTRICT
+#endif
+
 namespace gaborator {

 // The _naive versions are used when SSE is not available, and as
@ -40,12 +47,12 @@ static inline
 std::complex<double> complex_mul(std::complex<double> a_,
                                 std::complex<double> b_)
 {
-    __v2df a = _mm_setr_pd(a_.real(), a_.imag());
-    __v2df b = _mm_setr_pd(b_.real(), b_.imag());
-    __v2df as = (__v2df) _mm_shuffle_pd(a, a, 0x1);
-    __v2df t0 = _mm_mul_pd(a, _mm_shuffle_pd(b, b, 0x0));
-    __v2df t1 = _mm_mul_pd(as, _mm_shuffle_pd(b, b, 0x3));
-    __v2df c = __builtin_ia32_addsubpd(t0, t1); // SSE3
+    __m128d a = _mm_setr_pd(a_.real(), a_.imag());
+    __m128d b = _mm_setr_pd(b_.real(), b_.imag());
+    __m128d as = (__m128d) _mm_shuffle_pd(a, a, 0x1);
+    __m128d t0 = _mm_mul_pd(a, _mm_shuffle_pd(b, b, 0x0));
+    __m128d t1 = _mm_mul_pd(as, _mm_shuffle_pd(b, b, 0x3));
+    __m128d c = _mm_addsub_pd(t0, t1); // SSE3
    return std::complex<double>(c[0], c[1]);
 }
 #else
@ -66,9 +73,9 @@ std::complex<float> complex_mul(std::complex<float> a_,

 template <class T, class U, class V>
 static inline void
-elementwise_product_naive(T *r,
-                          U *a,
-                          V *b,
+elementwise_product_naive(T *GABORATOR_RESTRICT r,
+                          U *GABORATOR_RESTRICT a,
+                          V *GABORATOR_RESTRICT b,
                          size_t n)
 {
    for (size_t i = 0; i < n; i++)
@ -77,9 +84,9 @@ elementwise_product_naive(T *r,

 template <class T, class U, class V, class S>
 static inline void
-elementwise_product_times_scalar_naive(T *r,
-                                       U *a,
-                                       V *b,
+elementwise_product_times_scalar_naive(T *GABORATOR_RESTRICT r,
+                                       U *GABORATOR_RESTRICT a,
+                                       V *GABORATOR_RESTRICT b,
                                       S s,
                                       size_t n)
 {
@ -90,12 +97,13 @@ elementwise_product_times_scalar_naive(T *r,
 // I is the input complex data type, O is the output data type
 template <class I, class O>
 static inline void
-complex_magnitude_naive(I *inv,
-                        O *outv,
+complex_magnitude_naive(I *GABORATOR_RESTRICT inv,
+                        O *GABORATOR_RESTRICT outv,
                        size_t n)
 {
    for (size_t i = 0; i < n; i++)
-        outv[i] = std::sqrt(inv[i].real() * inv[i].real() + inv[i].imag() * inv[i].imag());
+        outv[i] = std::sqrt(inv[i].real() * inv[i].real() +
+                            inv[i].imag() * inv[i].imag());
 }

 #if GABORATOR_USE_SSE3_INTRINSICS
@ -105,11 +113,11 @@ complex_magnitude_naive(I *inv,
 // Perform two complex float multiplies in parallel

 static inline
-__v4sf complex_mul_vec2(__v4sf aa, __v4sf bb) {
-    __v4sf aas =_mm_shuffle_ps(aa, aa, 0xb1);
-    __v4sf t0 = _mm_mul_ps(aa, _mm_moveldup_ps(bb));
-    __v4sf t1 = _mm_mul_ps(aas, _mm_movehdup_ps(bb));
-    return __builtin_ia32_addsubps(t0, t1); // SSE3
+__m128 complex_mul_vec2(__m128 aa, __m128 bb) {
+    __m128 aas =_mm_shuffle_ps(aa, aa, 0xb1);
+    __m128 t0 = _mm_mul_ps(aa, _mm_moveldup_ps(bb));
+    __m128 t1 = _mm_mul_ps(aas, _mm_movehdup_ps(bb));
+    return _mm_addsub_ps(t0, t1); // SSE3
 }

 // Calculate the elementwise product of a complex float vector
@ -123,12 +131,12 @@ elementwise_product(std::complex<float> *cv,
 {
    assert((n & 1) == 0);
    n >>= 1;
-    __v4sf *c = (__v4sf *) cv;
-    const __v4sf *a = (const __v4sf *) av;
-    const __v4sf *b = (const __v4sf *) bv;
+    __m128 *c = (__m128 *) cv;
+    const __m128 *a = (const __m128 *) av;
+    const __m128 *b = (const __m128 *) bv;
    while (n--) {
-        __v4sf aa = *a++;
-        __v4sf bb = *b++;
+        __m128 aa = *a++;
+        __m128 bb = *b++;
         *c++ = complex_mul_vec2(aa, bb);
    }
 }
@ -146,13 +154,13 @@ elementwise_product(std::complex<float> *cv,
 {
    assert((n & 3) == 0);
    n >>= 2;
-    __v4sf *c = (__v4sf *) cv;
-    const __v4sf *a = (const __v4sf *) av;
-    const __v4sf *b = (const __v4sf *) bv;
+    __m128 *c = (__m128 *) cv;
+    const __m128 *a = (const __m128 *) av;
+    const __m128 *b = (const __m128 *) bv;
    while (n--) {
-        __v4sf a0 = (__v4sf) _mm_loadu_si128((const __m128i *) a++);
-        __v4sf a1 = (__v4sf) _mm_loadu_si128((const __m128i *) a++);
-        __v4sf bb = *b++;
+        __m128 a0 = (__m128) _mm_loadu_si128((const __m128i *) a++);
+        __m128 a1 = (__m128) _mm_loadu_si128((const __m128i *) a++);
+        __m128 bb = *b++;
         *c++ = _mm_mul_ps(a0, _mm_unpacklo_ps(bb, bb));
         *c++ = _mm_mul_ps(a1, _mm_unpackhi_ps(bb, bb));
    }
@ -167,13 +175,13 @@ elementwise_product_times_scalar(std::complex<float> *cv,
 {
    assert((n & 1) == 0);
    n >>= 1;
-    const __v4sf *a = (const __v4sf *) av;
-    const __v4sf *b = (const __v4sf *) bv;
-    const __v4sf dd = (__v4sf) { d.real(), d.imag(), d.real(), d.imag() };
-    __v4sf *c = (__v4sf *) cv;
+    const __m128 *a = (const __m128 *) av;
+    const __m128 *b = (const __m128 *) bv;
+    const __m128 dd = (__m128) { d.real(), d.imag(), d.real(), d.imag() };
+    __m128 *c = (__m128 *) cv;
    while (n--) {
-        __v4sf aa = *a++;
-        __v4sf bb = *b++;
+        __m128 aa = *a++;
+        __m128 bb = *b++;
         *c++ = complex_mul_vec2(complex_mul_vec2(aa, bb), dd);
    }
 }
@ -192,19 +200,19 @@ complex_magnitude(std::complex<float> *inv,
        *outv++ = std::sqrt(v.real() * v.real() + v.imag() * v.imag());
        n--;
    }
-    const __v4sf *in = (const __v4sf *) inv;
-    __v4sf *out = (__v4sf *) outv;
+    const __m128 *in = (const __m128 *) inv;
+    __m128 *out = (__m128 *) outv;
    while (n >= 4) {
-        __v4sf aa = *in++; // c0re c0im c1re c1im
-        __v4sf aa2 = _mm_mul_ps(aa, aa); // c0re^2 c0im^2 c1re^2 c1im^2
-        __v4sf bb = *in++; // c2re c2im c3re c3im
-        __v4sf bb2 = _mm_mul_ps(bb, bb); // etc
+        __m128 aa = *in++; // c0re c0im c1re c1im
+        __m128 aa2 = _mm_mul_ps(aa, aa); // c0re^2 c0im^2 c1re^2 c1im^2
+        __m128 bb = *in++; // c2re c2im c3re c3im
+        __m128 bb2 = _mm_mul_ps(bb, bb); // etc
        // Gather the real parts: x0 x2 y0 y2
        // 10 00 10 00 = 0x88
-        __v4sf re2 =_mm_shuffle_ps(aa2, bb2, 0x88);
-        __v4sf im2 =_mm_shuffle_ps(aa2, bb2, 0xdd);
-        __v4sf mag2 = _mm_add_ps(re2, im2);
-        __v4sf mag = __builtin_ia32_sqrtps(mag2);
+        __m128 re2 =_mm_shuffle_ps(aa2, bb2, 0x88);
+        __m128 im2 =_mm_shuffle_ps(aa2, bb2, 0xdd);
+        __m128 mag2 = _mm_add_ps(re2, im2);
+        __m128 mag = _mm_sqrt_ps(mag2);
        // Unaligned store
        _mm_storeu_si128((__m128i *)out, (__m128i)mag);
        out++;
@ -298,4 +306,6 @@ complex_magnitude(I *inv,

 } // namespace

+#undef GABORATOR_RESTRICT
+
 #endif
--- a/gaborator/version.h
+++ b/gaborator/version.h
@ -1,7 +1,7 @@
 //
 // Version number
 //
-// Copyright (C) 2015-2021 Andreas Gustafsson.  This file is part of
+// Copyright (C) 2015-2023 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.
 //
@ -9,7 +9,7 @@
 #ifndef _GABORATOR_VERSION_H
 #define _GABORATOR_VERSION_H

-#define GABORATOR_VERSION_MAJOR 1
-#define GABORATOR_VERSION_MINOR 7
+#define GABORATOR_VERSION_MAJOR 2
+#define GABORATOR_VERSION_MINOR 0

 #endif
Author	SHA1	Message	Date
DataHoarder	f86f6d213a	Use standard intrinsic types, adjust __builtin with intrinsics	2023-10-15 16:43:27 +02:00
DataHoarder	cfe987a51d	Release 2.0	2023-10-15 16:42:31 +02:00