libxmp/docs/modguide/interpolate.txt

Dear Amiga mod playback software author,

you are receiving this note because your software is implementing playback
capabilities for an important format, the Amiga protracker 1.x-3.x mod format.

As you probably are aware, most of these modules were composed on the Amiga
computer for playback on Amiga system. The Amiga hardware playback system is
unlike the PC environment, and needs some special care if authentic-seeming
Amiga mod playback is desired.

Typical playback software for the various sample tracker formats is architected
roughly like this:

    [ Player routine ] --> [ Mixer and resampler ] --> Output

The player routine component produces information about the samples, their
sampling frequencies and volume levels (and filtering, panning, etc. depending
on format) to the mixer and resampler component, which uses the player
routine's instructions to generate the final output by resampling the sample
stream from memory and mixing it to physical output channels.  I'm going to
focus exclusively to the Mixer and resampler component here.

A typical implementation of the resampler uses linear interpolation, or some
higher-order polynomial interpolation, splines, or even sinc-based
"theoretically correct" resampling with some 8-tap FIR or more. Unfortunately,
all of them miss the mark as far as Amiga is concerned. For typical mod player,
the most accurate emulation of Amiga output occurs when *no* interpolation is
used at all. However, there seems to be nothing special about the mixing --
simply adding the resampled channels together suffices.

Why the Amiga is different?
---------------------------

To understand what goes on with the Amiga sound chip, Paula, is to understand that
Paula does no interpolation of any kind. Paula's output is strictly a pulse wave,
produced on 3546895 Hz frequency, which is the Paula clock rate for PAL systems.

The Amiga period values, typically valued somewhere between 120 and 800, count Paula
ticks. For instance, a period value of 400 means that Paula waits 400 ticks holding
the output constant, and then changes to the next sample.

To simulate this perfectly, we would produce an output sample stream at that
3.55 MHz rate, one sample per tick, and then simply hold the output constant
while counting down to zero from the period value, then switch to next sample,
reset period clock back to original value, and so on.  This stream would be
very accurate representation of the Paula's output.

Aliasing in Paula output
------------------------

It is also worth it to understand that playing back a sample from Paula
replaces the original waveform with its smoothly undulating shapes with nothing
but hard-edged pulses. This will have an impact to the sound, of course, and is
called aliasing distortion. It's called aliasing, because copies of the
original sound's frequency spectrum repeat throughout the entire audible part
of audio spectrum and beyond, theoretically up to infinity.

This may surprise those who believe that because Amiga's maximum sampling rate
is some 28867 Hz, the maximum frequency produced by Amiga is threrefore limited
to 14.5 kHz or so. However, this is merely the highest frequency that may be
produced without aliasing.

In addition, some misguided people have tried to run their mod players at 28867
Hz mixing frequency based on some kind of confused idea of what the mixing rate
means for Amiga. This won't work either.

The analog filters -- more complications
----------------------------------------

After leaving the chip, the sound enters a 6 dB/oct resistor+capacitor (RC)
low-pass filter tuned at approximately 5 kHz. I will call this component "fixed
filter" in this document. The purpose of this filter is probably to remove some
of that aliasing which is inherent in the pulse synthesis. Additionally, it is
possible to engage a low-pass 12 dB/oct Butterworth filter tuned at
approximately 3.2 kHz by turning the Amiga power LED brighter with a special
protracker command.  I will call this component "LED filter".

These filters should probably be understood as reconstruction filters, because
their purpose is to attempt to "reconstruct" the original smoothly undulating
waveform after a pulse-based version of it has been generated by some D/A
converter.  However, proper reconstruction filters should probably feature
steeper cut-offs than just 6 or 12 dB, and they should be dynamic; their tuning
should follow the sampling rate, as it is the sampling rate that decides where
the first aliased frequency will occur. (As an example, some Atari ST models
have a DAC chip which is fed through DMA at any of the four possible sampling
rates, and the analog reconstruction filter's cutoff also is chosen based on
the sampling rate.)

Regardless, this is the hand we have been dealt with, and in order to
accurately simulate Amiga 500 audio playback, we thus need to emulate these two
filters in addition of producing the right kind of output stream from Paula
emulation.

