%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demodulate FSK data from the tor number station. It's slow (mostly the
% matched filter for the sync signal), and I only have a hobbyist level
% of understanding of DSP and SDR so there's definitely room for
% improvement, but it seems to accurately demodulate the signal.
%
% Data transmissions need to be manually clipped out of the broadcast, and
% I recommend downsampling to 2000 Hz to speed up demodulation. Then set
% the directory and file name (minus the extension) below and run, output
% will be written to [name].txt in the same directory.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% name = '2016-08-24_00-22-25';
% name = '2016-08-24_01-07-50';
% name = '2016-08-24_01-52-47';
% name = '2016-08-24_02-38-07';
name = '2016-08-24_03-23-18';
directory = '/Users/nfmbfmpofu/projects/tor num station/download/server/clips/';
source = [directory name '.wav'];
dest = [directory name '.txt'];

% read the audio file in, and open the output file
[y,Fs] = audioread(source);
fd = fopen(dest, 'w');

fprintf(fd, '# loaded "%s.wav": %d samples @ %d Hz\n', name, length(y), Fs);

data_baud = 20;
sync_baud = data_baud * 2;

window_size = 128;
overlap = window_size - 1;
bin_width = Fs / 2 / window_size;

% TODO something doesn't work when trying a smaller overlap, and that would
% speed things up a lot if it worked
fft_samples_per_bit = Fs / data_baud / (window_size / (overlap + 1));

% define the tone frequencies 
freq_data0 = 390;
freq_data1 = 515;
freq_sync0 = 640;
freq_sync1 = 940;

% +1 because matlab is 1-indexed
bin_data0 = round(freq_data0 / bin_width) + 1;
bin_data1 = round(freq_data1 / bin_width) + 1;
bin_sync0 = round(freq_sync0 / bin_width) + 1;
bin_sync1 = round(freq_sync1 / bin_width) + 1;

% define some arbitrary signal thresholds to ignore noise
sync_peak_threshold = 0.2;
data_bit_threshold = 0.1;

% spectrogram(y,window_size,overlap,[],Fs,'yaxis');
% colormap bone;
% return;

fprintf(fd, '# sliding FFT...\n');

fft = spectrogram(y, window_size, overlap);
fft_mag = abs(fft);

% split out the sync and data signals from the FFT
fprintf(fd, '# extracting data and sync signals...\n');
sync_signal = fft_mag(bin_sync1, :) - fft_mag(bin_sync0, :);
data_signal = fft_mag(bin_data1, :) - fft_mag(bin_data0, :);

% normalize the signals
sync_signal = sync_signal / max(max(sync_signal), -min(sync_signal));
data_signal = data_signal / max(max(data_signal), -min(data_signal));

% generate an ideal sync signal to match against
samples_per_sync_tone = round(Fs / sync_baud);
gen_tone = @(hz, samples) sin(2 * pi * hz * (0:samples - 1) / Fs);
sync_waveform_known = [...
    zeros(1, samples_per_sync_tone)...
    gen_tone(freq_sync1, samples_per_sync_tone)...
    gen_tone(freq_sync0, samples_per_sync_tone)...
    zeros(1, samples_per_sync_tone)...
];

sync_fft_known = spectrogram(sync_waveform_known, window_size, overlap);
sync_fft_mag_known = abs(sync_fft_known);
sync_signal_known = sync_fft_mag_known(bin_sync1, :) - sync_fft_mag_known(bin_sync0, :);
sync_signal_known = sync_signal_known / max(max(sync_signal_known), -min(sync_signal_known));

% matched filter for the received sync signal, matched to the known sync symbol generated above
fprintf(fd, '# running sync matched filter...\n');
sync_signal_known_length = length(sync_signal_known);
sync_match_len = length(sync_signal) - sync_signal_known_length;
sync_signal_matched = zeros(1, sync_match_len);
% TODO this is painfully slow
for i = 0:sync_match_len
    sync_signal_matched(i + 1) = xcorr(sync_signal(i+1:i+sync_signal_known_length), sync_signal_known, 0);
end

% normalize the matched signal
sync_signal_matched = sync_signal_matched / max(max(sync_signal_matched), -min(sync_signal_matched));

% find where the locations the sync symbol starts
fprintf(fd, '# finding sync peaks...\n');
[sync_peaks, sync_locs] = findpeaks(sync_signal_matched, 'MinPeakHeight', sync_peak_threshold);

% hold on;
% plot(sync_signal,'linewidth',2);
% plot(data_signal,'linewidth',2);
% plot(sync_signal_matched,'linewidth',2);
% plot(sync_locs, sync_peaks, 'v');
% hold off;

% start the clock after the first sync block (sync block 2 tones @ twice
% the baud rate, so sync block length = bit length)
clock = sync_locs(1) + fft_samples_per_bit * 2;
bits = '';
bit_index = 1;
next_tick_index = 2;
signal_length = length(data_signal);

fprintf(fd, '# running clock and slicing bits...\n');

% run the baud clock
while clock <= signal_length
    % the last symbol ends with a different tone, instead of a sync block,
    % so the clock keeps ticking and recording any bits above the threshold
    is_past_last_sync = next_tick_index > length(sync_locs);
    
    if is_past_last_sync || clock - sync_locs(next_tick_index) < fft_samples_per_bit / 2
       % slice the last bit
        bit = sum(data_signal(clock - fft_samples_per_bit:clock)) / fft_samples_per_bit; 
        
        if abs(bit) >= data_bit_threshold
            if (bit > 0)
                bits(bit_index) = '1';
            else 
                bits(bit_index) = '0';
            end

            bit_index = bit_index + 1;
        end
    end

    
    % output the symbol if we've hit the next sync tick
    if ~is_past_last_sync && clock >= sync_locs(next_tick_index)
        % write symbol to file
        if isempty(bits)
            fprintf(fd, '\n');
        else
            fprintf(fd, '%s ', bits);
        end
        
        % sync the next clock tick to the end of the sync bit, it will be
        % ticked further to the end of the first bit below
        clock = sync_locs(next_tick_index) + fft_samples_per_bit;
        next_tick_index = next_tick_index + 1;
        
        % reset symbol state for the next symbol
        bit_index = 1;
        bits = '';
    end
    
    % tick the clock to the next bit
    clock = clock + fft_samples_per_bit;
end

% if there is a symbol that wasn't ended by a sync, write it to disk as well
if bit_index > 0
    fprintf(fd, '%s ', bits);
end

fclose(fd);