~ubuntu-branches/ubuntu/raring/vice/raring

« back to all changes in this revision

Viewing changes to src/resid-fp/sid.cc

  • Committer: Bazaar Package Importer
  • Author(s): Laszlo Boszormenyi (GCS)
  • Date: 2009-03-31 00:37:15 UTC
  • mfrom: (1.1.7 upstream) (9.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20090331003715-i5yisvcfv7mgz3eh
Tags: 2.1.dfsg-1
* New major upstream release (closes: #495937).
* Add desktop files (closes: #501181).

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
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
}