No interpolation -- use resampling!
-----------------------------------

To properly resample audio, we need to achieve two things, often done simultaneously:

* remove all components in the audio stream above some selected frequency, usually
  about 20 kHz.
* decimate the signal by discarding samples until we obtain output at the correct rate.
  For 44.1 kHz output we need to discard approximately 79 of every 80 samples.

Decimation can be safely done only after the lowpass step has been completed.
If decimation is done without lowpass filtering, some extra aliasing shall be
introduced to the remaining audio. What happens is that all frequencies in the
spectrum now alias in the audible range, making the ultrasonic sounds produced by
the pulse-wave synthesis to appear as additional noise in the mod playback.

Fortunately, the previous paragraph contains the answer to the problem. It is
possible to make simulated Paula's output consist of "bandlimited steps", or
BLEPs as they are called. The output is done by simulating Paula at 3.55 MHz
clock rate and replacing the hard edges of Paula's pulses by BLEPs, which are
carefully constructed in such a fashion to not contain those difficult frequencies
above 21 kHz or so.

Because the input stream does not contain frequencies near or above nyquist
frequency of the final sampling rate, all we have to do is to build an abstract
version of Paula's output signal that we can evalute to real sample whenever it
is needed. The abstract version consists two parts: the output level of
"genuine" Paula with hard-edged pulses, and the states of various BLEPs which
which are then used to correct the "hard edges" of Paula's waveform to the
smooth curves of the BLEPs instead.

The mixing of one blep thus has the following procedure, roughly:

* adjust the output level of Paula to the new value.
* start a new blep, and scale it in such a fashion that its amplitude is the difference
  between new output level and the old level.

If we ask the synthesis function to produce a sample at this moment, it would still
report the old level. But as we clock Paula forward, and ask again, the blep gradually
changes the output level from the old value to the new value. At some point we are done
with the BLEP, and then the output is simply constant at the new output level.

This leaves the filters, but they are easy. Because the output is done by
mixing bleps, we can implement the analog filters in the bleps themselves, by
filtering them with digital models of the analog filters to produce versions of
BLEPs that have treble components reduced to match Amiga 500's operation with
LED filter on and off. This may seem somewhat miraculous, but this is how it
works.

The implementation
------------------

We need to compute the blep tables. Firstly we need a function that contains
all the frequencies in the world. One such function is called the delta
function. We define it as follows:

    delta(x) = 0 if x != 0
    delta(x) = 1 if x = 0.

If we low-pass filter the delta, it morphs into a sinc function. sinc function
is a special one, it has all the frequencies from 0 to some given cutoff value.
The cutoff is determined by the rate the sinc is oscillating in.

    f(x, rate) = sin(x * pi / rate) / (x * pi / rate) if x != 0
    f(0, rate) = 1 if x = 0

With rate = 1.0 the sinc does no filtering at all, and becomes a delta function
in effect. If we plot the sinc function by asking the value of sinc at each
integer value, it is zero everywhere except 1 at 0. To do no filtering at all
is the same as having the cutoff at exactly half of the sampling rate. If the
sinc is to remove everything above, say, 21 kHz, rate = 3.55e6 / 21e3 / 2.

An ideal sinc filter extends from negative infinity to the positive infinity.
Because a practical implementation is needed, the sinc function must be
truncated at some point.  The truncation affects the width of the stopband of
the sinc function -- short functions begin filtering earlier and are not as
steeply filtering as long filters.  I've found that 2048 point sinc functions
are sufficient to perform filtering for Paula, so we'd sample the sinc function
from f(-1024, rate) to f(1023, rate) and thus yield 2048 values.

Because 2048 point function is quite short, there is an issue at the edges,
where the sinc is still oscillating. Some windowing function is needed to
smooth the sinc properly at edges to zero. The Hann window, or the raised
cosine window is the simplest. It is defined as hann(x, length) = 0.5 - 0.5 *
cos(x / length * 2 * pi). For this, the length would be 2048, and the values of
x would go from 0 to 2047.

The result should look like a symmetric "water-drop" with few small side lobes.
If you have some plotting software, you can see the sinc function at this
point. gnuplot suffices, you can type something like:

