home *** CD-ROM | disk | FTP | other *** search
/ Large Pack of OldSkool DOS MOD Trackers / goattracker_2.68.zip / src / resid-fp / sid.cpp < prev    next >
C/C++ Source or Header  |  2009-01-03  |  29KB  |  937 lines
  1. //  ---------------------------------------------------------------------------
  2. //  This file is part of reSID, a MOS6581 SID emulator engine.
  3. //  Copyright (C) 2004  Dag Lem <resid@nimrod.no>
  4. //
  5. //  This program is free software; you can redistribute it and/or modify
  6. //  it under the terms of the GNU General Public License as published by
  7. //  the Free Software Foundation; either version 2 of the License, or
  8. //  (at your option) any later version.
  9. //
  10. //  This program is distributed in the hope that it will be useful,
  11. //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12. //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13. //  GNU General Public License for more details.
  14. //
  15. //  You should have received a copy of the GNU General Public License
  16. //  along with this program; if not, write to the Free Software
  17. //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18. //  ---------------------------------------------------------------------------
  19.  
  20. #include "sid.h"
  21.  
  22. #include <math.h>
  23.  
  24. extern float convolve(const float *a, const float *b, int n);
  25. extern float convolve_sse(const float *a, const float *b, int n);
  26.  
  27. enum host_cpu_feature {
  28.     HOST_CPU_MMX=1, HOST_CPU_SSE=2, HOST_CPU_SSE2=4, HOST_CPU_SSE3=8
  29. };
  30.  
  31. /* This code is appropriate for 32-bit and 64-bit x86 CPUs. */
  32. #if defined(__x86_64__) || defined(__i386__) || defined(_MSC_VER)
  33.  
  34. struct cpu_x86_regs_s {
  35.   unsigned int eax;
  36.   unsigned int ebx;
  37.   unsigned int ecx;
  38.   unsigned int edx;
  39. };
  40. typedef struct cpu_x86_regs_s cpu_x86_regs_t;
  41.  
  42. static cpu_x86_regs_t get_cpuid_regs(unsigned int index)
  43. {
  44.   cpu_x86_regs_t retval;
  45.  
  46. #if defined(_MSC_VER) /* MSVC assembly */
  47.   __asm {
  48.     mov eax, [index]
  49.     cpuid
  50.     mov [retval.eax], eax
  51.     mov [retval.ebx], ebx
  52.     mov [retval.ecx], ecx
  53.     mov [retval.edx], edx
  54.   }
  55. #else /* GNU assembly */
  56.   asm("movl %4, %%eax; cpuid; movl %%eax, %0; movl %%ebx, %1; movl %%ecx, %2; movl %%edx, %3;"
  57.       : "=m" (retval.eax),
  58.         "=m" (retval.ebx),
  59.         "=m" (retval.ecx),
  60.         "=m" (retval.edx)
  61.       : "r"  (index)
  62.       : "eax", "ebx", "ecx", "edx");
  63. #endif
  64.  
  65.   return retval;
  66. }
  67.  
  68. static int host_cpu_features_by_cpuid(void)
  69. {
  70.   cpu_x86_regs_t regs = get_cpuid_regs(1);
  71.  
  72.   int features = 0;
  73.   if (regs.edx & (1 << 23))
  74.     features |= HOST_CPU_MMX;
  75.   if (regs.edx & (1 << 25))
  76.     features |= HOST_CPU_SSE;
  77.   if (regs.edx & (1 << 26))
  78.     features |= HOST_CPU_SSE2;
  79.   if (regs.ecx & (1 << 0))
  80.     features |= HOST_CPU_SSE3;
  81.  
  82.   return features;
  83. }
  84.  
  85. static int host_cpu_features(void)
  86. {
  87.   static int features = 0;
  88.   static int features_detected = 0;
  89. /* 32-bit only */
  90. #if defined(__i386__) || (defined(_MSC_VER) && defined(_WIN32))
  91.   unsigned long temp1, temp2;
  92. #endif
  93.  
  94.   if (features_detected)
  95.     return features;
  96.   features_detected = 1;
  97.  
  98. #if defined(_MSC_VER) && defined(_WIN32) /* MSVC compatible assembly appropriate for 32-bit Windows */
  99.   /* see if we are dealing with a cpu that has the cpuid instruction */
  100.   __asm {
  101.     pushf
  102.     pop eax
  103.     mov [temp1], eax
  104.     xor eax, 0x200000
  105.     push eax
  106.     popf
  107.     pushf
  108.     pop eax
  109.     mov [temp2], eax
  110.     push [temp1]
  111.     popf
  112.   }
  113. #endif
  114. #if defined(__i386__) /* GNU assembly */
  115.   asm("pushfl; popl %%eax; movl %%eax, %0; xorl $0x200000, %%eax; pushl %%eax; popfl; pushfl; popl %%eax; movl %%eax, %1; pushl %0; popfl "
  116.       : "=r" (temp1),
  117.       "=r" (temp2)
  118.       :
  119.       : "eax");
  120. #endif
  121. #if defined(__i386__) || (defined(_MSC_VER) && defined(_WIN32))
  122.   temp1 &= 0x200000;
  123.   temp2 &= 0x200000;
  124.   if (temp1 == temp2) {
  125.     /* no cpuid support, so we can't test for SSE availability -> false */
  126.     return 0;
  127.   }
  128. #endif
  129.  
  130.   /* find the highest supported cpuid function, returned in %eax */
  131.   if (get_cpuid_regs(0).eax < 1) {
  132.     /* no cpuid 1 function, we can't test for features -> no features */
  133.     return 0;
  134.   }
  135.  
  136.   features = host_cpu_features_by_cpuid();
  137.   return features;
  138. }
  139.  
  140. #else /* !__x86_64__ && !__i386__ && !_MSC_VER */
  141. static int host_cpu_features(void)
  142. {
  143.   return 0;
  144. }
  145. #endif
  146.  
  147. float SIDFP::kinked_dac(const int x, const float nonlinearity, const int max)
  148. {
  149.     float value = 0.f;
  150.     
  151.     int bit = 1;
  152.     float weight = 1.f;
  153.     const float dir = 2.0f * nonlinearity;
  154.     for (int i = 0; i < max; i ++) {
  155.         if (x & bit)
  156.             value += weight;
  157.         bit <<= 1;
  158.         weight *= dir;
  159.     }
  160.  
  161.     return value / (weight / nonlinearity) * (1 << max);
  162. }
  163.  
  164. // ----------------------------------------------------------------------------
  165. // Constructor.
  166. // ----------------------------------------------------------------------------
  167. SIDFP::SIDFP()
  168. {
  169. #if (RESID_USE_SSE==1)
  170.   can_use_sse = (host_cpu_features() & HOST_CPU_SSE) != 0;
  171. #else
  172.   can_use_sse = false;
  173. #endif
  174.  
  175.   // Initialize pointers.
  176.   sample = 0;
  177.   fir = 0;
  178.  
  179.   voice[0].set_sync_source(&voice[2]);
  180.   voice[1].set_sync_source(&voice[0]);
  181.   voice[2].set_sync_source(&voice[1]);
  182.  
  183.   set_sampling_parameters(985248, SAMPLE_INTERPOLATE, 44100);
  184.  
  185.   bus_value = 0;
  186.   bus_value_ttl = 0;
  187.  
  188.   input(0);
  189. }
  190.  
  191.  
  192. // ----------------------------------------------------------------------------
  193. // Destructor.
  194. // ----------------------------------------------------------------------------
  195. SIDFP::~SIDFP()
  196. {
  197.   delete[] sample;
  198.   delete[] fir;
  199. }
  200.  
  201.  
  202. // ----------------------------------------------------------------------------
  203. // Set chip model.
  204. // ----------------------------------------------------------------------------
  205. void SIDFP::set_chip_model(chip_model model)
  206. {
  207.   for (int i = 0; i < 3; i++) {
  208.     voice[i].set_chip_model(model);
  209.   }
  210.  
  211.   filter.set_chip_model(model);
  212.   extfilt.set_chip_model(model);
  213. }
  214.  
  215. /* nonlinear DAC support, set 1 for 8580 / no effect, about 0.96 otherwise */
  216. void SIDFP::set_voice_nonlinearity(float nl)
  217. {
  218.   for (int i = 0; i < 3; i++) {
  219.     voice[i].set_nonlinearity(nl);
  220.   }
  221. }
  222.  
  223. // ----------------------------------------------------------------------------
  224. // SID reset.
  225. // ----------------------------------------------------------------------------
  226. void SIDFP::reset()
  227. {
  228.   for (int i = 0; i < 3; i++) {
  229.     voice[i].reset();
  230.   }
  231.   filter.reset();
  232.   extfilt.reset();
  233.  
  234.   bus_value = 0;
  235.   bus_value_ttl = 0;
  236. }
  237.  
  238.  
  239. // ----------------------------------------------------------------------------
  240. // Write 16-bit sample to audio input.
  241. // NB! The caller is responsible for keeping the value within 16 bits.
  242. // Note that to mix in an external audio signal, the signal should be
  243. // resampled to 1MHz first to avoid sampling noise.
  244. // ----------------------------------------------------------------------------
  245. void SIDFP::input(int sample)
  246. {
  247.   // Voice outputs are 20 bits. Scale up to match three voices in order
  248.   // to facilitate simulation of the MOS8580 "digi boost" hardware hack.
  249.   ext_in = (float) ( (sample << 4) * 3 );
  250. }
  251.  
  252. float SIDFP::output()
  253. {
  254.   const float range = 1 << 15;
  255.   return extfilt.output() / (4095.f * 255.f * 3.f * 1.5f / range);
  256. }
  257.  
  258. // ----------------------------------------------------------------------------
  259. // Read registers.
  260. //
  261. // Reading a write only register returns the last byte written to any SID
  262. // register. The individual bits in this value start to fade down towards
  263. // zero after a few cycles. All bits reach zero within approximately
  264. // $2000 - $4000 cycles.
  265. // It has been claimed that this fading happens in an orderly fashion, however
  266. // sampling of write only registers reveals that this is not the case.
  267. // NB! This is not correctly modeled.
  268. // The actual use of write only registers has largely been made in the belief
  269. // that all SID registers are readable. To support this belief the read
  270. // would have to be done immediately after a write to the same register
  271. // (remember that an intermediate write to another register would yield that
  272. // value instead). With this in mind we return the last value written to
  273. // any SID register for $2000 cycles without modeling the bit fading.
  274. // ----------------------------------------------------------------------------
  275. reg8 SIDFP::read(reg8 offset)
  276. {
  277.   switch (offset) {
  278.   case 0x19:
  279.     return potx.readPOT();
  280.   case 0x1a:
  281.     return poty.readPOT();
  282.   case 0x1b:
  283.     return voice[2].wave.readOSC();
  284.   case 0x1c:
  285.     return voice[2].envelope.readENV();
  286.   default:
  287.     return bus_value;
  288.   }
  289. }
  290.  
  291.  
  292. // ----------------------------------------------------------------------------
  293. // Write registers.
  294. // ----------------------------------------------------------------------------
  295. void SIDFP::write(reg8 offset, reg8 value)
  296. {
  297.   bus_value = value;
  298.   bus_value_ttl = 0x4000;
  299.  
  300.   switch (offset) {
  301.   case 0x00:
  302.     voice[0].wave.writeFREQ_LO(value);
  303.     break;
  304.   case 0x01:
  305.     voice[0].wave.writeFREQ_HI(value);
  306.     break;
  307.   case 0x02:
  308.     voice[0].wave.writePW_LO(value);
  309.     break;
  310.   case 0x03:
  311.     voice[0].wave.writePW_HI(value);
  312.     break;
  313.   case 0x04:
  314.     voice[0].writeCONTROL_REG(value);
  315.     break;
  316.   case 0x05:
  317.     voice[0].envelope.writeATTACK_DECAY(value);
  318.     break;
  319.   case 0x06:
  320.     voice[0].envelope.writeSUSTAIN_RELEASE(value);
  321.     break;
  322.   case 0x07:
  323.     voice[1].wave.writeFREQ_LO(value);
  324.     break;
  325.   case 0x08:
  326.     voice[1].wave.writeFREQ_HI(value);
  327.     break;
  328.   case 0x09:
  329.     voice[1].wave.writePW_LO(value);
  330.     break;
  331.   case 0x0a:
  332.     voice[1].wave.writePW_HI(value);
  333.     break;
  334.   case 0x0b:
  335.     voice[1].writeCONTROL_REG(value);
  336.     break;
  337.   case 0x0c:
  338.     voice[1].envelope.writeATTACK_DECAY(value);
  339.     break;
  340.   case 0x0d:
  341.     voice[1].envelope.writeSUSTAIN_RELEASE(value);
  342.     break;
  343.   case 0x0e:
  344.     voice[2].wave.writeFREQ_LO(value);
  345.     break;
  346.   case 0x0f:
  347.     voice[2].wave.writeFREQ_HI(value);
  348.     break;
  349.   case 0x10:
  350.     voice[2].wave.writePW_LO(value);
  351.     break;
  352.   case 0x11:
  353.     voice[2].wave.writePW_HI(value);
  354.     break;
  355.   case 0x12:
  356.     voice[2].writeCONTROL_REG(value);
  357.     break;
  358.   case 0x13:
  359.     voice[2].envelope.writeATTACK_DECAY(value);
  360.     break;
  361.   case 0x14:
  362.     voice[2].envelope.writeSUSTAIN_RELEASE(value);
  363.     break;
  364.   case 0x15:
  365.     filter.writeFC_LO(value);
  366.     break;
  367.   case 0x16:
  368.     filter.writeFC_HI(value);
  369.     break;
  370.   case 0x17:
  371.     filter.writeRES_FILT(value);
  372.     break;
  373.   case 0x18:
  374.     filter.writeMODE_VOL(value);
  375.     break;
  376.   default:
  377.     break;
  378.   }
  379. }
  380.  
  381.  
  382. // ----------------------------------------------------------------------------
  383. // Constructor.
  384. // ----------------------------------------------------------------------------
  385. SIDFP::State::State()
  386. {
  387.   int i;
  388.  
  389.   for (i = 0; i < 0x20; i++) {
  390.     sid_register[i] = 0;
  391.   }
  392.  
  393.   bus_value = 0;
  394.   bus_value_ttl = 0;
  395.  
  396.   for (i = 0; i < 3; i++) {
  397.     accumulator[i] = 0;
  398.     shift_register[i] = 0x7ffff8;
  399.     rate_counter[i] = 0;
  400.     rate_counter_period[i] = 9;
  401.     exponential_counter[i] = 0;
  402.     exponential_counter_period[i] = 1;
  403.     envelope_counter[i] = 0;
  404.     envelope_state[i] = EnvelopeGeneratorFP::RELEASE;
  405.     hold_zero[i] = true;
  406.   }
  407. }
  408.  
  409.  
  410. // ----------------------------------------------------------------------------
  411. // Read state.
  412. // ----------------------------------------------------------------------------
  413. SIDFP::State SIDFP::read_state()
  414. {
  415.   State state;
  416.   int i, j;
  417.  
  418.   for (i = 0, j = 0; i < 3; i++, j += 7) {
  419.     WaveformGeneratorFP& wave = voice[i].wave;
  420.     EnvelopeGeneratorFP& envelope = voice[i].envelope;
  421.     state.sid_register[j + 0] = wave.freq & 0xff;
  422.     state.sid_register[j + 1] = wave.freq >> 8;
  423.     state.sid_register[j + 2] = wave.pw & 0xff;
  424.     state.sid_register[j + 3] = wave.pw >> 8;
  425.     state.sid_register[j + 4] =
  426.       (wave.waveform << 4)
  427.       | (wave.test ? 0x08 : 0)
  428.       | (wave.ring_mod ? 0x04 : 0)
  429.       | (wave.sync ? 0x02 : 0)
  430.       | (envelope.gate ? 0x01 : 0);
  431.     state.sid_register[j + 5] = (envelope.attack << 4) | envelope.decay;
  432.     state.sid_register[j + 6] = (envelope.sustain << 4) | envelope.release;
  433.   }
  434.  
  435.   state.sid_register[j++] = filter.fc & 0x007;
  436.   state.sid_register[j++] = filter.fc >> 3;
  437.   state.sid_register[j++] = (filter.res << 4) | filter.filt;
  438.   state.sid_register[j++] =
  439.     (filter.voice3off ? 0x80 : 0)
  440.     | (filter.hp_bp_lp << 4)
  441.     | filter.vol;
  442.  
  443.   // These registers are superfluous, but included for completeness.
  444.   for (; j < 0x1d; j++) {
  445.     state.sid_register[j] = read(j);
  446.   }
  447.   for (; j < 0x20; j++) {
  448.     state.sid_register[j] = 0;
  449.   }
  450.  
  451.   state.bus_value = bus_value;
  452.   state.bus_value_ttl = bus_value_ttl;
  453.  
  454.   for (i = 0; i < 3; i++) {
  455.     state.accumulator[i] = voice[i].wave.accumulator;
  456.     state.shift_register[i] = voice[i].wave.shift_register;
  457.     state.rate_counter[i] = voice[i].envelope.rate_counter;
  458.     state.rate_counter_period[i] = voice[i].envelope.rate_period;
  459.     state.exponential_counter[i] = voice[i].envelope.exponential_counter;
  460.     state.exponential_counter_period[i] = voice[i].envelope.exponential_counter_period;
  461.     state.envelope_counter[i] = voice[i].envelope.envelope_counter;
  462.     state.envelope_state[i] = voice[i].envelope.state;
  463.     state.hold_zero[i] = voice[i].envelope.hold_zero;
  464.   }
  465.  
  466.   return state;
  467. }
  468.  
  469.  
  470. // ----------------------------------------------------------------------------
  471. // Write state.
  472. // ----------------------------------------------------------------------------
  473. void SIDFP::write_state(const State& state)
  474. {
  475.   int i;
  476.  
  477.   for (i = 0; i <= 0x18; i++) {
  478.     write(i, state.sid_register[i]);
  479.   }
  480.  
  481.   bus_value = state.bus_value;
  482.   bus_value_ttl = state.bus_value_ttl;
  483.  
  484.   for (i = 0; i < 3; i++) {
  485.     voice[i].wave.accumulator = state.accumulator[i];
  486.     voice[i].wave.shift_register = state.shift_register[i];
  487.     voice[i].envelope.rate_counter = state.rate_counter[i];
  488.     voice[i].envelope.rate_period = state.rate_counter_period[i];
  489.     voice[i].envelope.exponential_counter = state.exponential_counter[i];
  490.     voice[i].envelope.exponential_counter_period = state.exponential_counter_period[i];
  491.     voice[i].envelope.envelope_counter = state.envelope_counter[i];
  492.     voice[i].envelope.state = state.envelope_state[i];
  493.     voice[i].envelope.hold_zero = state.hold_zero[i];
  494.   }
  495. }
  496.  
  497.  
  498. // ----------------------------------------------------------------------------
  499. // Enable filter.
  500. // ----------------------------------------------------------------------------
  501. void SIDFP::enable_filter(bool enable)
  502. {
  503.   filter.enable_filter(enable);
  504. }
  505.  
  506.  
  507. // ----------------------------------------------------------------------------
  508. // Enable external filter.
  509. // ----------------------------------------------------------------------------
  510. void SIDFP::enable_external_filter(bool enable)
  511. {
  512.   extfilt.enable_filter(enable);
  513. }
  514.  
  515.  
  516. // ----------------------------------------------------------------------------
  517. // I0() computes the 0th order modified Bessel function of the first kind.
  518. // This function is originally from resample-1.5/filterkit.c by J. O. Smith.
  519. // ----------------------------------------------------------------------------
  520. double SIDFP::I0(double x)
  521. {
  522.   // Max error acceptable in I0 could be 1e-6, which gives that 96 dB already.
  523.   // I'm overspecify these errors to get a beautiful FFT dump of the FIR.
  524.   const double I0e = 1e-10;
  525.  
  526.   double sum, u, halfx, temp;
  527.   int n;
  528.  
  529.   sum = u = n = 1;
  530.   halfx = x/2.0;
  531.  
  532.   do {
  533.     temp = halfx/n++;
  534.     u *= temp*temp;
  535.     sum += u;
  536.   } while (u >= I0e*sum);
  537.  
  538.   return sum;
  539. }
  540.  
  541.  
  542. // ----------------------------------------------------------------------------
  543. // Setting of SID sampling parameters.
  544. //
  545. // Use a clock freqency of 985248Hz for PAL C64, 1022730Hz for NTSC C64.
  546. // The default end of passband frequency is pass_freq = 0.9*sample_freq/2
  547. // for sample frequencies up to ~ 44.1kHz, and 20kHz for higher sample
  548. // frequencies.
  549. //
  550. // For resampling, the ratio between the clock frequency and the sample
  551. // frequency is limited as follows:
  552. //   125*clock_freq/sample_freq < 16384
  553. // E.g. provided a clock frequency of ~ 1MHz, the sample frequency can not
  554. // be set lower than ~ 8kHz. A lower sample frequency would make the
  555. // resampling code overfill its 16k sample ring buffer.
  556. // 
  557. // The end of passband frequency is also limited:
  558. //   pass_freq <= 0.9*sample_freq/2
  559.  
  560. // E.g. for a 44.1kHz sampling rate the end of passband frequency is limited
  561. // to slightly below 20kHz. This constraint ensures that the FIR table is
  562. // not overfilled.
  563. // ----------------------------------------------------------------------------
  564. bool SIDFP::set_sampling_parameters(float clock_freq, sampling_method method,
  565.                                   float sample_freq, float pass_freq)
  566. {
  567.   clock_frequency = clock_freq;
  568.   sampling = method;
  569.  
  570.   filter.set_clock_frequency(clock_freq);
  571.   extfilt.set_clock_frequency(clock_freq);
  572.   adjust_sampling_frequency(sample_freq);
  573.  
  574.   sample_offset = 0;
  575.   sample_prev = 0;
  576.  
  577.   // FIR initialization is only necessary for resampling.
  578.   if (method != SAMPLE_RESAMPLE_INTERPOLATE)
  579.   {
  580.     delete[] sample;
  581.     delete[] fir;
  582.     sample = 0;
  583.     fir = 0;
  584.     return true;
  585.   }
  586.   
  587.   const int bits = 16;
  588.  
  589.   if (pass_freq > 20000)
  590.     pass_freq = 20000;  
  591.   if (2*pass_freq/sample_freq > 0.9)
  592.     pass_freq = 0.9f*sample_freq/2;
  593.  
  594.   // 16 bits -> -96dB stopband attenuation.
  595.   const double A = -20*log10(1.0/(1 << bits));
  596.  
  597.   // For calculation of beta and N see the reference for the kaiserord
  598.   // function in the MATLAB Signal Processing Toolbox:
  599.   // http://www.mathworks.com/access/helpdesk/help/toolbox/signal/kaiserord.html
  600.   const double beta = 0.1102*(A - 8.7);
  601.   const double I0beta = I0(beta);
  602.  
  603.   double f_samples_per_cycle = sample_freq/clock_freq;
  604.   double f_cycles_per_sample = clock_freq/sample_freq;
  605.  
  606.   /* This code utilizes the fact that aliasing back to 20 kHz from
  607.    * sample_freq/2 is inaudible. This allows us to define a passband
  608.    * wider than normally. We might also consider aliasing back to pass_freq,
  609.    * but as this can be less than 20 kHz, it might become audible... */
  610.   double aliasing_allowance = sample_freq / 2 - 20000;
  611.   if (aliasing_allowance < 0)
  612.     aliasing_allowance = 0;
  613.  
  614.   double transition_bandwidth = sample_freq/2 - pass_freq + aliasing_allowance;
  615.   {
  616.     /* Filter order according to Kaiser's paper. */
  617.  
  618.     int N = (int) ((A - 7.95)/(2 * M_PI * 2.285 * transition_bandwidth/sample_freq) + 0.5);
  619.     N += N & 1;
  620.  
  621.     // The filter length is equal to the filter order + 1.
  622.     // The filter length must be an odd number (sinc is symmetric about x = 0).
  623.     fir_N = int(N*f_cycles_per_sample) + 1;
  624.     fir_N |= 1;
  625.  
  626.     // Check whether the sample ring buffer would overfill.
  627.     if (fir_N > RINGSIZE - 1)
  628.       return false;
  629.  
  630.     /* Error is bound by 1.234 / L^2 */
  631.     fir_RES = (int) (sqrt(1.234 * (1 << bits)) / f_cycles_per_sample + 0.5);
  632.   }
  633.  
  634.   // Allocate memory for FIR tables.
  635.   delete[] fir;
  636.   fir = new float[fir_N*fir_RES];
  637.  
  638.   // The cutoff frequency is midway through the transition band.
  639.   double wc = (pass_freq + transition_bandwidth/2) / sample_freq * M_PI * 2;
  640.  
  641.   // Calculate fir_RES FIR tables for linear interpolation.
  642.   for (int i = 0; i < fir_RES; i++) {
  643.     double j_offset = double(i)/fir_RES;
  644.     // Calculate FIR table. This is the sinc function, weighted by the
  645.     // Kaiser window.
  646.     for (int j = 0; j < fir_N; j ++) {
  647.       double jx = j - fir_N/2. - j_offset;
  648.       double wt = wc*jx/f_cycles_per_sample;
  649.       double temp = jx/(fir_N/2);
  650.       double Kaiser =
  651.         fabs(temp) <= 1 ? I0(beta*sqrt(1 - temp*temp))/I0beta : 0;
  652.       double sincwt =
  653.         fabs(wt) >= 1e-8 ? sin(wt)/wt : 1;
  654.       fir[i * fir_N + j] = (float) (f_samples_per_cycle*wc/M_PI*sincwt*Kaiser);
  655.     }
  656.   }
  657.  
  658.   // Allocate sample buffer.
  659.   if (!sample) {
  660.     sample = new float[RINGSIZE*2];
  661.   }
  662.   // Clear sample buffer.
  663.   for (int j = 0; j < RINGSIZE*2; j++) {
  664.     sample[j] = 0;
  665.   }
  666.   sample_index = 0;
  667.  
  668.   return true;
  669. }
  670.  
  671. // ----------------------------------------------------------------------------
  672. // Adjustment of SID sampling frequency.
  673. //
  674. // In some applications, e.g. a C64 emulator, it can be desirable to
  675. // synchronize sound with a timer source. This is supported by adjustment of
  676. // the SID sampling frequency.
  677. //
  678. // NB! Adjustment of the sampling frequency may lead to noticeable shifts in
  679. // frequency, and should only be used for interactive applications. Note also
  680. // that any adjustment of the sampling frequency will change the
  681. // characteristics of the resampling filter, since the filter is not rebuilt.
  682. // ----------------------------------------------------------------------------
  683. void SIDFP::adjust_sampling_frequency(float sample_freq)
  684. {
  685.   cycles_per_sample = clock_frequency/sample_freq;
  686. }
  687.  
  688. void SIDFP::age_bus_value(cycle_count n) {
  689.   if (bus_value_ttl != 0) {
  690.     bus_value_ttl -= n;
  691.     if (bus_value_ttl <= 0) {
  692.         bus_value = 0;
  693.         bus_value_ttl = 0;
  694.     }
  695.   }
  696. }
  697.  
  698. // ----------------------------------------------------------------------------
  699. // SID clocking - 1 cycle.
  700. // ----------------------------------------------------------------------------
  701. void SIDFP::clock()
  702. {
  703.   int i;
  704.  
  705.   // Clock amplitude modulators.
  706.   for (i = 0; i < 3; i++) {
  707.     voice[i].envelope.clock();
  708.   }
  709.  
  710.   // Clock oscillators.
  711.   for (i = 0; i < 3; i++) {
  712.     voice[i].wave.clock();
  713.   }
  714.  
  715.   // Synchronize oscillators.
  716.   for (i = 0; i < 3; i++) {
  717.     voice[i].wave.synchronize();
  718.   }
  719.  
  720.   // Clock filter.
  721.   extfilt.clock(filter.clock(voice[0].output(), voice[1].output(), voice[2].output(), ext_in));
  722. }
  723.  
  724. // ----------------------------------------------------------------------------
  725. // SID clocking with audio sampling.
  726. // Fixpoint arithmetics is used.
  727. //
  728. // The example below shows how to clock the SID a specified amount of cycles
  729. // while producing audio output:
  730. //
  731. // while (delta_t) {
  732. //   bufindex += sid.clock(delta_t, buf + bufindex, buflength - bufindex);
  733. //   write(dsp, buf, bufindex*2);
  734. //   bufindex = 0;
  735. // }
  736. // 
  737. // ----------------------------------------------------------------------------
  738. int SIDFP::clock(cycle_count& delta_t, short* buf, int n, int interleave)
  739. {
  740.   /* XXX I assume n is generally large enough for delta_t here... */
  741.   age_bus_value(delta_t);
  742.   int res;
  743.   switch (sampling) {
  744.   default:
  745.   case SAMPLE_INTERPOLATE:
  746.     res = clock_interpolate(delta_t, buf, n, interleave);
  747.     break;
  748.   case SAMPLE_RESAMPLE_INTERPOLATE:
  749.     res = clock_resample_interpolate(delta_t, buf, n, interleave);
  750.     break;
  751.   }
  752.  
  753.   filter.nuke_denormals();
  754.   extfilt.nuke_denormals();
  755.  
  756.   return res;
  757. }
  758.  
  759. // ----------------------------------------------------------------------------
  760. // SID clocking with audio sampling - cycle based with linear sample
  761. // interpolation.
  762. //
  763. // Here the chip is clocked every cycle. This yields higher quality
  764. // sound since the samples are linearly interpolated, and since the
  765. // external filter attenuates frequencies above 16kHz, thus reducing
  766. // sampling noise.
  767. // ----------------------------------------------------------------------------
  768. RESID_INLINE
  769. int SIDFP::clock_interpolate(cycle_count& delta_t, short* buf, int n,
  770.                            int interleave)
  771. {
  772.   int s = 0;
  773.   int i;
  774.  
  775.   for (;;) {
  776.     float next_sample_offset = sample_offset + cycles_per_sample;
  777.     int delta_t_sample = (int) next_sample_offset;
  778.     if (delta_t_sample > delta_t) {
  779.       break;
  780.     }
  781.     if (s >= n) {
  782.       return s;
  783.     }
  784.     for (i = 0; i < delta_t_sample - 1; i++) {
  785.       clock();
  786.     }
  787.     if (i < delta_t_sample) {
  788.       sample_prev = output();
  789.       clock();
  790.     }
  791.  
  792.     delta_t -= delta_t_sample;
  793.     sample_offset = next_sample_offset - delta_t_sample;
  794.  
  795.     float sample_now = output();
  796.     int v = (int)(sample_prev + (sample_offset * (sample_now - sample_prev)));
  797.     // Saturated arithmetics to guard against 16 bit sample overflow.
  798.     const int half = 1 << 15;
  799.     if (v >= half) {
  800.       v = half - 1;
  801.     }
  802.     else if (v < -half) {
  803.       v = -half;
  804.     }
  805.     buf[s++*interleave] = v;
  806.  
  807.     sample_prev = sample_now;
  808.   }
  809.  
  810.   for (i = 0; i < delta_t - 1; i++) {
  811.     clock();
  812.   }
  813.   if (i < delta_t) {
  814.     sample_prev = output();
  815.     clock();
  816.   }
  817.   sample_offset -= delta_t;
  818.   delta_t = 0;
  819.   return s;
  820. }
  821.  
  822. // ----------------------------------------------------------------------------
  823. // SID clocking with audio sampling - cycle based with audio resampling.
  824. //
  825. // This is the theoretically correct (and computationally intensive) audio
  826. // sample generation. The samples are generated by resampling to the specified
  827. // sampling frequency. The work rate is inversely proportional to the
  828. // percentage of the bandwidth allocated to the filter transition band.
  829. //
  830. // This implementation is based on the paper "A Flexible Sampling-Rate
  831. // Conversion Method", by J. O. Smith and P. Gosset, or rather on the
  832. // expanded tutorial on the "Digital Audio Resampling Home Page":
  833. // http://www-ccrma.stanford.edu/~jos/resample/
  834. //
  835. // By building shifted FIR tables with samples according to the
  836. // sampling frequency, this implementation dramatically reduces the
  837. // computational effort in the filter convolutions, without any loss
  838. // of accuracy. The filter convolutions are also vectorizable on
  839. // current hardware.
  840. //
  841. // Further possible optimizations are:
  842. // * An equiripple filter design could yield a lower filter order, see
  843. //   http://www.mwrf.com/Articles/ArticleID/7229/7229.html
  844. // * The Convolution Theorem could be used to bring the complexity of
  845. //   convolution down from O(n*n) to O(n*log(n)) using the Fast Fourier
  846. //   Transform, see http://en.wikipedia.org/wiki/Convolution_theorem
  847. // * Simply resampling in two steps can also yield computational
  848. //   savings, since the transition band will be wider in the first step
  849. //   and the required filter order is thus lower in this step.
  850. //   Laurent Ganier has found the optimal intermediate sampling frequency
  851. //   to be (via derivation of sum of two steps):
  852. //     2 * pass_freq + sqrt [ 2 * pass_freq * orig_sample_freq
  853. //       * (dest_sample_freq - 2 * pass_freq) / dest_sample_freq ]
  854. //
  855. // NB! the result of right shifting negative numbers is really
  856. // implementation dependent in the C++ standard.
  857. // ----------------------------------------------------------------------------
  858. RESID_INLINE
  859. int SIDFP::clock_resample_interpolate(cycle_count& delta_t, short* buf, int n,
  860.                                     int interleave)
  861. {
  862.   int s = 0;
  863.  
  864.   for (;;) {
  865.     float next_sample_offset = sample_offset + cycles_per_sample;
  866.     /* full clocks left to next sample */
  867.     int delta_t_sample = (int) next_sample_offset;
  868.     if (delta_t_sample > delta_t || s >= n)
  869.       break;
  870.  
  871.     /* clock forward delta_t_sample samples */
  872.     for (int i = 0; i < delta_t_sample; i++) {
  873.       clock();
  874.       sample[sample_index] = sample[sample_index + RINGSIZE] = output();
  875.       ++ sample_index;
  876.       sample_index &= RINGSIZE - 1;
  877.     }
  878.     delta_t -= delta_t_sample;
  879.  
  880.     /* Phase of the sample in terms of clock, [0 .. 1[. */
  881.     sample_offset = next_sample_offset - (float) delta_t_sample;
  882.  
  883.     /* find the first of the nearest fir tables close to the phase */
  884.     float fir_offset_rmd = sample_offset * fir_RES;
  885.     int fir_offset = (int) fir_offset_rmd;
  886.     /* [0 .. 1[ */
  887.     fir_offset_rmd -= (float) fir_offset;
  888.  
  889.     /* find fir_N most recent samples, plus one extra in case the FIR wraps. */
  890.     float* sample_start = sample + sample_index - fir_N + RINGSIZE - 1;
  891.  
  892.     float v1 =
  893. #if (RESID_USE_SSE==1)
  894.       can_use_sse ? convolve_sse(sample_start, fir + fir_offset*fir_N, fir_N) :
  895. #endif
  896.         convolve(sample_start, fir + fir_offset*fir_N, fir_N);
  897.  
  898.     // Use next FIR table, wrap around to first FIR table using
  899.     // previous sample.
  900.     if (++ fir_offset == fir_RES) {
  901.       fir_offset = 0;
  902.       ++ sample_start;
  903.     }
  904.     float v2 =
  905. #if (RESID_USE_SSE==1)
  906.       can_use_sse ? convolve_sse(sample_start, fir + fir_offset*fir_N, fir_N) :
  907. #endif
  908.         convolve(sample_start, fir + fir_offset*fir_N, fir_N);
  909.  
  910.     // Linear interpolation between the sinc tables yields good approximation
  911.     // for the exact value.
  912.     int v = (int) (v1 + fir_offset_rmd * (v2 - v1));
  913.  
  914.     // Saturated arithmetics to guard against 16 bit sample overflow.
  915.     const int half = 1 << 15;
  916.     if (v >= half) {
  917.       v = half - 1;
  918.     }
  919.     else if (v < -half) {
  920.       v = -half;
  921.     }
  922.  
  923.     buf[s ++ * interleave] = v;
  924.   }
  925.  
  926.   /* clock forward delta_t samples */
  927.   for (int i = 0; i < delta_t; i++) {
  928.     clock();
  929.     sample[sample_index] = sample[sample_index + RINGSIZE] = output();
  930.     ++ sample_index;
  931.     sample_index &= RINGSIZE - 1;
  932.   }
  933.   sample_offset -= (float) delta_t;
  934.   delta_t = 0;
  935.   return s;
  936. }
  937.