I decided to give some of the ideas discussed in this thread a try. Below is a decimator implemented in Java. It uses a low-pass Butterworth filter to strip out frequencies above the Nyquist frequency (half the output sampling frequency) and it uses cubic splines to generate the output samples. Its constructor allows you to configure the input sampling frequency, the output sampling frequency and the filter order.

This image shows its frequency response curve for an order 16 filter with an output sampling frequency of 48000 Hz. The decimator uses 64-bit doubles, providing a reasonable range of stability. And, timing tests suggest that it does not consume a lot of CPU.

Samples to be decimated are added at the input sampling frequency via the addInputSample() method. The processOutputSample() method of the DecimatorListener, which is passed into the constructor, is called back at the output sampling rate with the decimated samples.

Per the discussion on this thread, I recommend generating audio samples every 3 CPU cycles. I.e. for NTSC, the input sampling rate would be 19687500.0 / (3.0 * 11.0) and for PAL, it would be 53203425.0 / (3.0 * 32.0). Unfortunately, if you try a higher rate, such as every other CPU cycle or every CPU cycle, at a filter order of 16, the Butterworth filter becomes unstable. For the output sampling frequency, I suggest either 44100.0 or 48000.0. The filter order needs to be a base-2 number and as mentioned, certain combinations of these values will result in an unusable filter.

Code: Select all

```
public class Decimator {
private static final double I2 = 1.0 / 2.0;
private static final double I3 = 1.0 / 3.0;
private static final double I6 = 1.0 / 6.0;
private final double[] inputSamples;
private final double[] outputSamples;
private final double[] A;
private final double[] B;
private final int lengthMask;
private final int filterOrder;
private final double cutoff;
private final double inputSamplingPeriod;
private final double ouputSamplingPeriod;
private final DecimatorListener listener;
private double time;
private int index;
// frequencies in Hz
public Decimator(
double inputSamplingFrequency,
double outputSamplingFrequency,
int filterOrder,
DecimatorListener listener) {
this.filterOrder = filterOrder;
this.listener = listener;
cutoff = outputSamplingFrequency / inputSamplingFrequency;
inputSamplingPeriod = 1.0 / inputSamplingFrequency;
ouputSamplingPeriod = 1.0 / outputSamplingFrequency;
inputSamples = new double[filterOrder << 1];
outputSamples = new double[inputSamples.length];
lengthMask = inputSamples.length - 1;
A = computeAs();
B = computeBs();
}
public void addInputSample(double value) {
inputSamples[index] = value;
double outputSample = B[0] * value;
for (int i = filterOrder; i >= 1; i--) {
outputSample += B[i] * inputSamples[(index - i) & lengthMask]
- A[i] * outputSamples[(index - i) & lengthMask];
}
outputSamples[index] = outputSample;
time += inputSamplingPeriod;
if (time >= ouputSamplingPeriod) {
listener.processOutputSample(computeSpline(
1.0 - (time - ouputSamplingPeriod) / inputSamplingPeriod,
outputSamples[(index - 3) & lengthMask],
outputSamples[(index - 2) & lengthMask],
outputSamples[(index - 1) & lengthMask],
outputSamples[index]));
time -= ouputSamplingPeriod;
}
index = (index + 1) & lengthMask;
}
private double computeSpline(double t, double y0, double y1, double y2,
double y3) {
double c0 = y1;
double c1 = y2 - I3 * y0 - I2 * y1 - I6 * y3;
double c2 = I2 * (y0 + y2) - y1;
double c3 = I6 * (y3 - y0) + I2 * (y1 - y2);
return ((c3 * t + c2) * t + c1) * t + c0;
}
private double[] computeAs() {
final double theta = Math.PI * cutoff;
final double sinTheta = Math.sin(theta);
final double cosTheta = Math.cos(theta);
double[] binomials = new double[filterOrder << 1];
for (int i = filterOrder - 1; i >= 0; i--) {
int i2 = i << 1;
double poleAngle = Math.PI * (i2 + 1.0) / (filterOrder << 1);
double sinPoleAngle = Math.sin(poleAngle);
double cosPoleAngle = Math.cos(poleAngle);
double v = 1.0 + sinTheta * sinPoleAngle;
binomials[i2] = -cosTheta / v;
binomials[i2 + 1] = -sinTheta * cosPoleAngle / v;
}
binomials = computeBinomials(binomials);
double[] as = new double[filterOrder + 1];
as[0] = 1.0;
as[1] = binomials[0];
as[2] = binomials[2];
for (int i = filterOrder; i >= 3; i--) {
as[i] = binomials[(i << 1) - 2];
}
return as;
}
private double[] computeBs() {
final double[] bs = new double[filterOrder + 1];
bs[0] = 1.0;
bs[1] = filterOrder;
for (int i = 2; i <= filterOrder >> 1; i++) {
bs[i] = (filterOrder - i + 1) * bs[i - 1] / i;
bs[filterOrder - i] = bs[i];
}
bs[filterOrder - 1] = filterOrder;
bs[filterOrder] = 1;
double scale = computeScale();
for (int i = 0; i < bs.length; ++i) {
bs[i] *= scale;
}
return bs;
}
private double computeScale() {
double omega = Math.PI * cutoff;
double sinOmega = Math.sin(omega);
double poleAngle = Math.PI / (filterOrder << 1);
double scale = 1.0;
for (int i = (filterOrder >> 1) - 1; i >= 0; i--) {
scale *= 1.0 + sinOmega * Math.sin(((i << 1) + 1.0) * poleAngle);
}
omega /= 2.0;
sinOmega = Math.sin(omega);
if ((filterOrder & 1) == 1) {
scale *= sinOmega + Math.cos(omega);
}
scale = Math.pow(sinOmega, filterOrder) / scale;
return scale;
}
private double[] computeBinomials(final double[] C) {
final int N = C.length >> 1;
final double[] D = new double[N << 1];
for (int i = 0; i < N; i++) {
int i2 = i << 1;
for (int j = i; j > 0; j--) {
int j2 = j << 1;
D[j2] += C[i2] * D[j2 - 2] - C[i2 + 1] * D[j2 - 1];
D[j2 + 1] += C[i2] * D[j2 - 1] + C[i2 + 1] * D[j2 - 2];
}
D[0] += C[i2];
D[1] += C[i2 + 1];
}
return D;
}
}
```

Code: Select all

```
public interface DecimatorListener {
void processOutputSample(double sample);
}
```