set xrange [-1024:1023]
plot sin(x * pi / 84) / (x * pi / 84) * (0.5 - 0.5 * cos((x + 1024) / 2048 * 2 * pi))

The sinc table is final mixing frequency independent
-----------------------------------------------------

The only things that go into calculating the sinc table are the Paula clock
frequency and the frequency of the cutoff. As long as calls to obtaining output
occur at rates greater or equal to 44.1 kHz, the table works. With rates below
about 20 kHz, the blep table itself will produce aliasing.

The function is sampled at Paula clock rate, which is constant; it therefore
also has a fixed duration in seconds. For 2048 points, it's 2048 paula ticks or
roughly 0.6 milliseconds. From this it follows that at 44.1 kHz sampling
frequency, about 25 samples of the table make it into the output.

The table must however be longer than merely 25 samples because depending on
the phase where the blep starts, ie. paula pulse occurs, it's always different
values of the table that need to be chosen, so all of them become used at some
point. (However, in truth the table does not need to be that precise, it could
well be sampled at a lower frequency without compromising output quality
noticeably. Perhaps only about 16 different phase values are needed, but now it
has 80 (as many phases as there are clocks between samples).)

Filtering the sinc
------------------

The minimum-phase reconstruction is implemented in cepstral domain in the
included program.  My knowledge of the operations involved is too limited to
explain why or how it works, but after we do this, the filter shape changes in
such a way that as much as possible of the filter occurs "early on" and the
filter "settles" for longer period after an initial large transition. So,
before minimum phase reconstruction, the filter is a symmetric waterdrop-like
wave, with pre- and postringing at the peak, but after the reconstruction the
peak moves to the start of filter, and the filter displays decaying oscillation
for most of its duration.  The upshot of the reconstruction is that it allows
our IIR filters the longest period to settle after the initial excitation by
the large part of the peak, thus allowing us to truncate them earlier.

There has been some talk that minimum-phase reconstructions are in general
better audio filters than linear-phase filters (which are symmetric by the
midpoint), because they do not display the preringing that is perceptible to
ear if it occurs for too long. The increased postringing is not an issue, as it
often occurs in nature, and thus the ear "expects" it.  However, the issue is
made complicated by the fact that phase relationships are by definition
modified by minimum-phase reconstructions. At any case, this choice makes no
audible difference in this case, because in both cases the ringing is tuned to
occur at frequencies above 20 kHz, and the phase delays are almost linear
across the entire audible band.

The fixed filter
----------------

We should apply the IIR filters to the sinc function before proceeding to build the
band-limited step.  The 6 dB/oct filter is an IIR filter that is computed
according to the following formula:

    output(n) = b0 * input(n) + (1 - b0) * output(n - 1)

that is, the output is a linear combination of input and the value of output in the
last iteration.

The b0 determines the cutoff value and is computed as follows:

    omega = 2 * pi * 5000 / sampling_rate
    b0 = 1 / (1 + 1 / omega)

As a reminder, the cutoff frequency would be 5.0 kHz and sampling rate is 3.55
MHz because we will operate the blep that is "sampled" at 3.55 MHz. The rest of
the problem here is just applying the previous formula to filtering the BLEP.

Amiga 1200 current leakage simulation
-------------------------------------

It seems that Amiga 1200 output response attenuated by some 0.2 - 0.3 dB near
20 kHz.  This would seem to be due to the LED filter even though it isn't
enabled, because this effect does not manifest when turning the LED filter on.
It is difficult to be sure, though, because of limited resolution and the
extreme attenuation, noise floor gets very close to signal. Regardless, the
effect can be modelled with a 32 kHz lowpass RC filter.

A word of caution
-----------------

To correctly apply the filter, we first need to extend the sinc we computed by
imagining that it is preceded by, say, 10 000 constant values of the first
value of the filter, and these are run through the filter to ensure that the
filter is properly settled for the base level of the blep. This is especially
important if the first value of the sinc is not zero, or if the filter is
reused between computations for any reason.

Additionally, to validate that the length of the blep is sufficient, it would
be helpful to examine how much the IIR filter oscillates at the end of the
table by computing the value of

    filter(last_output_value) - last_output_value

