Trigger on down-flanks instead of up-flanks.
[c64tapwav] / decode.cpp
1 #include <stdio.h>
2 #include <string.h>
3 #include <math.h>
4 #include <unistd.h>
5 #include <assert.h>
6 #include <vector>
7 #include <algorithm>
8
9 #include "interpolate.h"
10
11 #define BUFSIZE 4096
12 #define HYSTERESIS_LIMIT 3000
13 #define SAMPLE_RATE 44100
14 #define C64_FREQUENCY 985248
15 #define TAP_RESOLUTION 8
16
17 #define SYNC_PULSE_START 1000
18 #define SYNC_PULSE_END 15000
19 #define SYNC_PULSE_LENGTH 378.0
20 #define SYNC_TEST_TOLERANCE 1.10
21
22 struct tap_header {
23         char identifier[12];
24         char version;
25         char reserved[3];
26         unsigned int data_len;
27 };
28
29 // between [x,x+1]
30 double find_zerocrossing(const std::vector<short> &pcm, int x)
31 {
32         if (pcm[x] == 0) {
33                 return x;
34         }
35         if (pcm[x + 1] == 0) {
36                 return x + 1;
37         }
38
39         assert(pcm[x + 1] < 0);
40         assert(pcm[x] > 0);
41
42         double upper = x;
43         double lower = x + 1;
44         while (upper - lower > 1e-6) {
45                 double mid = 0.5f * (upper + lower);
46                 if (lanczos_interpolate(pcm, mid) > 0) {
47                         upper = mid;
48                 } else {
49                         lower = mid;
50                 }
51         }
52
53         return 0.5f * (upper + lower);
54 }
55
56 struct pulse {
57         double time;  // in seconds from start
58         double len;   // in seconds
59 };
60         
61 int main(int argc, char **argv)
62 {
63         std::vector<short> pcm;
64
65         while (!feof(stdin)) {
66                 short buf[BUFSIZE];
67                 ssize_t ret = fread(buf, 2, BUFSIZE, stdin);
68                 if (ret >= 0) {
69                         pcm.insert(pcm.end(), buf, buf + ret);
70                 }
71         }       
72
73 #if 0
74         for (int i = 0; i < LEN; ++i) {
75                 in[i] += rand() % 10000;
76         }
77 #endif
78
79 #if 0
80         for (int i = 0; i < LEN; ++i) {
81                 printf("%d\n", in[i]);
82         }
83 #endif
84
85         std::vector<pulse> pulses;  // in seconds
86
87         // Find the flanks.
88         int last_bit = -1;
89         double last_downflank = -1;
90         for (unsigned i = 0; i < pcm.size(); ++i) {
91                 int bit = (pcm[i] > 0) ? 1 : 0;
92                 if (bit == 0 && last_bit == 1) {
93                         // Check if we ever go up above HYSTERESIS_LIMIT before we dip down again.
94                         bool true_pulse = false;
95                         unsigned j;
96                         int min_level_after = 32767;
97                         for (j = i; j < pcm.size(); ++j) {
98                                 min_level_after = std::min<int>(min_level_after, pcm[j]);
99                                 if (pcm[j] > 0) break;
100                                 if (pcm[j] < -HYSTERESIS_LIMIT) {
101                                         true_pulse = true;
102                                         break;
103                                 }
104                         }
105
106                         if (!true_pulse) {
107 #if 0
108                                 fprintf(stderr, "Ignored down-flank at %.6f seconds due to hysteresis (%d < %d).\n",
109                                         double(i) / SAMPLE_RATE, -min_level_after, HYSTERESIS_LIMIT);
110 #endif
111                                 i = j;
112                                 continue;
113                         } 
114
115                         // down-flank!
116                         double t = find_zerocrossing(pcm, i - 1) * (1.0 / SAMPLE_RATE);
117                         if (last_downflank > 0) {
118                                 pulse p;
119                                 p.time = t;
120                                 p.len = t - last_downflank;
121                                 pulses.push_back(p);
122                         }
123                         last_downflank = t;
124                 }
125                 last_bit = bit;
126         }
127
128         // Calibrate on the first ~25k pulses (skip a few, just to be sure).
129         double calibration_factor = 1.0f;
130         if (pulses.size() < SYNC_PULSE_END) {
131                 fprintf(stderr, "Too few pulses, not calibrating!\n");
132         } else {
133                 double sum = 0.0;
134                 for (int i = SYNC_PULSE_START; i < SYNC_PULSE_END; ++i) {
135                         sum += pulses[i].len;
136                 }
137                 double mean_length = C64_FREQUENCY * sum / (SYNC_PULSE_END - SYNC_PULSE_START);
138                 calibration_factor = SYNC_PULSE_LENGTH / mean_length;
139                 fprintf(stderr, "Calibrated sync pulse length: %.2f -> %.2f (change %+.2f%%)\n",
140                         mean_length, SYNC_PULSE_LENGTH, 100.0 * (calibration_factor - 1.0));
141
142                 // Check for pulses outside +/- 10% (sign of misdetection).
143                 for (int i = SYNC_PULSE_START; i < SYNC_PULSE_END; ++i) {
144                         double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
145                         if (cycles < SYNC_PULSE_LENGTH / SYNC_TEST_TOLERANCE || cycles > SYNC_PULSE_LENGTH * SYNC_TEST_TOLERANCE) {
146                                 fprintf(stderr, "Sync cycle with downflank at %.6f was detected at %.0f cycles; misdetect?\n",
147                                         pulses[i].time, cycles);
148                         }
149                 }
150
151                 // Compute the standard deviation (to check for uneven speeds).
152                 double sum2 = 0.0;
153                 for (int i = SYNC_PULSE_START; i < SYNC_PULSE_END; ++i) {
154                         double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
155                         sum2 += (cycles - SYNC_PULSE_LENGTH) * (cycles - SYNC_PULSE_LENGTH);
156                 }
157                 double stddev = sqrt(sum2 / (SYNC_PULSE_END - SYNC_PULSE_START - 1));
158                 fprintf(stderr, "Sync pulse length standard deviation: %.2f cycles\n",
159                         stddev);
160         }
161
162         FILE *fp = fopen("cycles.plot", "w");
163         std::vector<char> tap_data;
164         for (unsigned i = 0; i < pulses.size(); ++i) {
165                 double cycles = pulses[i].len * calibration_factor * C64_FREQUENCY;
166                 fprintf(fp, "%f %f\n", pulses[i].time, cycles);
167                 int len = lrintf(cycles / TAP_RESOLUTION);
168                 if (i > SYNC_PULSE_END && (cycles < 100 || cycles > 800)) {
169                         fprintf(stderr, "Cycle with downflank at %.6f was detected at %.0f cycles; misdetect?\n",
170                                         pulses[i].time, cycles);
171                 }
172                 if (len <= 255) {
173                         tap_data.push_back(len);
174                 } else {
175                         int overflow_len = lrintf(cycles);
176                         tap_data.push_back(0);
177                         tap_data.push_back(overflow_len & 0xff);
178                         tap_data.push_back((overflow_len >> 8) & 0xff);
179                         tap_data.push_back(overflow_len >> 16);
180                 }
181         }
182         fclose(fp);
183
184         tap_header hdr;
185         memcpy(hdr.identifier, "C64-TAPE-RAW", 12);
186         hdr.version = 1;
187         hdr.reserved[0] = hdr.reserved[1] = hdr.reserved[2] = 0;
188         hdr.data_len = tap_data.size();
189
190         fwrite(&hdr, sizeof(hdr), 1, stdout);
191         fwrite(tap_data.data(), tap_data.size(), 1, stdout);
192 }