SHARE

TWEET

# Untitled

a guest
Jan 12th, 2019
61
Never

**Not a member of Pastebin yet?**

**, it unlocks many cool features!**

__Sign Up__- 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.

RAW Paste Data

We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.