If the filter has not yet damped, the difference will not be zero, and depending of
the oscillation magnitude, artifacts could be introduced.

The LED filter
--------------

If no LED filtering emulation is desired, the next BLEP can be left uncomputed.
To simulate the LED filter, we clone the previous sinc function for applying
the LED filter on.  The model for LED filter requires application of a biquad
filter that simulates the 2nd order Butterwoth lowpass filter.

The formula for biquad is as follows:

    output(n) = b0 * input(n) + b1 * input(n-1) + b2 * input(n-2)
        - a1 * output(n - 1) - a2 * output(n - 2)

That is, two samples of past input history and two samples of past output history and
the current sample are required. The coefficients are derived using bilinear transform
in the Python source; feel free to look at the gory details for computing b0, b1, b2,
a1 and a2 from there.

So, anyway, we should go and filter our copy of the previous sinc function
with the biquad formula, and then proceed to figuring out how to convert the
filtered sincs to bandlimited steps.

The step-function
-----------------

Step is our representation of Paula's pulses. We could define it as follows:

    f(x) = 0 if x < 0
    f(x) = 1 if x >= 0

The band-limited step
---------------------

Convolution in DSP is the application of a filter in the time domain. It is
basically like computing a dot product between two vectors. The other "vector"
is the tabularised version of the convolution function, and the other is the
sample data around the point we wish to compute convolution at. We would be
talking about 2048-dimensional vector here, if vectors were otherwise an
appropriate metaphor.

In order to bandlimit a step, we can calculate a representation of the step
sample stream.  We start with the entire sample buffer at zero and calculate
the convolution. It is all zero as expected.

Then we introduce the edge, or the first sample of the step into the sample
buffer. The convolution's value is now equal to the first value in the
convolution function. As more samples are introduced, the value of convolution
is the sum of the values in the convolution function over the area of the step
function. Thus, the bandlimited step is in effect just a numeric integral of
the convolution function.

We should take care to scale the step accordingly. If you are dealing with
16-bit sample data, you would probably want to scale the step in such a fashion
that it starts at, say, 0 and ends at 65536. Or, as I often prefer, it starts
at -65536 and ends at 0. This is just to make the blep evaluation code simple.
The downside with starting with negative number is that the sign takes 1 byte
to encode in the resultant C source. For the example program, the blep
resolution is 17 bits and the blep runs from a positive value to zero.

Sample code in C (but for mono output)
--------------------------------------

/* 131072 to 0, 2048 entries */
#define BLEP_SCALE 17
#define BLEP_SIZE 2048

/* the step function */
int32_t sine_integral[BLEP_SIZE] = { precomputed table };

/* the instantenous value of Paula output */
static int16_t global_output_level;

/* count of simultaneous bleps to keep track of */
static unsigned int active_bleps;

/* the structure that holds data of bleps */
typedef struct {
    int16_t level;
    int16_t age;
} blep_state_t;

/* place to keep our bleps in. MAX_BLEPS should be
   defined as a BLEP_SIZE / MINIMUM_EVENT_INTERVAL.
   For Paula, minimum event interval could be even 1, but it makes sense to
   limit it to some higher value such as 16. */
static blep_state_t    blepstate[MAX_BLEPS];

/* return output simulated as series of bleps */
int16_t
output_sample(int filter_table) {
    int i;
    int32_t output;

    output = global_output_level << BLEP_SCALE;
    for (i = 0; i < active_bleps; i += 1)
        output -= winsinc_integral[filter_table][blepstate[i].age] * blepstate[i].level;
    output >>= BLEP_SCALE;
    
    if (output < -32768)
        output = -32768;
    else if (output > 32767)
        output = 32767;
    
    return output;
}

In the code where new Paula samples get added to the system, you would so
something like this:

function input_sample(int16_t sample) {
    if (sample != global_output_level) {
        /* Start a new blep: level is the difference, age (or phase) is 0 clocks. */
        if (active_bleps > MAX_BLEPS - 1) {
            fprintf(stderr, "warning: active blep list truncated!\n");
            active_bleps = MAX_BLEPS - 1;
        }
        /* Make room for new blep */
        memmove(&blepstate[1], &blepstate[0], sizeof(blepstate[0]) * active_bleps);
        /* Update state to account for the new blep */
        active_bleps += 1;
        blepstate[0].age = 0;
        blepstate[0].level = sample - global_output_level;
        global_output_level = sample;
    }
}

