Some advice of DSP here...

// Precise method which guarantees v = v1 when t = 1.
float lerp(float v0, float v1, float t) {
  return (1-t)*v0 + t*v1;
}

void linearinterpolation(float * src, int len, float * dest, int f)
{
   int i;
   float sample;
   for (i = 0; i < len - 1; i++)
   {
      sample = lerp(src[i], src[i+1], f);
      dest[i] = sample;
   }
}

for (i = 0; i < len_out; ++i)
{
    // convert output position to input position
    float pos_in = float(i) * float(len_in) / float(len_out);
    // calculate blend value for linear interpolation
    float alpha = pos_in - floorf(pos_in);
    // find the two samples to interpolate
    int p0 = int(pos_in);
    int p1 = p0 + 1;
    dest[i] = lerp(src[p0], src[p1], alpha);
}

The emulator merely sets the input clock rate and output sample rate, adds waveforms by specifying the clock times where their amplitude changes, then reads the resulting output samples.

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;
  } 
}

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

Fs = 1786840;
cutoff = 20000;
%[z, p, k] = ellip(16, .1, 90, cutoff/(Fs/2));
[z, p, k] = butter(16, cutoff/(Fs/2));
[sos, g] = zp2sos(z, p, k);

Bs = sos(:,1:3);
As = sos(:,4:6);
[nsec,temp]=size(sos);
nsamps = 16384;
x = g*[1,zeros(1,nsamps-1)];
for i=1:nsec
 x = filter(Bs(i,:),As(i,:),x);
end
X = fft(x);
f=[0:nsamps-1]*Fs/nsamps;
figure(1);
grid('on');
axis([0 Fs/2 -100 5]);
legend('off');
plot(f(1:nsamps/2),20*log10(X(1:nsamps/2)));

%[z, p, k] = ellip(16, .1, 90, cutoff/(Fs/2));
[z, p, k] = butter(16, cutoff/(Fs/2));

fc = 24000; % Cut-off frequency (Hz)
fs = 19687500/11; % Sampling rate (Hz)
order = 8; % Filter order
[B,A] = butter(order,2*fc/fs); % [0:pi] maps to [0:1] here
[sos,g] = tf2sos(B,A)
Bs = sos(:,1:3); % Section numerator polynomials
As = sos(:,4:6); % Section denominator polynomials
[nsec,temp] = size(sos);
nsamps = 4096;

x = g*[1,zeros(1,nsamps-1)];
for i=1:nsec
  x = filter(Bs(i,:),As(i,:),x);
end
figure(2);
X=fft(x);
f = [0:nsamps-1]*fs/nsamps; grid('on');
axis([0 50000 -100 5], "manual"); legend('off');
plot(f(1:nsamps/2),20*log10(X(1:nsamps/2)));

Multiply an array of input samples with the array of weights (FIR filter) and sum them up to produce an output sample.

clear
fc = 20000; % Cut-off frequency (Hz)
fs = 19687500.0 / 11.0; % NTSC Sampling rate (Hz)
%fs = 53203425.0 / 32.0; % PAL Sampling rate (Hz)
order = 13; % Filter order

[z, p, k] = ellip(order, 0.1, 100, 2*fc/fs);

[sos, g] = zp2sos(z, p, k);

Bs = sos(:,1:3)
num2hex(Bs)
As = sos(:,4:6)
num2hex(As)
[nsec,temp]=size(sos);
nsamps = 100000;
x = g*[1,zeros(1,nsamps-1)];
for i=1:nsec
 x = filter(Bs(i,:),As(i,:),x);
