2 // dsp.cpp: Digital Signal Processing module
4 // Part of the WOZ Maker project
6 // (C) 2018 Underground Software
19 // Really, there shouldn't be more than 5 waveforms per track
21 double initSyncTime[6];
23 static const double SIGMA_THRESHOLD = 2.5;
26 double StdDeviation(double * a, uint32_t num, double * pMean/*= NULL*/)
28 double sigma = 0, mean = 0;
30 for(uint32_t i=0; i<num; i++)
33 mean = mean / (double)num;
38 for(uint32_t i=0; i<num; i++)
39 sigma += (a[i] - mean) * (a[i] - mean);
41 return sqrt(sigma / (double)num);
45 double StdDeviationWave(uint32_t * sync, uint32_t num, double * pMean/*= NULL*/)
47 double sigma = 0, mean = 0;
50 for(uint32_t i=0; i<num; i++)
52 a[i] = (double)Global::wave[sNum[i]][sync[i]];
56 mean = mean / (double)num;
58 // Capture the mean, if a pointer was passed in
62 for(uint32_t i=0; i<num; i++)
63 sigma += (a[i] - mean) * (a[i] - mean);
65 return sqrt(sigma / (double)num);
69 double Largest(double * a, uint32_t num)
71 double largest = a[0];
73 for(uint32_t i=1; i<num; i++)
83 uint32_t Largest(uint32_t * a, uint32_t num)
85 uint32_t largest = a[0];
87 for(uint32_t i=1; i<num; i++)
97 uint32_t LargestIndex(double * a, uint32_t num)
99 double largest = a[0];
102 for(uint32_t i=1; i<num; i++)
115 uint32_t LargestIndex(uint32_t * a, uint32_t num)
117 uint32_t largest = a[0];
120 for(uint32_t i=1; i<num; i++)
133 double Smallest(double * a, uint32_t num)
135 double smallest = a[0];
137 for(uint32_t i=1; i<num; i++)
147 uint32_t Smallest(uint32_t * a, uint32_t num)
149 uint32_t smallest = a[0];
151 for(uint32_t i=1; i<num; i++)
161 uint32_t SmallestIndex(double * a, uint32_t num)
163 double smallest = a[0];
166 for(uint32_t i=1; i<num; i++)
179 uint32_t SmallestIndex(uint32_t * a, uint32_t num)
181 uint32_t smallest = a[0];
184 for(uint32_t i=1; i<num; i++)
198 // Compare two arrays to see if they are equal
200 bool Equal(uint32_t * a1, uint32_t * a2, uint32_t num)
202 for(uint32_t i=0; i<num; i++)
212 uint32_t IndexForTime(uint32_t trackNum, uint32_t time)
214 uint32_t curTime = 0, pos = 0;
216 while (curTime < time)
217 curTime += Global::synWave[trackNum][pos++];
224 // This is where we attempt to use the stream data to collapse them down to a single track. Taken as wholes, each individual track will drift out of sync with respect to the others, and fairly quickly as the drive mechanism does not read consistently at a set RPM as you might expect, but instead it fluctuates all over the place--which throws a big ol' monkey wrench into the works. (Yes, the overall *average* stays around a set RPM, but locally--nope.)
226 // The big problem is that the disk spins via a motor which isn't perfect in its rotation that is connected to the hub that spins the disk via a belt, which further adds more uncertainty to the speed of the disk as it rotates. The end result is that there is no consistency vis-a-vis the speed of the disk, and so the results of one read to the next can and will vary wildly. It's a miracle that these things are readable at all. O_o
228 // So how do we cope with it? I'm sure there's a well known paper or set of papers which addresses noise rejection from multiple signals read from a single source, but I have yet to find them. So what we have below is a set of ad-hoc assumptions and heuristics that attempt to deal with this problem.
230 // The basic assumption is that from any given pulse to the next one, the drift will tend to be small and thus as we move from pulse to pulse, we can stay more or less in sync by resetting the base from which we calculate the next time interval to the previous step. This works remarkably well until we run into a long period of zeroes, which seems to throw off the Applesauce hardware which read these streams in the first place; or when some spurious one pulses pop up in one stream that doesn't in one or more of the others.
233 void FindInitialSyncForStreams2(uint32_t * sync, uint32_t num)
235 bool lookAt[6] = { true, true, false, false, false, false };
237 // while (FindSyncBetweenFirstTwo(sync) == false)
238 while (FindSyncBetweenFirstTwo(sync) < 1000)
240 for(uint32_t i=0; i<num; i++)
244 if (sync[i] >= Global::waveLen[sNum[i]])
246 printf("Could not find initial sync!!\n");
252 uint32_t lookahead = LookaheadWave(sync, 2);
255 printf("Found sync for streams 1 & 2 at %d and %d (la=%d)\n", sync[0], sync[1], lookahead);
257 // Add in the others, one at a time and try to synchronize them
258 uint32_t bestSync[6];
260 for(uint32_t i=2; i<num; i++)
265 for(uint32_t j=0; j<200; j++)
267 uint32_t syncSave = sync[i];
269 // And check again for a match
270 for(uint32_t k=0; k<200; k++)
272 uint32_t la = LookaheadWave(sync, i + 1);
277 memcpy(bestSync, sync, sizeof(bestSync));
285 // Kick synchronized streams forward to check the rest...
286 for(uint32_t k=0; k<i; k++)
290 memcpy(sync, bestSync, sizeof(bestSync));
293 for(uint32_t i=0; i<num; i++)
294 printf("Sync for stream %d: %d (%s)\n", i, sync[i], (lookAt[i] ? "yes" : "no"));
295 printf("Lookahead: %d\n", LookaheadWave(sync, num));
300 uint32_t Synthesize2(uint32_t * sync, uint32_t num, uint32_t * synthWave, float * swAmplitude)
303 uint32_t syncPt[7], syncSave[7];
304 double a[7], b[7], currentRow[7];
307 memcpy(syncPt, sync, sizeof(syncPt));
311 for(uint32_t i=0; i<num; i++)
313 if (syncPt[i] >= Global::streamLen[i])
317 // Save these *before* we step forward
318 memcpy(syncSave, syncPt, sizeof(syncSave));
319 // StepForward(syncPt, a, num);
320 for(uint32_t i=0; i<num; i++)
322 a[i] = (double)Global::wave[sNum[i]][syncPt[i]];
327 double sigma = StdDeviation(a, num, &mean);
329 if (sigma <= SIGMA_THRESHOLD)
331 synthWave[swLen] = (uint32_t)(mean + 0.5);
332 swAmplitude[swLen] = 1.0f;
338 // The basic idea is to find the extent of this anomalous region and then use what's there as an amplitude by summing the signals where they correspond.
339 uint32_t syncSteps = 0;
340 uint32_t lookahead = LookaheadWave(syncPt, num);
342 // Basically we move forward thru the area until either 1) the sigma is OK, or 2) we've pulled everything as far as it can go and thus this becomes our new local start. 3) will go in at some point which is, do a lookahead and see what the sigma looks like; if it's good and the current one is stinko, (not *too* stinko) then we've found a new local start point, otherwise keep going.
344 Note that 3) works OK to get out of sticky situations, but sometimes the signals are fairly noisy and you will end up in situations where the lookahead is zero and the next sigma isn't great (& the current isn't any good either), but using that as the new local start is OK and will permit walking forward/backward through the stream. Some of the time you will find the average of the current sigma and the following will be below the threshold with a good lookahead value (10, 100, 1000), this is also a good time to set the new local start.
348 for(uint32_t i=0; i<num; i++)
349 printf("[%.2lf]", a[i]);
350 //printf(" => %.2lf {drift=%.2lf, la=%d}\n", sigma, driftSigma, lookahead);
351 printf(" => %.2lf {la=%d}\n", sigma, lookahead);
353 for(uint32_t j=0; j<num; j++) { printf("{%d}%s", syncPt[j], (j < num - 1 ? "-" : "")); }
356 bool keepWorking = true;
362 // If we've been in this loop too long, something's very wrong
366 // Pull the lesser streams forward to catch up to the largest
367 uint32_t largestIdx = LargestIndex(a, num);
370 for(uint32_t i=0; i<num; i++)
381 // if ((a[largestIdx] - a[i]) > 7.0)
382 if ((a[largestIdx] - a[i]) > 9.0)
385 a[i] += (double)Global::wave[sNum[i]][syncPt[i]];
387 printf(" %6d ", Global::wave[sNum[i]][syncPt[i]]);
395 // See if we've run out of stream to synthesize...
396 if (syncPt[i] >= Global::waveLen[sNum[i]])
403 for(uint32_t i=0; i<num; i++)
404 printf("<%6.2lf>", a[i]);
406 // Calculate the standard deviation at this new point
407 sigma = StdDeviation(a, num);//, &mean);
409 // Calculate the following sigma & lookahead too
410 for(uint32_t i=0; i<num; i++)
411 b[i] = (double)Global::wave[sNum[i]][syncPt[i]];
413 double nextSigma = StdDeviation(b, num);
414 lookahead = LookaheadWave(syncPt, num);
416 printf(" ==> %.2lf%s", sigma, (sigma <= SIGMA_THRESHOLD ? " (OK)" : (!pulled ? "(STEP)" : "")));
417 printf(" (lookahead=%.2lf, %d)\n", nextSigma, lookahead);
419 for(uint32_t j=0; j<num; j++) { printf("{%d}%s", syncPt[j], (j < num - 1 ? "-" : "")); }
422 // A handful of heuristics. If the sigma is OK, we're done
426 // If we couldn't pull any streams forward for some reason, punt
430 // If the average of the current & next sigmas is OK, go ahead
431 if (((sigma + nextSigma) / 2.0) < 2.5)
434 // If the next sigma is OK and the lookahead too, move on
435 if ((nextSigma < 2.5) && (lookahead > 10))
440 // Now, we've re-established sync, now fill the appropriate spots
443 So the algorithm to do this proppaly is like so:
445 Step 1: Find the smallest step forward through the section; set amplitude to 0.2.
446 Step 2: Loop thru the other streams and see if they correspond. If so, advance their stream start & add 0.2 to the amplitude for each one seen.
447 Step 3: Add the mean value of the corresponding values to the synthesized waveform + wave amplitude
448 Step 4: If all of the stream starts are at the end of the section, we are done.
449 [maybe. how do you handle the case where you have streams that look like so:
457 .......................................|
458 .............................|.........|
459 .............................|.........|
460 ..................................|....|
461 .......................................|
463 this will fail to add in the correct amount for the times that had no part in the anomaly. so to fix this, maybe have a baseline value that bumps upward to show where the "cursor" (so to speak) is, so at the end we can subtract the mean that we got above from this baseline to get the remaining time to synthesize, like so:
465 synthWave[swLen] = (uint32_t)((mean - baseline) + 0.5)
466 swAmplitude[swLen] = amplitude;
468 no need for that. the times that step directly to the end can safely be ignored.
471 150 -> 150 (ignores 175)
472 150 goes in, ampl. .4, we now have:
480 only one is left that we can grab:
484 last pulse was at 150 (40%), so we need to subtract that out
488 so maybe we'd do like so:
489 -1) create a list containing the start values for the start of the anomalous section and replace them as the counters advance.
490 0) are streams at the end? if so, stuff in the final value to the synthesized wave & stop, otherwise:
491 1) find the smallest in the current row and use that as a basis.
492 2) find streams that correlate, if any
493 3) advance the streams that correlate & subtract out the minimum from the other rows
502 smallest is 150: pick 150 -> picks up 150
503 advance the 150 rows & subtract 150 from the rest:
511 pick smallest: 25 -> picks up nothing
512 advance the 25 & subtract it from the rest:
520 We've reached the end of the section, add the final amount
522 --------------------------------------------------------------------------------
523 syncSave before: [2->4][2->6][2->4][2->6][2->6]
524 syncSave after: [3->4][3->6][3->4][3->6][3->6]
526 currentRow: [120.89][49.97][46.98][49.98][48.98] => 28.78
528 2: [3->4][3->6][3->4][3->6][3->6]
530 a[0]=46.98: [120.89] ==> 36.95
531 [49.97]{17.00} ==> 1.50
532 <106.00>[49.97][49.98]{28.00} ==> 1.41
533 [49.97][49.98][48.98]{26.00} ==> 1.22
535 3: [3->4][4->6][4->4][4->6][4->6]
537 a[0]=17.00: [73.91] ==> 28.45
538 <56.00>[28.00]{45.00} ==> 5.50
539 [28.00][26.00]{47.00} ==> 4.78
541 4: [3->4][5->6][4->4][5->6][5->6]
543 a[0]=45.00: [56.91] ==> 5.95
544 [56.00]{32.00} ==> 5.50
545 <32.00>[56.00][47.00]{32.00} ==> 4.78
547 5: [3->4][6->6][4->4][6->6][6->6]
551 6: [4->4][6->6][4->4][6->6][6->6]
553 2->4 2->6 2->4 2->6 2->6
554 [120.89][ 49.97][ 46.98][ 49.98][ 48.98] => 28.78
555 ( 73.91) 17.00 106.00 28.00 26.00
556 ( 56.91) 56.00 45.00 47.00
557 32.00 32.00 32.00 32.00
561 /*printf("syncSave before: ");
562 for(uint32_t i=0; i<num; i++)
563 printf("[%d->%d]", syncSave[i], syncPt[i]);
566 for(uint32_t i=0; i<num; i++)
568 printf("Sync %d: %d -> %d (", sNum[i], syncSave[i], syncPt[i]);
569 for(uint32_t j=syncSave[i]; j<=syncPt[i]; j++)
570 printf("%d%s", Global::wave[sNum[i]][j], (j == syncPt[i] ? "" : ", "));
574 // Collect the current row
575 for(uint32_t i=0; i<num; i++)
577 currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
581 /*printf("syncSave after: ");
582 for(uint32_t i=0; i<num; i++)
583 printf("[%d->%d]", syncSave[i], syncPt[i]);
585 printf("currentRow: ");
586 for(uint32_t i=0; i<num; i++)
587 printf("[%.2lf]", currentRow[i]);
588 printf(" => %.2lf\n", StdDeviation(currentRow, num));//*/
590 //printf("%d: ", swLen);
592 float ampStep = 1.0f / (float)num;
596 if (Equal(syncSave, syncPt, num))
599 float amplitude = ampStep;
600 uint32_t smallestIdx = SmallestIndex(currentRow, num);
602 a[0] = currentRow[smallestIdx];
603 //printf("a[0]=%.2lf: ", a[0]);
605 for(uint32_t i=0; i<num; i++)
607 // Make sure we haven't gone beyond the end of any of the streams we're looking at
608 if (syncSave[i] >= Global::streamLen[i])
611 if (i == smallestIdx)
613 // We always pick up the smallest in the current set, since that's our basis
614 currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
616 //printf("<%.2lf>", currentRow[i]);
620 // Save the current row & subtract out the smallest as a basis
621 double crSave = currentRow[i];
622 currentRow[i] -= a[0];
624 // Don't bother with stuff that's out of range (@ current sync point!) :-P
625 if (syncSave[i] == syncPt[i])
628 // Try to find correlates
631 sigma = StdDeviation(a, cnt);
632 //for(uint32_t j=1; j<cnt; j++)
633 // printf("[%.2lf]", a[j]);
636 if (sigma <= 5.5)//SIGMA_THRESHOLD)
638 // If the sigma is in range for this item, pick it up
639 currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
642 //printf("{%.2lf}", currentRow[i]);
643 amplitude += ampStep;
647 // Otherwise reject the last item; it's too far out
650 //printf(" ==> %.2lf\n\t", sigma);
653 // Get the average of the items picked up and add it to the wave
654 sigma = StdDeviation(a, cnt, &mean);
655 synthWave[swLen] = (uint32_t)(mean + 0.5);
656 swAmplitude[swLen] = amplitude;
660 sigma = StdDeviation(currentRow, num, &mean);
661 synthWave[swLen] = mean;
662 swAmplitude[swLen] = 1.0f;
672 bool ResyncWave(uint32_t * sync, uint32_t num, int8_t dir/*= 1*/)
676 for(uint32_t i=0; i<num; i++)
678 a[i] = (double)Global::wave[sNum[i]][sync[i]];
683 for(uint32_t i=0; i<num; i++)
684 printf("<%4d>", (uint32_t)a[i]);
686 printf(" {la: %d}\n", LookaheadWave(sync, num, dir));
689 uint32_t syncCount = 0;
690 bool keepGoing = true;
696 printf("Could not synchronize streams...\n");
700 for(uint32_t i=0; i<num; i++)
702 if (((dir == 1) && (sync[i] >= Global::waveLen[sNum[i]]))
703 || ((dir == -1) && (sync[i] == 0)))
705 printf("End of stream...\n");
710 // Deal with the problem. The heuristic is find the largest value in the set, then advance the other streams by one step until they are "close" to the largest. If the std. deviation is too large, try again; eventually it should be within good limits.
711 uint32_t largestIdx = LargestIndex(a, num);
712 bool stepped = false;
714 for(uint32_t i=0; i<num; i++)
726 // If the difference between the highest and this is over 7, add the next time interval and let the sync counter, uh, slip :-)
727 // if ((a[largestIdx] - a[i]) > 7.0)
728 if ((a[largestIdx] - a[i]) > 9.0)
731 double ttn = (double)Global::wave[sNum[i]][sync[i]];
735 printf(" %4d ", (uint32_t)ttn);
744 double sigma = StdDeviation(a, num);
745 double nextSigma = StdDeviationWave(sync, num);
746 uint32_t lookahead = LookaheadWave(sync, num, dir);
750 for(uint32_t i=0; i<num; i++)
751 printf("<%4d>", (uint32_t)a[i]);
753 printf(" ==> %.2lf ", sigma);
754 printf("(next: %.2lf, la: %d)\n", nextSigma, lookahead);
758 // A handful of heuristics... 1st: If sigma is good, we're done
762 // Couldn't pull, so punt for now
766 // If the next sigma + current / 2 is OK, move on
767 if (((sigma + nextSigma) / 2.0) < 2.50)
770 // If the next sigma is good, and the lookahead is good, move on
771 // [Slight tweak: only do this if the current sigma isn't terrible]
772 if ((nextSigma < 2.50) && (lookahead > 10) && (sigma < 20.0))
781 uint32_t LookaheadWave(uint32_t * passedInSync, uint32_t num, int8_t dir/*= 1*/, bool * perfect/*= NULL*/)
785 memcpy(sync, passedInSync, sizeof(sync));
792 if (StdDeviationWave(sync, num) > 2.50)
800 for(uint32_t i=0; i<num; i++)
804 if (((dir == 1) && (sync[i] >= Global::waveLen[sNum[i]]))
805 || ((dir == -1) && (sync[i] == 0)))
816 uint32_t FindSyncBetweenFirstTwo(uint32_t * sync)
818 uint32_t bestSync[2], bestLA = 0, syncSave = sync[0];
821 // Move the first stream first...
822 for(uint32_t i=0; i<200; i++)
824 uint32_t la = LookaheadWave(sync, 2, 1, &perfect);
826 // See if the sync we've found is the best so far...
830 memcpy(bestSync, sync, sizeof(bestSync));
839 // Otherwise, reset the first stream and move the 2nd
842 for(uint32_t i=0; i<200; i++)
844 // We can bump this first since we already checked the [0][0] case above
846 uint32_t la = LookaheadWave(sync, 2, 1, &perfect);
848 // See if the sync we've found is the best so far...
852 memcpy(bestSync, sync, sizeof(bestSync));
859 // Set sync to the best we found & return the best lookahead count
860 memcpy(sync, bestSync, sizeof(bestSync));
867 // Walk backwards through the streams to try to find the starting point
869 bool AttemptToFindStart(uint32_t * sync, uint32_t num)
872 for(uint32_t i=0; i<num; i++)
878 bool keepGoing = true;
882 double sigma = StdDeviationWave(sync, num);
886 bool result = ResyncWave(sync, num, -1);
893 for(uint32_t i=0; i<num; i++)
897 for(uint32_t i=0; i<num; i++)
907 for(uint32_t i=0; i<num; i++)
908 printf("Sync for stream %d: %d\n", i, sync[i]);
909 printf("Lookahead: %d\n", LookaheadWave(sync, num));
916 bool InitialSync(uint32_t * sync, uint32_t num)
918 // Let's try a different approach to this...
919 bool lookAt[6] = { true, true, false, false, false, false };
920 uint32_t syncSave[2];
922 // while (FindSyncBetweenFirstTwo(sync) == false)
923 // while (FindSyncBetweenFirstTwo(sync) < 1000)
926 memcpy(syncSave, sync, sizeof(syncSave));
927 uint32_t bestSBFT = FindSyncBetweenFirstTwo(sync);
929 if (bestSBFT >= 1000)
932 memcpy(sync, syncSave, sizeof(syncSave));
934 for(uint32_t i=0; i<num; i++)
938 if (sync[i] >= Global::waveLen[sNum[i]])
940 printf("Could not find initial sync!!\n");
946 uint32_t lookahead = LookaheadWave(sync, 2);
949 printf("Found sync for streams 1 & 2 at %d and %d (la=%d)\n", sync[0], sync[1], lookahead);
951 // Add in the others, one at a time and try to synchronize them
952 uint32_t syncBest[6];
955 for(uint32_t i=2; i<num; i++)
960 for(uint32_t j=0; j<200; j++)
962 uint32_t syncSave = sync[i];
964 for(uint32_t k=0; k<200; k++)
966 uint32_t la = LookaheadWave(sync, i + 1, 1, &perfect);
971 memcpy(syncBest, sync, sizeof(syncBest));
974 if (perfect || ((bestLA + j) == lookahead))
985 // Kick all the other streams forward to check the rest...
986 for(uint32_t k=0; k<i; k++)
990 memcpy(sync, syncBest, sizeof(syncBest));
993 for(uint32_t i=0; i<num; i++)
994 printf("Sync for stream %d: %d (%s)\n", i, sync[i], (lookAt[i] ? "yes" : "no"));
995 printf("Lookahead: %d\n", LookaheadWave(sync, num));
1001 uint32_t LoopLookahead(uint32_t trackNum, uint32_t start, uint32_t loopPoint, bool * test/*= NULL*/)
1007 for(uint32_t i=loopPoint, pos=start; i<Global::synWaveLen[trackNum]; i++, pos++)
1009 a[0] = Global::synWave[trackNum][pos];
1010 a[1] = Global::synWave[trackNum][i];
1011 double sigma = StdDeviation(a, 2);
1030 //#define DEBUGSYNTH
1031 void SynthesizeTrack(uint32_t trackNum)
1034 uint32_t sync[6] = { 0 };
1036 for(uint32_t i=0; i<Global::numStreams; i++)
1038 // Skip streams that aren't on the current track & BITS captures
1039 if ((Global::stream[i]->location != trackNum)
1040 || (Global::stream[i]->captureType == 2))
1044 printf("sNum[%d] = %d\n", num, i);
1046 // Save the stream # so we can find it later
1050 // We should just bypass the synthesis part of this, not the loop point finding... :-P !!! FIX !!!
1053 printf("%.2f: Single stream, length: %d\n", (float) trackNum / 4.0f, Uint32LE(Global::stream[sNum[0]]->dataLength));
1055 memcpy(Global::synWave[trackNum], Global::wave[sNum[0]], Global::waveLen[sNum[0]] * 4);
1056 Global::synWaveLen[trackNum] = Global::waveLen[sNum[0]];
1060 This is failing for track 9.00: it finds an initial sync @ 202, 200, la=1257 but real sync point is 0,2,2 (it says 5,0,0, which is wrong)
1064 InitialSync(sync, num);
1065 AttemptToFindStart(sync, num);
1067 uint32_t count = 0, bloops = 0;
1068 bool keepGoing = true;
1069 Global::synWaveLen[trackNum] = 0;
1074 double sigma = StdDeviationWave(sync, num, &mean);
1079 Global::synWave[trackNum][Global::synWaveLen[trackNum]] = (uint32_t)(mean + 0.5);
1080 Global::synWaveLen[trackNum]++;
1082 for(uint32_t i=0; i<num; i++)
1088 uint32_t startOfSynth = Global::synWaveLen[trackNum];
1093 printf("Bloop seen at pos ");
1094 for(uint32_t k=0; k<num; k++)
1095 printf("%d%s", sync[k], (k == num - 1 ? "" : ", "));
1099 bloops++; // Count the # of times we've had to punt
1100 uint32_t syncSave[6] = { 0 };
1101 memcpy(syncSave, sync, sizeof(syncSave));
1102 /*bool result =*/ ResyncWave(sync, num);
1104 // Now go through the catty wampus crap and zero out the parts that aren't all ones across the board
1105 double timeToOnes = 0;
1106 double a[6], currentRow[6];
1108 // Collect the current row
1109 for(uint32_t i=0; i<num; i++)
1111 currentRow[i] = (double)Global::wave[sNum[i]][syncSave[i]];
1115 for(uint32_t i=0; i<num; i++)
1116 printf("syncSave[%d]=%d (sync: %d)\n", i, syncSave[i], sync[i]);
1121 if (Equal(syncSave, sync, num))
1124 printf("syncSave: ");
1125 for(uint32_t i=0; i<num; i++)
1126 printf("%d%s", syncSave[i], (i == num - 1 ? "" : ", "));
1127 printf("\n sync: ");
1128 for(uint32_t i=0; i<num; i++)
1129 printf("%d%s", sync[i], (i == num - 1 ? "" : ", "));
1133 uint32_t amplitude = 1;
1134 uint32_t smallestIdx = SmallestIndex(currentRow, num);
1136 a[0] = currentRow[smallestIdx];
1138 printf("a[0]=%.2lf:\n", a[0]);
1141 for(uint32_t i=0; i<num; i++)
1143 // Make sure we haven't gone beyond the end of any of the streams we're looking at
1144 if (syncSave[i] >= Global::waveLen[sNum[i]])
1146 printf("Ran past end of stream %d (stream %d): syncSave=%d, waveLen=%d\n", i, sNum[i], syncSave[i], Global::waveLen[sNum[i]]);
1150 if (i == smallestIdx)
1152 // Maybe we should chuck this above the if()?
1153 // if (syncSave[i] == sync[i])
1156 // We always pick up the smallest in the current set, since that's our basis
1157 //the reason this works is because this simulates the "add next and subtract smallest", since the subtraction of the smallest would always be zero. :-P
1158 currentRow[i] = Global::wave[sNum[i]][syncSave[i]];
1159 // currentRow[i] += Global::wave[sNum[i]][syncSave[i]];
1162 printf("\t<%.2lf> (i=%d)\n", currentRow[i], i);
1167 // Save the current row & subtract out the smallest as a basis
1168 double crSave = currentRow[i];
1169 currentRow[i] -= a[0];
1171 // Don't bother with stuff that's out of range (@ current sync point!) :-P
1172 if (syncSave[i] == sync[i])
1175 // Try to find correlates
1178 sigma = StdDeviation(a, cnt);
1181 for(uint32_t j=1; j<cnt; j++)
1182 printf("[%.2lf]", a[j]);
1186 if (sigma <= 5.5)//SIGMA_THRESHOLD)
1188 // If the sigma is in range for this item, pick it up
1189 currentRow[i] += (double)Global::wave[sNum[i]][syncSave[i]];
1193 printf("{%.2lf}", currentRow[i]);
1199 // Otherwise reject the last item; it's too far out
1203 printf(" ==> %.2lf\n", sigma);
1207 // Get the average of the items picked up and add it to the wave
1208 sigma = StdDeviation(a, cnt, &mean);
1211 // Did we find the line of ones yet? If so, chuck it in the wave
1212 if (amplitude == num)
1214 Global::synWave[trackNum][Global::synWaveLen[trackNum]] = (uint32_t)(timeToOnes + 0.5);
1215 Global::synWaveLen[trackNum]++;
1220 for(uint32_t j=0; j<num; j++)
1221 printf("_%.2lf_", currentRow[j]);
1226 // Finally, chuck the final line of ones in
1227 sigma = StdDeviation(currentRow, num, &mean);
1229 Global::synWave[trackNum][Global::synWaveLen[trackNum]] = (uint32_t)(timeToOnes + 0.5);
1230 Global::synWaveLen[trackNum]++;
1232 printf("Blip @ %d: ", startOfSynth);
1233 for(uint32_t w=startOfSynth; w<Global::synWaveLen[trackNum]; w++)
1234 printf("%d%s", Global::synWave[trackNum][w], (w == Global::synWaveLen[trackNum] - 1 ? "" : ", "));
1239 for(uint32_t i=0; i<num; i++)
1241 if (sync[i] >= Global::waveLen[sNum[i]])
1250 printf("\n%.2f: Synthesis done. Bit count: %d, bloop count: %d\n", (float)trackNum / 4.0f, count, bloops);
1253 printf("%.2f: Synthesized wave length: %d\n", (float)trackNum / 4.0f , Global::synWaveLen[trackNum]);
1255 // Now, find the loop point and make BITS from it
1257 Now, we take the track and make an ouroborus out of it by finding the sync point. We can also possibly use it to get rid of more false ones like we did above with the resynchronization code. So how do we do this?
1259 We could use an approach similar to the original sync between 1 & 2, but there's no guarantee that the lookahead will be good right away. So we need to do a combination where we count the places where it coincides well and plow through spots where it doesn't (and possibly mask any ones in those areas if it's a good match).
1261 We can start with the estimated loop point that the Applesauce hardware found, by walking it a bit before that point we can make sure we don't miss it.
1263 uint32_t loopStart, loopEnd;
1264 bool perfect = FindLoopPoint(trackNum, loopStart, loopEnd);
1266 // We use the loop point from stream #1 because we're lazy
1267 uint32_t estLoopTime = Uint32LE(Global::stream[sNum[0]]->estLoopPoint);
1268 uint32_t curTime = 0;
1269 uint32_t pos = 0, startOfLoop = 0;
1271 uint32_t maxLookahead = 0;
1272 uint32_t maxPos = 0;
1273 bool /*full, full2,*/ full0, full1;
1274 uint32_t blips, worstBlips = 0, longestStretch, stretch, maxLongestStretch = 0, worstBlipsPos = 0, mlsPos = 0;
1277 // Get the index to the correct time in "pos"
1278 while (curTime < estLoopTime)
1279 curTime += Global::synWave[trackNum][pos++];
1281 printf("%.2f: estLoopTime=%d, curTime=%d, pos=%d\n", (float)trackNum / 4.0f, estLoopTime, curTime, pos);
1282 uint32_t posSave = pos;
1284 // First, see if we have basically a perfect loop out of the box (typically, all $AAs)
1287 worstBlipsPos = pos;
1289 maxLongestStretch = LoopLookahead(trackNum, 0, pos, &perfect);
1295 // "pos" now holds the start of the search for our loop point.
1296 for(uint32_t i=0; i<200; i++)
1298 uint32_t lookahead = 0;
1302 longestStretch = stretch = 0;
1303 startOfLoop = 0;//shouldn't be necessary...
1305 for(j=(pos+i-100), pos0=0; j<Global::synWaveLen[trackNum]; j++, pos0++)
1307 uint32_t lla0 = LoopLookahead(trackNum, pos0 + 0, j, &full0);
1308 uint32_t lla1 = LoopLookahead(trackNum, pos0 + 1, j, &full1);
1313 if (stretch > longestStretch)
1314 longestStretch = stretch;
1320 else if (full1 && first) // only check 1st time (for now)
1334 a[0] = Global::synWave[trackNum][pos0++];
1335 a[1] = Global::synWave[trackNum][j++];
1338 double sigma = StdDeviation(a, 2);
1343 // printf("First blip: %.0lf, %.0lf (sigma=%.2lf) [pos=%d, la=%d]\n", a[0], a[1], sigma, i + pos - 100, pos0);
1345 if (trackNum == 137)
1347 uint32_t lla = LoopLookahead(trackNum, 0, j, &full);
1348 uint32_t llap1 = LoopLookahead(trackNum, 1, j, &full2);
1349 printf("%.2f: pos0=%d, j=%d, sigma=%lf, la=%d%s, la+1=%d%s", (float)trackNum / 4.0f, pos0, j, sigma, lla, (full ? " *equal" : ""), llap1, (full2 ? " *equal" : ""));
1354 // Hmm, not in sync. See if we can re-sync
1361 a[0] += Global::synWave[trackNum][pos0++];
1363 a[1] += Global::synWave[trackNum][j++];
1365 // Maybe put all the other heuristics in too?
1366 if (StdDeviation(a, 2) <= 2.50)
1369 // Now you've gone *too* far...
1370 if (j >= Global::synWaveLen[trackNum])
1377 if (trackNum == 137)
1378 printf(" blips=%d, la=%d\n", blips, lookahead);
1381 //maybe we should keep an array of the analyses so we can decide afterward which one is best...
1382 if (blips < worstBlips)
1385 worstBlipsPos = pos + i - 100;
1388 if (lookahead > maxLookahead)
1390 maxLookahead = lookahead;
1391 maxPos = pos + i - 100;
1394 if (longestStretch > maxLongestStretch)
1396 maxLongestStretch = longestStretch;
1397 mlsPos = pos + i - 100;
1407 bitLen[trackNum] = RenderBitstream(bits[trackNum], &byteLen[trackNum], Global::synWave[trackNum], startOfLoop, maxPos);
1412 uint32_t estLoopTime = Uint32LE(Global::stream[sNum[0]]->estLoopPoint);
1413 uint32_t pos = IndexForTime(trackNum, estLoopTime);
1414 printf("Perfect loop point found at position %d (elp%+d, start=%d, end=%d)\n----------------------------------------------------------------------\n", loopEnd + 1, loopEnd + 1 - pos, loopStart, loopEnd);
1417 printf("%s loop point found at position %d (elp%+d, la=%d, startPos=%d, blips=%d, blipsPos=%d, longest stretch=%d, mlsPos=%d, sol=%d)\n----------------------------------------------------------------------\n", (blips == 0 ? "Perfect" : "Possible"), maxPos, maxPos - posSave, maxLookahead, pos, worstBlips, worstBlipsPos, maxLongestStretch, mlsPos, startOfLoop);
1421 36110 is the loop point. But it fails to find it for some reason. [Because it jumped on the first one it found instead of the longest...]
1424 /* for(uint32_t i=0; i<maxLookahead+5; i++)
1426 a[0] = Global::synWave[trackNum][i];
1427 a[1] = Global::synWave[trackNum][i + maxPos];
1428 double sigma = StdDeviation(a, 2);
1429 printf("[%.2lf]", sigma);
1437 // Attempt to find where the track loop is. Passes back what it thinks they
1438 // are in start and end.
1440 bool FindLoopPoint(uint32_t trackNum, uint32_t & start, uint32_t & end)
1443 Now, we take the track and make an ouroborus out of it by finding the sync point. We can also possibly use it to get rid of more false ones like we did with the resynchronization code. So how do we do this?
1445 We could use an approach similar to the original sync between 1 & 2, but there's no guarantee that the lookahead will be good right away. So we need to do a combination where we count the places where it coincides well and plow through spots where it doesn't (and possibly mask any ones in those areas if it's a good match).
1447 We can start with the estimated loop point that the Applesauce hardware found, by walking it a bit before that point we can make sure we don't miss it.
1449 // We use the loop point from stream #1 because we're lazy (better would prolly be to take the average of all of them)
1450 uint32_t estLoopTime = Uint32LE(Global::stream[sNum[0]]->estLoopPoint);
1451 uint32_t pos = IndexForTime(trackNum, estLoopTime);
1452 uint32_t posSave = pos;
1454 bool full0, full1, perfect;
1456 printf("%.2f: estLoopTime=%d, pos=%d\n", (float)trackNum / 4.0f, estLoopTime, pos);
1458 // First, see if we have basically a perfect loop out of the box
1459 // (typically, a track of all $AAs)
1460 LoopLookahead(trackNum, 0, pos, &perfect);
1469 uint32_t bestMatchCount = 0, bestBlips = 1e9, longestStretch = 0;
1470 uint32_t bestMatchCountPos = 0, bestBlipsPos = 0, longestStretchPos = 0;
1472 // "pos" now holds the start of the search for our loop point.
1473 // Step from -99 to +99 around "pos" to find the loop point. If the input wave is very noisy, this might fail!
1474 for(int32_t i=-99; i<100; i++)
1476 uint32_t matchCount = 0;
1482 for(j=pos+i, pos0=0; j<Global::synWaveLen[trackNum]; j++, pos0++)
1484 uint32_t lla0 = LoopLookahead(trackNum, pos0 + 0, j, &full0);
1486 We do this because it may be the case that the start of the wave has an extraneous one that doesn't exist near the loop point; this way, we can still get a perfect loop point by skipping that bad one at the start.
1488 uint32_t lla1 = LoopLookahead(trackNum, pos0 + 1, j, &full1);
1490 // We count the # of successful matches, since there might be spots where it doesn't sync up perfectly--the combination of good matches plus least amount of desyncs ("blips") is what we're looking for
1493 if (lla0 > longestStretch)
1495 longestStretch = lla0;
1496 longestStretchPos = pos + i;
1501 // If we found perfect sync and it's the first time, return
1508 // Otherwise, just break out of the loop since it might not be the best
1511 else if (full1 && first) // only check 1st time (for now)
1513 // Since we found perfect sync, return
1526 a[0] = Global::synWave[trackNum][pos0++];
1527 a[1] = Global::synWave[trackNum][j++];
1530 double sigma = StdDeviation(a, 2);
1535 // printf("First blip: %.0lf, %.0lf (sigma=%.2lf) [pos=%d, la=%d]\n", a[0], a[1], sigma, i + pos - 100, pos0);
1537 if (trackNum == 137)
1539 uint32_t lla = LoopLookahead(trackNum, 0, j, &full);
1540 uint32_t llap1 = LoopLookahead(trackNum, 1, j, &full2);
1541 printf("%.2f: pos0=%d, j=%d, sigma=%lf, la=%d%s, la+1=%d%s", (float)trackNum / 4.0f, pos0, j, sigma, lla, (full ? " *equal" : ""), llap1, (full2 ? " *equal" : ""));
1546 // Hmm, not in sync. See if we can re-sync
1552 a[0] += Global::synWave[trackNum][pos0++];
1554 a[1] += Global::synWave[trackNum][j++];
1556 // Maybe put all the other heuristics in too? :-/
1557 if ((StdDeviation(a, 2) <= 2.50)
1558 || (j >= Global::synWaveLen[trackNum]))
1565 if (trackNum == 137)
1566 printf(" blips=%d, la=%d\n", blips, matchCount);
1569 //maybe we should keep an array of the analyses so we can decide afterward which one is best... probably not necessary...
1570 if (blips < bestBlips)
1573 bestBlipsPos = pos + i;
1576 if (matchCount > bestMatchCount)
1578 bestMatchCount = matchCount;
1579 bestMatchCountPos = pos + i;
1582 // Did we find a perfect lookahead? If so, we're done.
1587 // Now we're done trying to fit the loop, and we didn't find a perfect one, so set up the best we found:
1588 end = bestMatchCountPos - 1;
1590 uint32_t possibleCoverage = Global::synWaveLen[trackNum] - pos;
1591 printf("Possible loop point found at position %d (elp%+d, bestMatchCount=%d (out of %d, %.1f%%), startPos=%d, blips=%d, blipsPos=%d, longest stretch=%d, mlsPos=%d)\n----------------------------------------------------------------------\n", end + 1, end + 1 - posSave, bestMatchCount, possibleCoverage, ((float)bestMatchCount / (float)possibleCoverage) * 100.0f, pos, bestBlips, bestBlipsPos, longestStretch, longestStretchPos);
1598 Using a histogram, we can get an idea of the drift from the ones:
1616 -------------------------------------------------
1617 0000|0584|1792|9611|2627|5382|0482|0058|0013|0000
1618 28 29 30 31 32 33 34 35 36 37
1620 Peak of the ones occurs at 31, with a secondary at 33. Twos drift:
1631 -------------------------------------------------
1632 0000|0545|0277|5482|1223|5168|0252|0082|0102|0000
1633 57 58 59 60 61 62 63 64 65 66
1635 This looks like stretched out ones (they might all be spurious...):
1645 ---------------------------------------
1646 0000|0015|0007|0002|0011|0001|0002|0000
1647 43 44 45 46 47 48 49 50