|
| 1 | +#!/usr/bin/python |
| 2 | + |
| 3 | +import math |
| 4 | +import sys |
| 5 | +from scipy import array, exp |
| 6 | +from scipy.signal import firwin |
| 7 | +from scipy.fftpack import fft, ifft |
| 8 | + |
| 9 | +def fir_minphase(table, pad_size=8): |
| 10 | + table = list(table) |
| 11 | + # table should be a real-valued table of FIR coefficients |
| 12 | + convolution_size = len(table) |
| 13 | + table += [0] * (convolution_size * (pad_size - 1)) |
| 14 | + |
| 15 | + # compute the real cepstrum |
| 16 | + # fft -> abs + ln -> ifft -> real |
| 17 | + cepstrum = ifft(map(lambda x: math.log(x), abs(fft(table)))) |
| 18 | + # because the positive and negative freqs were equal, imaginary content is neglible |
| 19 | + # cepstrum = map(lambda x: x.real, cepstrum) |
| 20 | + |
| 21 | + # window the cepstrum in such a way that anticausal components become rejected |
| 22 | + cepstrum[1 :len(cepstrum)/2] *= 2; |
| 23 | + cepstrum[len(cepstrum)/2+1:len(cepstrum) ] *= 0; |
| 24 | + |
| 25 | + # now cancel the previous steps: |
| 26 | + # fft -> exp -> ifft -> real |
| 27 | + cepstrum = ifft(map(exp, fft(cepstrum))) |
| 28 | + return map(lambda x: x.real, cepstrum[0:convolution_size]) |
| 29 | + |
| 30 | +class BiquadFilter(object): |
| 31 | + __slots__ = ['b0', 'b1', 'b2', 'a1', 'a2', 'x1', 'x2', 'y1', 'y2'] |
| 32 | + |
| 33 | + def __init__(self, b0, b1, b2, a1, a2): |
| 34 | + self.b0 = b0 |
| 35 | + self.b1 = b1 |
| 36 | + self.b2 = b2 |
| 37 | + self.a1 = a1 |
| 38 | + self.a2 = a2 |
| 39 | + self.reset() |
| 40 | + |
| 41 | + def reset(self): |
| 42 | + self.x1 = 0.0 |
| 43 | + self.x2 = 0.0 |
| 44 | + self.y1 = 0.0 |
| 45 | + self.y2 = 0.0 |
| 46 | + |
| 47 | + def filter(self, x0): |
| 48 | + y0 = self.b0*x0 + self.b1*self.x1 + self.b2*self.x2 - self.a1*self.y1 - self.a2*self.y2 |
| 49 | + self.x2 = self.x1 |
| 50 | + self.x1 = x0 |
| 51 | + |
| 52 | + self.y2 = self.y1 |
| 53 | + self.y1 = y0 |
| 54 | + |
| 55 | + return y0 |
| 56 | + |
| 57 | +def make_rc_lopass(sample_rate, freq): |
| 58 | + omega = 2 * math.pi * freq / sample_rate; |
| 59 | + term = 1 + 1/omega; |
| 60 | + return BiquadFilter(1/term, 0.0, 0.0, -1.0 + 1/term, 0.0); |
| 61 | + |
| 62 | +def make_butterworth(fs, fc, res_db=0): |
| 63 | + # 2nd-order Butterworth s-domain coefficients are: |
| 64 | + # |
| 65 | + # b0 = 1.0 b1 = 0 b2 = 0 |
| 66 | + # a0 = 1 a1 = sqrt(2) a2 = 1 |
| 67 | + # |
| 68 | + # by tweaking the a1 parameter, some resonance can be produced. |
| 69 | + |
| 70 | + res = 10.0 ** (-res_db / 10.0 / 2) |
| 71 | + ar = z_transform(1, 0, 0, 1, math.sqrt(2) * res, 1, fc, fs) |
| 72 | + return BiquadFilter(*ar) |
| 73 | + |
| 74 | +# observe: a and b are reversed here. To be absolutely clear: |
| 75 | +# a is the nominator and b is the denominator. :-/ |
| 76 | +def z_transform(a0, a1, a2, b0, b1, b2, fc, fs): |
| 77 | + # prewarp s-domain coefficients |
| 78 | + wp = 2.0 * fs * math.tan(math.pi * fc / fs) |
| 79 | + a2 /= wp * wp |
| 80 | + a1 /= wp |
| 81 | + b2 /= wp * wp |
| 82 | + b1 /= wp |
| 83 | + |
| 84 | + # compute bilinear transform and return it |
| 85 | + bd = 4 * b2 * fs * fs + 2 * b1 * fs + b0 |
| 86 | + return [ |
| 87 | + (4 * a2 * fs ** 2 + 2 * a1 * fs + a0) / bd, |
| 88 | + (2 * a0 - 8 * a2 * fs ** 2) / bd, |
| 89 | + (4 * a2 * fs ** 2 - 2 * a1 * fs + a0) / bd, |
| 90 | + (2 * b0 - 8 * b2 * fs ** 2) / bd, |
| 91 | + (4 * b2 * fs ** 2 - 2 * b1 * fs + b0) / bd, |
| 92 | + ] |
| 93 | + |
| 94 | +def quantize(x, bits, scale=False): |
| 95 | + x = list(x) |
| 96 | + fact = 2 ** bits |
| 97 | + |
| 98 | + # this adjusts range precisely between -65536 and 0 so that our bleps look right. |
| 99 | + # we should only do this if the table is integrated; in fact this code probably |
| 100 | + # belonged to the integration function. |
| 101 | + correction_factor = 1.0 |
| 102 | + if scale: |
| 103 | + correction_factor = x[-1] - x[0] |
| 104 | + |
| 105 | + for _ in range(len(x)): |
| 106 | + val = x[_] * fact / correction_factor; |
| 107 | + # correct rounding |
| 108 | + if val < 0: |
| 109 | + val = int(val - 0.5) |
| 110 | + else: |
| 111 | + val = int(val + 0.5) |
| 112 | + # leave scaled? |
| 113 | + if not scale: |
| 114 | + val /= float(fact) |
| 115 | + x[_] = val * -1 |
| 116 | + return x |
| 117 | + |
| 118 | +def lin2db(lin): |
| 119 | + return 20 * (math.log(lin) / math.log(10)) |
| 120 | + |
| 121 | +def print_spectrum(table, sample_rate): |
| 122 | + for _ in range(len(table) / 2): |
| 123 | + mag = lin2db(abs(table[_])) |
| 124 | + pha = math.atan2(table[_].real, table[_].imag) |
| 125 | + print "%s %s %s" % (float(_) / len(table) * sample_rate, mag, pha) |
| 126 | + |
| 127 | +def print_fir(table, format='gnuplot'): |
| 128 | + if format == 'gnuplot': |
| 129 | + for _ in range(len(table)): |
| 130 | + print "%s %s" % (_, table[_]) |
| 131 | + elif format == 'c': |
| 132 | + col = 0 |
| 133 | + print " {" |
| 134 | + for _ in range(len(table)): |
| 135 | + col += len(str(table[_])) + 1 |
| 136 | + if col >= 80: |
| 137 | + print |
| 138 | + col = 0 |
| 139 | + sys.stdout.write("%s," % table[_]) |
| 140 | + if col != 0: |
| 141 | + print |
| 142 | + print " }," |
| 143 | + |
| 144 | +def integrate(table): |
| 145 | + total = 0 |
| 146 | + for _ in table: |
| 147 | + total += _ |
| 148 | + startval = -total |
| 149 | + new = [] |
| 150 | + for _ in table: |
| 151 | + startval += _ |
| 152 | + new.append(startval) |
| 153 | + |
| 154 | + return new |
| 155 | + |
| 156 | +def run_filter(flt, table): |
| 157 | + flt.reset() |
| 158 | + |
| 159 | + # initialize filter to stable state |
| 160 | + for _ in range(10000): |
| 161 | + flt.filter(table[0]) |
| 162 | + |
| 163 | + # now run the filter |
| 164 | + newtable = [] |
| 165 | + for _ in range(len(table)): |
| 166 | + newtable.append(flt.filter(table[_])) |
| 167 | + |
| 168 | + return newtable, abs(table[-1] - flt.filter(table[-1])) |
| 169 | + |
| 170 | +AMIGA_PAL_CLOCK = 3546895 |
| 171 | +def main(): |
| 172 | + spectrum = len(sys.argv) > 1 |
| 173 | + |
| 174 | + # Because Amiga only has 84 dB SNR, the noise floor is low enough with -90 dB. |
| 175 | + # A500 model uses slightly lower-quality kaiser window to obtain slightly |
| 176 | + # steeper stopband attenuation. The fixed filters attenuates the sidelobes by |
| 177 | + # 12 dB, compensating for the worse performance of the kaiser window. |
| 178 | + |
| 179 | + # 21 kHz stopband is not fully attenuated by 22 kHz. If the sampling frequency |
| 180 | + # is 44.1 kHz, all frequencies above 22 kHz will alias over 20 kHz, thus inaudible. |
| 181 | + # The output should be aliasingless for 48 kHz sampling frequency. |
| 182 | + unfiltered_a500 = firwin(2048, 21000.0 / AMIGA_PAL_CLOCK * 2, window=('kaiser', 8.0)) |
| 183 | + unfiltered_a1200 = firwin(2048, 21000.0 / AMIGA_PAL_CLOCK * 2, window=('kaiser', 9.0)) |
| 184 | + # move filtering effects to start to allow IIRs more time to settle |
| 185 | + unfiltered_a500 = fir_minphase(unfiltered_a500) |
| 186 | + unfiltered_a1200 = fir_minphase(unfiltered_a1200) |
| 187 | + |
| 188 | + # make digital models for the filters on Amiga 500 and 1200. |
| 189 | + filter_fixed5khz = make_rc_lopass(AMIGA_PAL_CLOCK, 4900.0) |
| 190 | + # the leakage filter seems to reduce treble in both models a bit |
| 191 | + # the A500 filter seems to be well modelled only with a 4.9 kHz |
| 192 | + # filter although the component values would suggest 5 kHz filter. |
| 193 | + filter_leakage = make_rc_lopass(AMIGA_PAL_CLOCK, 32000.0) |
| 194 | + filter_led = make_butterworth(AMIGA_PAL_CLOCK, 3275.0, res_db=-0.70) |
| 195 | + |
| 196 | + # apply fixed filter to A500 |
| 197 | + amiga500_off, error500_off = run_filter(filter_fixed5khz, unfiltered_a500) |
| 198 | + # produce the filtered outputs |
| 199 | + amiga1200_off, error1200_off = run_filter(filter_leakage, unfiltered_a1200) |
| 200 | + |
| 201 | + # produce LED filters |
| 202 | + amiga500_on, error500_on = run_filter(filter_led, amiga500_off) |
| 203 | + amiga1200_on, error1200_on = run_filter(filter_led, amiga1200_off) |
| 204 | + |
| 205 | + # these values tell the error from truncating the run_filter() result at end. |
| 206 | + # print "error term magnitudes: %s" % map(lin2db, (error500_off, error500_on, error1200_off, error1200_on)) |
| 207 | + |
| 208 | + if not spectrum: |
| 209 | + # integrate to produce blep |
| 210 | + amiga500_off = integrate(amiga500_off) |
| 211 | + amiga500_on = integrate(amiga500_on) |
| 212 | + amiga1200_off = integrate(amiga1200_off) |
| 213 | + amiga1200_on = integrate(amiga1200_on) |
| 214 | + unfiltered_a1200 = integrate(unfiltered_a1200) |
| 215 | + |
| 216 | + # quantize and scale |
| 217 | + amiga500_off = quantize(amiga500_off, bits=17, scale=(not spectrum)) |
| 218 | + amiga500_on = quantize(amiga500_on, bits=17, scale=(not spectrum)) |
| 219 | + amiga1200_off = quantize(amiga1200_off, bits=17, scale=(not spectrum)) |
| 220 | + amiga1200_on = quantize(amiga1200_on, bits=17, scale=(not spectrum)) |
| 221 | + unfiltered_a1200 = quantize(unfiltered_a1200, bits=17, scale=(not spectrum)) |
| 222 | + |
| 223 | + if spectrum: |
| 224 | + spec = int(sys.argv[1]) |
| 225 | + if spec == -1: |
| 226 | + table = unfiltered_a1200 |
| 227 | + if spec == 0: |
| 228 | + table = amiga500_off |
| 229 | + if spec == 1: |
| 230 | + table = amiga500_on |
| 231 | + if spec == 2: |
| 232 | + table = amiga1200_off |
| 233 | + if spec == 3: |
| 234 | + table = amiga1200_on |
| 235 | + |
| 236 | + table = list(table); |
| 237 | + table += [0] * (16384 - len(table)) |
| 238 | + print_spectrum(fft(table), sample_rate=AMIGA_PAL_CLOCK) |
| 239 | + else: |
| 240 | + print " /*" |
| 241 | + print " * Table generated by contrib/sinc-integral.py." |
| 242 | + print " */" |
| 243 | + print |
| 244 | + print '#include "sinctable.h"' |
| 245 | + print |
| 246 | + print "/* tables are: a500 off, a500 on, a1200 off, a1200 on, vanilla. */" |
| 247 | + print "const int winsinc_integral[5][%d] = {" % len(unfiltered_a1200) |
| 248 | + print_fir(amiga500_off, format='c') |
| 249 | + print_fir(amiga500_on, format='c') |
| 250 | + print_fir(amiga1200_off, format='c') |
| 251 | + print_fir(amiga1200_on, format='c') |
| 252 | + print_fir(unfiltered_a1200, format='c') |
| 253 | + print "};" |
| 254 | + |
| 255 | +if __name__ == '__main__': |
| 256 | + main() |
| 257 | + |
0 commit comments