end
X = fft(x);
f=[0:nsamps-1]*fs/nsamps;
figure(1);
grid('on');
axis([0 fc*2 -100 5]);
legend('off');
pscale = 48
plot(f(1:nsamps/pscale),20*log10(X(1:nsamps/pscale)));

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 static final double NTSC_SAMPLING_FREQUENCY = 19687500.0 / 11.0;
  private static final double PAL_SAMPLING_FREQUENCY = 53203425.0 / 32.0;
 
  private static final long[] LONG_NTSC_A1 = {
    0xbfff67701c84a3deL,
    0xbfff8b62242e56e4L,
    0xbfffaf25694786a4L,
    0xbfffc959eca889e2L,
    0xbfffda7652351a06L,
    0xbfffe6495d699ef1L,
    0xbfef576a39439bdaL,
  };
 
  private static final long[] LONG_NTSC_A2 = {
    0x3feed6d89d961d3bL,
    0x3fef28cd096bf76aL,
    0x3fef7a3a53cd4f67L,
    0x3fefb5a2213e3abbL,
    0x3fefdbe39ed2f06eL,
    0x3feff554ab2f0006L,   
  };
 
  private static final long[] LONG_NTSC_B1 = {
    0xbfff2a2df4fd2ad2L,
    0xbfffbe3738d86ac8L,
    0xbfffd94c2258d55cL,
    0xbfffe22ba0cf6805L,
    0xbfffe5b14c11db72L,
    0xbfffe6ffa75a14ebL,
  };
 
  private static final long LONG_NTSC_G = 0x3ecc208a2d079678L;
 
  private static final double[] NTSC_A1;
  private static final double[] NTSC_A2;
  private static final double[] NTSC_B1;
  private static final double NTSC_G;
 
  static {
    NTSC_A1 = toDoubleArray(LONG_NTSC_A1);
    NTSC_A2 = toDoubleArray(LONG_NTSC_A2);
    NTSC_B1 = toDoubleArray(LONG_NTSC_B1);
    NTSC_G = Double.longBitsToDouble(LONG_NTSC_G);
  }
 
  private static final long[] LONG_PAL_A1 = {
    0xbfff5baa1b81a309L,
    0xbfff81d9b2eb15feL,
    0xbfffa7df44761c9dL,
    0xbfffc3c2600f5556L,
    0xbfffd5ff6fbab417L,
    0xbfffe2a551b417d0L,
    0xbfef4aa6ceb6c246L,
  };
 
  private static final long[] LONG_PAL_A2 = {
    0x3feec08d9ac7eb51L,
    0x3fef18943f900329L,
    0x3fef7018a2a55ed4L,
    0x3fefaffb4199e44aL,
    0x3fefd9237788b71bL,
    0x3feff4844bf56456L, 
  };
 
  private static final long[] LONG_PAL_B1 = {
    0xbfff08b370e2ff76L,
    0xbfffb3ce724731b7L,
    0xbfffd3295d0112b5L,
    0xbfffdd7036a20286L,
    0xbfffe184a77ef413L,
    0xbfffe307f8bb95fdL,
  };
 
  private static final long LONG_PAL_G = 0x3ece3d98b822795aL;
 
  private static final double[] PAL_A1;
  private static final double[] PAL_A2;
  private static final double[] PAL_B1;
  private static final double PAL_G;
 
  static {
    PAL_A1 = toDoubleArray(LONG_PAL_A1);
    PAL_A2 = toDoubleArray(LONG_PAL_A2);
    PAL_B1 = toDoubleArray(LONG_PAL_B1);
    PAL_G = Double.longBitsToDouble(LONG_PAL_G);
  } 
 
  private final double g;
  private final double[] a1;
  private final double[] a2;
  private final double[] b1;
  private final double[] d1 = new double[7];
  private final double[] d2 = new double[6];
  private final double[] ys = new double[4];
  private final double inputSamplingFrequency;
  private final double inputSamplingPeriod;
  private final double outputSamplingPeriod;
  private final DecimatorListener listener;
 
  private double time;
  private int index;
 
  public Decimator(boolean ntsc, double outputSamplingFrequency,
      DecimatorListener listener) {
   
    if (ntsc) {
      a1 = NTSC_A1;
      a2 = NTSC_A2;
      b1 = NTSC_B1;
      g = NTSC_G;
      inputSamplingFrequency = NTSC_SAMPLING_FREQUENCY;     
    } else {
      a1 = PAL_A1;
      a2 = PAL_A2;
      b1 = PAL_B1;
      g = PAL_G;
      inputSamplingFrequency = PAL_SAMPLING_FREQUENCY;
    }
   
    this.listener = listener;
   
    inputSamplingPeriod = 1.0 / inputSamplingFrequency;
    outputSamplingPeriod = 1.0 / outputSamplingFrequency;
  }
 
  private static double[] toDoubleArray(long[] values) {
    double[] ds = new double[values.length];
    for(int i = values.length - 1; i >= 0; i--) {
      ds[i] = Double.longBitsToDouble(values[i]);
    }
    return ds;
  }
 
  public void addInputSample(double x) {
   
    double y;
    for(int i = 0; i < 6; i++) {
      y = x + d1[i];   
      d1[i] = b1[i] * x - a1[i] * y + d2[i];
      d2[i] = x - a2[i] * y;
      x = y;
    }
    y = x + d1[6];
    d1[6] = x - a1[6] * y;       
    ys[index] = g * y;
   
    time += inputSamplingPeriod;
    if (time >= outputSamplingPeriod) {
      listener.processOutputSample(computeSpline(
          1.0 - (time - outputSamplingPeriod) * inputSamplingFrequency,
          ys[(index - 3) & 3],
          ys[(index - 2) & 3],
          ys[(index - 1) & 3],
          ys[index]));
      time -= outputSamplingPeriod;
    }
   
    index = (index + 1) & 3;
  }
   
  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;
  }   
}

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

