summaryrefslogtreecommitdiffstatshomepage
path: root/src/emu/resampler.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/emu/resampler.cpp')
-rw-r--r--src/emu/resampler.cpp317
1 files changed, 317 insertions, 0 deletions
diff --git a/src/emu/resampler.cpp b/src/emu/resampler.cpp
new file mode 100644
index 00000000000..b95fe80dfba
--- /dev/null
+++ b/src/emu/resampler.cpp
@@ -0,0 +1,317 @@
+// license:BSD-3-Clause
+// copyright-holders:Olivier Galibert
+
+// Audio resampler
+
+#include "emu.h"
+#include "resampler.h"
+
+// How an accurate resampler works ?
+
+// Resampling uses a number of well-known theorems we are not trying
+// to prove here.
+
+// Samping theorem. A digital signal sampled at frequency fs is
+// equivalent to an analog signal where all frequencies are between 0
+// and fs/2. Equivalent here means that the samples are unique given
+// the analog signal, the analog sugnal is unique given the samples,
+// and going analog -> digital -> analog is perfect.
+
+// That gives us point one: resampling from fs to ft is, semantically,
+// reconstructing the analog signal from the fs sampling, removing all
+// frequencies over ft/2, then sampling at ft.
+
+
+// Up-sampling theorem. Take a digital signal at frequency fs, and k
+// an integer > 1. Create a new digital signal at frequency fs*k by
+// alternatively taking one sample from the original signal and adding
+// k-1 zeroes. If one recreates the corresponding analog signal and
+// removes all frequencies over fs/2, then it will be identical to the
+// original analog signal, up to a constant multiplier on the
+// amplitude. For the curious the frequencies over fs/2 get copies of
+// the original spectrum with inversions, e.g. the frequency fs/2-a is
+// copied at fs/2+a, then it's not inverted at fs..fs*1.5, inverted
+// again between fs*1.5 and fs*2, etc.
+
+// A corollary is that if one starts for an analog signal with no
+// frequencies over fs/2, samples it at fs, then up-samples to fs*k by
+// adding zeroes, remove (filter) from the upsampled signal all
+// frequencies over fs/2 then reconstruct the analog signal you get a
+// result identical to the original signal. It's a perfect
+// upsampling, assuming the filtering is perfect.
+
+
+// Down-sampling theorem. Take a digital signal at frequency ft*k,
+// with k and integer > 1. Create a new digital signal at frequency
+// ft by alternatively taking one sample from the original signal and
+// dropping k-1 samples. If the original signal had no frequency over
+// ft/2, then the reconstructed analog signal is identical to the
+// original one, up to a constant multiplier on the amplitude. So it
+// is a perfect downsampling assuming the original signal has nothing
+// over ft/2. For the curious if there are frequencies over ft/2,
+// they end up added to the lower frequencies with inversions. The
+// frequency ft/2+a is added to ft/2-a, etc (signal to upsampling,
+// only the other way around).
+
+// The corollary there is that if one starts with a ft*k digital
+// signal, filters out everything over ft/2, then keeps only one
+// sample every k, then reconstruct the analog signal, you get the
+// original analog signal with frequencies over ft/2 removed, which is
+// reasonable given they are not representable at sampling frequency
+// ft anyway. As such it is called perfect because it's the best
+// possible result in any case.
+
+// Incidentally, the parasite audible frequencies added with the
+// wrapping when the original is insufficiently filtered before
+// dropping the samples are called aliasing, as in the high barely
+// audible frequencies that was there but not noticed gets aliased to
+// a very audible and annoying lower frequency.
+
+
+// As a result, the recipe to go from frequency fs to ft for a digital
+// signal is:
+
+// - find a frequency fm = ks*fs = kt*ft with ks and kt integers.
+// When fs and ft are integers (our case), the easy solution is
+// fm = fs * ft / gcd(fs, ft)
+
+// - up-sample the original signal x(t) into xm(t) with:
+// xm(ks*t) = x(t)
+// xm(other) = 0
+
+// - filter the resulting fm Hz signal to remove all frequencies above
+// fs/2. This is also called "lowpass at fs/2"
+
+// - lowpass at ft/2
+
+// - down-sample the fm signal into the resulting y(t) signal by:
+// y(t) = xm(kt*t)
+
+// And, assuming the filtering is perfect (it isn't, of course), the
+// result is a perfect resampling.
+
+// Now to optimize all that. The first point is that an ideal lowpass
+// at fs/2 followed by an ideal lowpass at ft/2 is strictly equivalent
+// to an ideal lowpass at min(fs/2, ft/2). So only one filter is
+// needed.
+
+// The second point depends on the type of filter used. In our case
+// the filter type known as FIR has a big advantage. A FIR filter
+// computes the output signal as a finite ponderated sum on the values
+// of the input signal only (also called a convolution). E.g.
+// y(t) = sum(k=0, n-1) a[k] * x[t-k]
+// where a[0..n-1] are constants called the coefficients of the filter.
+
+// Why this type of filter is pertinent shows up when building the
+// complete computation:
+
+// y(t) = filter(xm)[kt*t]
+// = sum(k=0, n-1) a[k] * xm[kt*t - k]
+// = sum(k=0, n-1) a[k] * | x[(kt*t-k)/ks] when kt*t-k is divisible by ks
+// | 0 otherwise
+// = sum(k=(kt*t) mod ks, n-1, step=ks) a[k] * x[(kt*t-k)/ks]
+
+// (noting p = (kt*t) mode ks, and a // b integer divide of a by b)
+// = sum(k=0, (n-1 - p))//ks) a[k*ks + p] x[(kt*t) // ks) - k]
+
+// Splitting the filter coefficients in ks phases ap[0..ks-1] where
+// ap[p][k] = a[p + ks*k], and noting t0 = (k*kt) // ks:
+
+// y(t) = sum(k=0, len(ap[p])-1) ap[p][k] * x[t0-k]
+
+// So we can take a big FIR filter and split it into ks interpolation
+// filters and just apply the correct one at each sample. We can make
+// things even easier by ensuring that the size of every interpolation
+// filter is the same.
+
+// The art of creating the big FIR filter so that it doesn't change
+// the signal too much is complicated enough that entire books have
+// been written on the topic. We use here a simple solution which is
+// to use a so-called zero-phase filter, which is a symmetrical filter
+// which looks into the future to filter out the frequencies without
+// changing the phases, and shift it in the past by half its length,
+// making it causal (e.g. not looking into the future anymore). It is
+// then called linear-phase, and has a latency of exactly half its
+// length. The filter itself is made very traditionally, by
+// multiplying a sinc by a Hann window.
+
+// The filter size is selected by maximizing the latency to 5ms and
+// capping the length at 400, which experimentally seems to ensure a
+// sharp rejection of more than 100dB in every case.
+
+// Finally, remember that up and downsampling steps multiply the
+// amplitude by a constant (upsampling divides by k, downsamply
+// multiply by k in fact). To compensate for that and numerical
+// errors the easiest way to to normalize each phase-filter
+// independently to ensure the sum of their coefficients is 1. It is
+// easy to see why it works: a constant input signal must be
+// transformed into a constant output signal at the exact same level.
+// Having the sum of coefficients being 1 ensures that.
+
+
+audio_resampler::audio_resampler(u32 fs, u32 ft)
+{
+ m_ft = ft;
+ m_fs = fs;
+
+ // Compute the multiplier for fs and ft to reach the common frequency
+ u32 gcd = compute_gcd(fs, ft);
+ m_ftm = fs / gcd;
+ m_fsm = ft / gcd;
+
+ // Compute the per-phase filter length to limit the latency to 5ms and capping it
+ m_order_per_lane = u32(fs * 0.005 * 2);
+ if(m_order_per_lane > 400)
+ m_order_per_lane = 400;
+
+ // Reduce the number of phases to be less than 200
+ m_phase_shift = 0;
+ while(((m_fsm - 1) >> m_phase_shift) >= 200)
+ m_phase_shift ++;
+
+ m_phases = ((m_fsm - 1) >> m_phase_shift) + 1;
+
+ // Compute the global filter length
+ u32 filter_length = m_order_per_lane * m_phases;
+ if((filter_length & 1) == 0)
+ filter_length --;
+ u32 hlen = filter_length / 2;
+
+ // Prepare the per-phase filters
+ m_coefficients.resize(m_phases);
+ for(u32 i = 0; i != m_phases; i++)
+ m_coefficients[i].resize(m_order_per_lane, 0.0);
+
+ // Select the filter cutoff. Keep it in audible range.
+ double cutoff = std::min(fs/2.0, ft/2.0);
+ if(cutoff > 20000)
+ cutoff = 20000;
+
+ // Compute the filter and send the coefficients to the appropriate phase
+ auto set_filter = [this](u32 i, float v) { m_coefficients[i % m_phases][i / m_phases] = v; };
+
+ double wc = 2 * M_PI * cutoff / (double(fs) * m_fsm / (1 << m_phase_shift));
+ double a = wc / M_PI;
+ for(u32 i = 1; i != hlen; i++) {
+ double win = cos(i*M_PI/hlen/2);
+ win = win*win;
+ double s = a * sin(i*wc)/(i*wc) * win;
+
+ set_filter(hlen-1+i, s);
+ set_filter(hlen-1-i, s);
+ }
+ set_filter(hlen-1, a);
+
+ // Normalize the per-phase filters
+ for(u32 i = 0; i != m_phases; i++) {
+ float s = 0;
+ for(u32 j = 0; j != m_order_per_lane; j++)
+ s += m_coefficients[i][j];
+ s = 1/s;
+ for(u32 j = 0; j != m_order_per_lane; j++)
+ m_coefficients[i][j] *= s;
+ }
+
+ // Compute the phase shift from one sample to the next
+ m_delta = m_ftm % m_fsm;
+ m_skip = m_ftm / m_fsm;
+}
+
+u32 audio_resampler::compute_gcd(u32 fs, u32 ft)
+{
+ u32 v1 = fs > ft ? fs : ft;
+ u32 v2 = fs > ft ? ft : fs;
+ while(v2) {
+ u32 v3 = v1 % v2;
+ v1 = v2;
+ v2 = v3;
+ }
+ return v1;
+}
+
+void audio_resampler::apply(const emu::detail::output_buffer_flat<sample_t> &src, std::vector<sample_t> &dest, u64 dest_sample, u32 srcc, float gain, u32 samples) const
+{
+ u32 seconds = dest_sample / m_ft;
+ u32 dsamp = dest_sample % m_ft;
+ u32 ssamp = (u64(dsamp) * m_fs) / m_ft;
+ u64 ssample = ssamp + u64(m_fs) * seconds;
+ u32 phase = (dsamp * m_ftm) % m_fsm;
+
+ const sample_t *s = src.ptrs(srcc, ssample - src.sync_sample());
+ sample_t *d = dest.data();
+ for(u32 sample = 0; sample != samples; sample++) {
+ sample_t acc = 0;
+ const sample_t *s1 = s;
+ const float *filter = m_coefficients[phase >> m_phase_shift].data();
+ for(u32 k = 0; k != m_order_per_lane; k++)
+ acc += *filter++ * *s1--;
+ *d++ += acc * gain;
+ phase += m_delta;
+ s += m_skip;
+ while(phase >= m_fsm) {
+ phase -= m_fsm;
+ s ++;
+ }
+ }
+}
+
+void audio_resampler::apply(const emu::detail::output_buffer_interleaved<s16> &src, std::vector<sample_t> &dest, u64 dest_sample, u32 srcc, float gain, u32 samples) const
+{
+ u32 seconds = dest_sample / m_ft;
+ u32 dsamp = dest_sample % m_ft;
+ u32 ssamp = (u64(dsamp) * m_fs) / m_ft;
+ u64 ssample = ssamp + u64(m_fs) * seconds;
+ u32 phase = (dsamp * m_ftm) % m_fsm;
+
+ gain /= 32768;
+
+ const s16 *s = src.ptrs(srcc, ssample - src.sync_sample());
+ sample_t *d = dest.data();
+ int step = src.channels();
+ for(u32 sample = 0; sample != samples; sample++) {
+ sample_t acc = 0;
+ const s16 *s1 = s;
+ const float *filter = m_coefficients[phase >> m_phase_shift].data();
+ for(u32 k = 0; k != m_order_per_lane; k++) {
+ acc += *filter++ * *s1;
+ s1 -= step;
+ }
+ *d++ += acc * gain;
+ phase += m_delta;
+ s += m_skip * step;
+ while(phase >= m_fsm) {
+ phase -= m_fsm;
+ s += step;
+ }
+ }
+}
+
+
+void audio_resampler::apply(const emu::detail::output_buffer_flat<sample_t> &src, std::vector<s16> &dest, u32 destc, int dchannels, u64 dest_sample, u32 srcc, float gain, u32 samples) const
+{
+ u32 seconds = dest_sample / m_ft;
+ u32 dsamp = dest_sample % m_ft;
+ u32 ssamp = (u64(dsamp) * m_fs) / m_ft;
+ u64 ssample = ssamp + u64(m_fs) * seconds;
+ u32 phase = (dsamp * m_ftm) % m_fsm;
+
+ gain *= 32768;
+
+ const sample_t *s = src.ptrs(srcc, ssample - src.sync_sample());
+ s16 *d = dest.data() + destc;
+ for(u32 sample = 0; sample != samples; sample++) {
+ sample_t acc = 0;
+ const sample_t *s1 = s;
+ const float *filter = m_coefficients[phase >> m_phase_shift].data();
+ for(u32 k = 0; k != m_order_per_lane; k++)
+ acc += *filter++ * *s1--;
+ *d += acc * gain;
+ d += dchannels;
+ phase += m_delta;
+ s += m_skip;
+ while(phase >= m_fsm) {
+ phase -= m_fsm;
+ s ++;
+ }
+ }
+}