From b1a1d1edc6208c186fd6f6eb18c0c3f91ad7e0c6 Mon Sep 17 00:00:00 2001 From: Karlis Goba Date: Tue, 13 Nov 2018 13:04:36 +0200 Subject: [PATCH] Tested sync candidate search --- .gitignore | 1 + decode_ft8.cpp | 173 +++++++++++++++++++++++++++++++------------------ 2 files changed, 111 insertions(+), 63 deletions(-) diff --git a/.gitignore b/.gitignore index a5515b3..ebfa2e2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ *.o gen_ft8 +decode_ft8 diff --git a/decode_ft8.cpp b/decode_ft8.cpp index dc19f0e..97de1ab 100644 --- a/decode_ft8.cpp +++ b/decode_ft8.cpp @@ -28,28 +28,71 @@ struct Candidate { int16_t score; uint16_t time_offset; uint16_t freq_offset; - uint8_t time_alt; - uint8_t freq_alt; + uint8_t time_sub; + uint8_t freq_sub; }; -void find_candidates(int num_blocks, int num_bins, const uint8_t * power, - int num_candidates, Candidate heap[num_candidates]) { +void heapify_down(Candidate * heap, int heap_size) { + // heapify from the root down + int current = 0; + while (true) { + int largest = current; + int left = 2 * current + 1; + int right = left + 1; + + if (left < heap_size && heap[left].score < heap[largest].score) { + largest = left; + } + if (right < heap_size && heap[right].score < heap[largest].score) { + largest = right; + } + if (largest == current) { + break; + } + + Candidate tmp = heap[largest]; + heap[largest] = heap[current]; + heap[current] = tmp; + current = largest; + } +} + + +void heapify_up(Candidate * heap, int heap_size) { + // heapify from the last node up + int current = heap_size - 1; + while (current > 0) { + int parent = (current - 1) / 2; + if (heap[current].score >= heap[parent].score) { + break; + } + + Candidate tmp = heap[parent]; + heap[parent] = heap[current]; + heap[current] = tmp; + current = parent; + } +} + + +// Find top N candidates in frequency and time according to their sync strength (looking at Costas symbols) +void find_sync(const uint8_t * power, int num_blocks, int num_bins, int num_candidates, Candidate * heap) { // Costas 7x7 tone pattern const uint8_t ICOS7[] = { 2,5,6,0,4,1,3 }; int heap_size = 0; for (int alt = 0; alt < 4; ++alt) { - for (int i = 0; i < num_blocks - NN; ++i) { - for (int j = 0; j < num_bins - 8; ++j) { + for (int time_offset = 0; time_offset < num_blocks - NN; ++time_offset) { + for (int freq_offset = 0; freq_offset < num_bins - 8; ++freq_offset) { int score = 0; // Compute score over bins 0-7, 36-43, 72-79 for (int m = 0; m <= 72; m += 36) { for (int k = 0; k < 7; ++k) { - int offset = ((i + k + m) * 4 + alt) * num_bins + j; - // score += 8 * (int)power[i + k + m][alt][j + ICOS7[k]] - + int offset = ((time_offset + k + m) * 4 + alt) * num_bins + freq_offset; + // score += 8 * (int)power[time_offset + k + m][alt][freq_offset + ICOS7[k]] - score += 8 * (int)power[offset + ICOS7[k]] - power[offset + 0] - power[offset + 1] - power[offset + 2] - power[offset + 3] - @@ -60,56 +103,23 @@ void find_candidates(int num_blocks, int num_bins, const uint8_t * power, // update the candidate list if (heap_size == num_candidates && score > heap[0].score) { - //printf("Removing score %d\n", heap[0].score); // extract the least promising candidate heap[0] = heap[heap_size - 1]; --heap_size; - // heapify from the root down - int current = 0; - while (true) { - int largest = current; - int left = 2 * current + 1; - int right = left + 1; - - if (left < heap_size && heap[left].score < heap[largest].score) { - largest = left; - } - if (right < heap_size && heap[right].score < heap[largest].score) { - largest = right; - } - if (largest == current) { - break; - } - - Candidate tmp = heap[largest]; - heap[largest] = heap[current]; - heap[current] = tmp; - current = largest; - } + heapify_down(heap, heap_size); } if (heap_size < num_candidates) { // add the current candidate - //printf("Adding score %d\n", score); heap[heap_size].score = score; - heap[heap_size].time_offset = i; - heap[heap_size].freq_offset = j; + heap[heap_size].time_offset = time_offset; + heap[heap_size].freq_offset = freq_offset; + heap[heap_size].time_sub = alt / 2; + heap[heap_size].freq_sub = alt % 2; ++heap_size; - // heapify from the last node up - int current = heap_size - 1; - while (current > 0) { - int parent = (current - 1) / 2; - if (heap[current].score >= heap[parent].score) { - break; - } - - Candidate tmp = heap[parent]; - heap[parent] = heap[current]; - heap[current] = tmp; - current = parent; - } + heapify_up(heap, heap_size); } } } @@ -117,10 +127,10 @@ void find_candidates(int num_blocks, int num_bins, const uint8_t * power, } -void extract_power(const float *signal, int num_samples, int num_bins, uint8_t * power) { +// Compute FFT magnitudes (log power) for each timeslot in the signal +void extract_power(const float * signal, int num_blocks, int num_bins, uint8_t * power) { const int block_size = 2 * num_bins; // Average over 2 bins per FSK tone const int nfft = 2 * block_size; // We take FFT of two blocks, advancing by one - const int num_blocks = (num_samples - (block_size/2) - block_size) / block_size; float window[nfft]; for (int i = 0; i < nfft; ++i) { @@ -137,31 +147,32 @@ void extract_power(const float *signal, int num_samples, int num_bins, uint8_t * kiss_fftr_cfg fft_cfg = kiss_fftr_alloc(nfft, 0, fft_work, &fft_work_size); int offset = 0; + float fft_norm = 1.0f / nfft; for (int i = 0; i < num_blocks; ++i) { // Loop over two possible time offsets (0 and block_size/2) - for (int time_offset = 0; time_offset <= block_size/2; time_offset += block_size/2) { + for (int time_sub = 0; time_sub <= block_size/2; time_sub += block_size/2) { kiss_fft_scalar timedata[nfft]; kiss_fft_cpx freqdata[nfft/2 + 1]; float mag_db[nfft/2 + 1]; // Extract windowed signal block for (int j = 0; j < nfft; ++j) { - timedata[j] = window[j] * signal[i * block_size + j + time_offset]; + timedata[j] = window[j] * signal[(i * block_size) + (j + time_sub)]; } kiss_fftr(fft_cfg, timedata, freqdata); // Compute log magnitude in decibels for (int j = 0; j < nfft/2 + 1; ++j) { - float mag2 = (freqdata[j].i * freqdata[j].i + freqdata[j].r * freqdata[j].r); - mag_db[j] = 10.0f * logf(1.0E-10f + mag2); + float mag2 = fft_norm * (freqdata[j].i * freqdata[j].i + freqdata[j].r * freqdata[j].r); + mag_db[j] = 10.0f * log10f(1.0E-10f + mag2); } // Loop over two possible frequency bin offsets (for averaging) - for (int freq_offset = 0; freq_offset <= 1; ++freq_offset) { + for (int freq_sub = 0; freq_sub < 2; ++freq_sub) { for (int j = 0; j < num_bins; ++j) { - float db1 = mag_db[j * 2 + freq_offset]; - float db2 = mag_db[j * 2 + freq_offset + 1]; + float db1 = mag_db[j * 2 + freq_sub]; + float db2 = mag_db[j * 2 + freq_sub + 1]; float db = (db1 + db2) / 2; // Scale decibels to unsigned 8-bit range @@ -177,14 +188,14 @@ void extract_power(const float *signal, int num_samples, int num_bins, uint8_t * } -int main(int argc, char **argv) { +int main(int argc, char ** argv) { // Expect one command-line argument if (argc < 2) { usage(); return -1; } - const char *wav_path = argv[1]; + const char * wav_path = argv[1]; int sample_rate = 12000; int num_samples = 15 * sample_rate; @@ -195,18 +206,54 @@ int main(int argc, char **argv) { return -1; } - const int num_bins = (int)(sample_rate / 2 / 6.25); + const float fsk_dev = 6.25f; + + const int num_bins = (int)(sample_rate / (2 * fsk_dev)); const int block_size = 2 * num_bins; const int num_blocks = (num_samples - (block_size/2) - block_size) / block_size; - uint8_t power[num_blocks * 4 * num_bins]; // [num_blocks][4][num_bins] ~ 200 KB - extract_power(signal, num_samples, num_bins, power); + printf("%d blocks, %d bins\n", num_blocks, num_bins); - const int num_candidates = 200; + extract_power(signal, num_blocks, num_bins, power); + + const int num_candidates = 250; Candidate heap[num_candidates]; - find_candidates(num_blocks, num_bins, power, num_candidates, heap); + find_sync(power, num_blocks, num_bins, num_candidates, heap); + + for (int i = 0; i < num_candidates; ++i) { + float freq_offset = (heap[i].freq_offset + heap[i].freq_sub / 2.0f) * fsk_dev; + float time_offset = (heap[i].time_offset + heap[i].time_sub / 2.0f) / fsk_dev; + // int offset = (heap[i].time_offset * 4 + heap[i].time_sub * 2 + heap[i].freq_sub) * num_bins + heap[i].freq_offset; + printf("%03d: score = %.1f freq = %.1f time = %.2f\n", i, heap[i].score / 7.0f / 2, freq_offset, time_offset); + } + + /* + + // take absolute magnitude + s2(0:7,k)=abs(csymb(1:8))/1e3 + + // skip Costas sync symbols + s1(0:7,j)=s2(0:7,k) + + // Normalize by median magnitude + s1=s1/xmeds1 + + // Extract bit significance + ps=s1(0:7,j) + bmeta(i4)=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3)) + bmeta(i2)=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5)) + bmeta(i1)=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6)) + + // Normalize by std. deviation + call normalizebmet(bmeta,3*ND) + + // Magical fudge/scale factor + scalefac=2.83 + llr0=scalefac*bmeta + + */ return 0; } \ No newline at end of file