public class AveragingDecimator {

  private final double inputSamplingPeriod;
  private final double outputSamplingPeriod;
  private final DecimatorListener listener;
 
  private double time;
  private double sumOfSamples;
  private int numberOfSamples;
 
  public AveragingDecimator(
      double inputSamplingFrequency,
      double outputSamplingFrequency,
      DecimatorListener listener) {

    this.inputSamplingPeriod = 1.0 / inputSamplingFrequency;
    this.outputSamplingPeriod = 1.0 / outputSamplingFrequency;
    this.listener = listener;
  }
 
  public void addInputSample(double x) {
   
    sumOfSamples += x;
    numberOfSamples++;
   
    time += inputSamplingPeriod;
    if (time >= outputSamplingPeriod) {
      listener.processOutputSample(sumOfSamples / numberOfSamples);
      sumOfSamples = 0;
      numberOfSamples = 0;
      time -= outputSamplingPeriod;
    }
  }
}

if (Sound.cc_counter >= samples_div)
{
   float sample;
   sample = avgsamples(Sound.X, Sound.cc_counter);
   sample *= amp_scale;
   sample *= 0x07FFF;
   sample = clip_sample(sample);
   Sound.buffer[Sound.buffer_index] = sample;
   Sound.buffer_index++;
   Sound.amp = 0;
   Sound.cc_counter -= samples_div;
}

float avgsamples(float * out_buffer, int len)
{
   float sample = 0;
   for (int i = 0; i < len; i++)
      sample += out_buffer[i];
   return sample / (float) len;
}

#define ChangeScale()(amp_scale = pow(1.122018454, g_db));

static inline float clip_sample(float sample_value) {
  if (sample_value < -32000) return -32000;
  if (sample_value > 32000) return 32000;
  return sample_value;
}

y[i] = a * x[i] + (1-a) * y[i-1];

y[n] = b0 * x[n] + d1;       
d1 = b1 * x[n] - a1 * y[n] + d2;       
d2 = b2 * x[n] - a2 * y[n];

// Return RC low-pass filter output samples, given input samples,
 // time interval dt, and time constant RC
 function lowpass(real[0..n] x, real dt, real RC)
   var real[0..n] y
   var real α := dt / (RC + dt)
   y[0] := x[0]
   for i from 1 to n
       y[i] := α * x[i] + (1-α) * y[i-1]
   return y

   var real α := dt / (RC + dt)

for i from 1 to n
    y[i] := α * x[i] + (1-α) * y[i-1]