And while no events occur, you would call the clock() function to advance the bleps.

void clock(unsigned int cycles) {
    for (i = 0; i < active_bleps; i += 1) {
        blepstate[i].age += cycles;
        if (blepstate[i].age >= BLEP_SIZE) {
            active_bleps = i;
            break;
        }
    }
}

This code is not very optimized in that it constantly keeps blepstate as a
sorted list. It is possible to implement it more efficiently, for instance by
keeping occurence times instead of ages, but that complicates the output sample
function.  The right optimization of the code might depend on whether periods
are small or large compared to sampling intervals. For Paula, periods are
typically much larger than sampling intervals, so code like this might be
better tradeoff. For modern computers, it takes quite little time, so I leave
the optimizing for people who are that way inclined. Typical length of the blep
list can be short, only 5-6 entries if the sample rates are low or occur
conveniently otherwise. At the maximum, Paula is constantly producing pulse
waves but it's probably sensible to limit maximum period to 124 unless Paula
DMA fetch starvation is simulated.

A typical sequence of events would look like this. Assuming that we sample output
at exactly 80 paula clock ticks apart, and that playback period is, say 140, we'd do
calls like this:

    input_sample(samplebuf[0]);
    clock(80);
    *output++ = output_sample();
    clock(60);
    input_sample(samplebuf[1]);
    clock(20);
    *output++ = output_sample();
    clock(80);
    *output++ = output_sample();
    clock(40);
    input_Sample(samplebuf[2]);
    ...

Results and analysis
--------------------

I have hi-fi recordings of Amiga 500 and 1200 playing Sainahi Circles with both
LED filter on and off. This analysis was largerly based on the sampling work by
pyksy on #amigaexotic at IRCnet. Big thanks to pyksy for his work! In addition,
Melkor (kohina board) sent me a hi-fi sampling of his Amiga 4000 playing Sainahi
Circles.

The results are summarised as number of png pictures that display the simulated
and real responses side-by-side and their residual (difference) spectrums. The
units are in Hz and dB.  There is a constant offset that follows from the fact
that the hi-fi samplings used 32-bit floating point samples, whose value range
is different from the integer-coded PCM data, yielding the largest part of the
difference. The rest is explained by the different recording levels used
between setups.

The files called a{1200,500,4000}_{on,off}.png shows the simulated response of
UADE and real Amiga superimposed. comb{1200,500,4000}_{on,off}.png show the
difference between two in dB units. These are the same pictures as before, but
with the spectra substracted, allowing for more precise evaluation of the
differences between simulated and real output.

We can see that UADE and Amiga produce same features closely, thus essentially
validating the model used for output synthesis.

We can observe a good cutoff starting around 20 kHz and quickly extinguishing
response by approximately 23 kHz. This should yield satisfactory results for
all sampling frequencies above and including 44.1 kHz. However, I personally
recommend 48 kHz as default output frequency because it works without
resampling on majority of cards, while 44.1 kHz output often may need a
resampling step on either software or hardware.

For LED off cases, the residual spectrum is largerly just noise. There is,
however a small 0.3 dB boost around 0 to 3 kHz. This can be seen on Amiga 500,
but not on Amiga 1200. In both Amiga 500 and Amiga 1200, the LED filter also
displayed this boost, and it was compensated for by adjusting the Butterworth
filter with a small resonance term.

Despite adjustments, the LED filter case is not perfect. The residual spectrum
is roughly flat to 3 kHz, thanks to manually fitting the filter, but after that
there is some extra attenuation which is then compensated by something else on
the analog filter.  For Amiga 500 case, the simulation hits noise floor and can
not yield anything meaningful above 15 kHz, but those measurements also show
the same general pattern.

Overall, I think the model yields accuracy to about 0.5 dB for all cases, which
should be sufficient for all practical purposes. Especially the Amiga 1200 case
seems well modelled.

Author
------

Antti S. Lankila <alankila@bel.fi>

feel free to send feedback, comments and spam.