1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
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 ++;
}
}
}
|