Guest User

tor_number_station.m

a guest
Aug 26th, 2016
87
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2. % Demodulate FSK data from the tor number station. It's slow (mostly the
  3. % matched filter for the sync signal), and I only have a hobbyist level
  4. % of understanding of DSP and SDR so there's definitely room for
  5. % improvement, but it seems to accurately demodulate the signal.
  6. %
  7. % Data transmissions need to be manually clipped out of the broadcast, and
  8. % I recommend downsampling to 2000 Hz to speed up demodulation. Then set
  9. % the directory and file name (minus the extension) below and run, output
  10. % will be written to [name].txt in the same directory.
  11. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  12.  
  13. % name = '2016-08-24_00-22-25';
  14. % name = '2016-08-24_01-07-50';
  15. % name = '2016-08-24_01-52-47';
  16. % name = '2016-08-24_02-38-07';
  17. name = '2016-08-24_03-23-18';
  18. directory = '/Users/nfmbfmpofu/projects/tor num station/download/server/clips/';
  19. source = [directory name '.wav'];
  20. dest = [directory name '.txt'];
  21.  
  22. % read the audio file in, and open the output file
  23. [y,Fs] = audioread(source);
  24. fd = fopen(dest, 'w');
  25.  
  26. fprintf(fd, '# loaded "%s.wav": %d samples @ %d Hz\n', name, length(y), Fs);
  27.  
  28. data_baud = 20;
  29. sync_baud = data_baud * 2;
  30.  
  31. window_size = 128;
  32. overlap = window_size - 1;
  33. bin_width = Fs / 2 / window_size;
  34.  
  35. % TODO something doesn't work when trying a smaller overlap, and that would
  36. % speed things up a lot if it worked
  37. fft_samples_per_bit = Fs / data_baud / (window_size / (overlap + 1));
  38.  
  39. % define the tone frequencies
  40. freq_data0 = 390;
  41. freq_data1 = 515;
  42. freq_sync0 = 640;
  43. freq_sync1 = 940;
  44.  
  45. % +1 because matlab is 1-indexed
  46. bin_data0 = round(freq_data0 / bin_width) + 1;
  47. bin_data1 = round(freq_data1 / bin_width) + 1;
  48. bin_sync0 = round(freq_sync0 / bin_width) + 1;
  49. bin_sync1 = round(freq_sync1 / bin_width) + 1;
  50.  
  51. % define some arbitrary signal thresholds to ignore noise
  52. sync_peak_threshold = 0.2;
  53. data_bit_threshold = 0.1;
  54.  
  55. % spectrogram(y,window_size,overlap,[],Fs,'yaxis');
  56. % colormap bone;
  57. % return;
  58.  
  59. fprintf(fd, '# sliding FFT...\n');
  60.  
  61. fft = spectrogram(y, window_size, overlap);
  62. fft_mag = abs(fft);
  63.  
  64. % split out the sync and data signals from the FFT
  65. fprintf(fd, '# extracting data and sync signals...\n');
  66. sync_signal = fft_mag(bin_sync1, :) - fft_mag(bin_sync0, :);
  67. data_signal = fft_mag(bin_data1, :) - fft_mag(bin_data0, :);
  68.  
  69. % normalize the signals
  70. sync_signal = sync_signal / max(max(sync_signal), -min(sync_signal));
  71. data_signal = data_signal / max(max(data_signal), -min(data_signal));
  72.  
  73. % generate an ideal sync signal to match against
  74. samples_per_sync_tone = round(Fs / sync_baud);
  75. gen_tone = @(hz, samples) sin(2 * pi * hz * (0:samples - 1) / Fs);
  76. sync_waveform_known = [...
  77.     zeros(1, samples_per_sync_tone)...
  78.     gen_tone(freq_sync1, samples_per_sync_tone)...
  79.     gen_tone(freq_sync0, samples_per_sync_tone)...
  80.     zeros(1, samples_per_sync_tone)...
  81. ];
  82.  
  83. sync_fft_known = spectrogram(sync_waveform_known, window_size, overlap);
  84. sync_fft_mag_known = abs(sync_fft_known);
  85. sync_signal_known = sync_fft_mag_known(bin_sync1, :) - sync_fft_mag_known(bin_sync0, :);
  86. sync_signal_known = sync_signal_known / max(max(sync_signal_known), -min(sync_signal_known));
  87.  
  88. % matched filter for the received sync signal, matched to the known sync symbol generated above
  89. fprintf(fd, '# running sync matched filter...\n');
  90. sync_signal_known_length = length(sync_signal_known);
  91. sync_match_len = length(sync_signal) - sync_signal_known_length;
  92. sync_signal_matched = zeros(1, sync_match_len);
  93. % TODO this is painfully slow
  94. for i = 0:sync_match_len
  95.     sync_signal_matched(i + 1) = xcorr(sync_signal(i+1:i+sync_signal_known_length), sync_signal_known, 0);
  96. end
  97.  
  98. % normalize the matched signal
  99. sync_signal_matched = sync_signal_matched / max(max(sync_signal_matched), -min(sync_signal_matched));
  100.  
  101. % find where the locations the sync symbol starts
  102. fprintf(fd, '# finding sync peaks...\n');
  103. [sync_peaks, sync_locs] = findpeaks(sync_signal_matched, 'MinPeakHeight', sync_peak_threshold);
  104.  
  105. % hold on;
  106. % plot(sync_signal,'linewidth',2);
  107. % plot(data_signal,'linewidth',2);
  108. % plot(sync_signal_matched,'linewidth',2);
  109. % plot(sync_locs, sync_peaks, 'v');
  110. % hold off;
  111.  
  112. % start the clock after the first sync block (sync block 2 tones @ twice
  113. % the baud rate, so sync block length = bit length)
  114. clock = sync_locs(1) + fft_samples_per_bit * 2;
  115. bits = '';
  116. bit_index = 1;
  117. next_tick_index = 2;
  118. signal_length = length(data_signal);
  119.  
  120. fprintf(fd, '# running clock and slicing bits...\n');
  121.  
  122. % run the baud clock
  123. while clock <= signal_length
  124.     % the last symbol ends with a different tone, instead of a sync block,
  125.     % so the clock keeps ticking and recording any bits above the threshold
  126.     is_past_last_sync = next_tick_index > length(sync_locs);
  127.    
  128.     if is_past_last_sync || clock - sync_locs(next_tick_index) < fft_samples_per_bit / 2
  129.        % slice the last bit
  130.         bit = sum(data_signal(clock - fft_samples_per_bit:clock)) / fft_samples_per_bit;
  131.        
  132.         if abs(bit) >= data_bit_threshold
  133.             if (bit > 0)
  134.                 bits(bit_index) = '1';
  135.             else
  136.                 bits(bit_index) = '0';
  137.             end
  138.  
  139.             bit_index = bit_index + 1;
  140.         end
  141.     end
  142.  
  143.    
  144.     % output the symbol if we've hit the next sync tick
  145.     if ~is_past_last_sync && clock >= sync_locs(next_tick_index)
  146.         % write symbol to file
  147.         if isempty(bits)
  148.             fprintf(fd, '\n');
  149.         else
  150.             fprintf(fd, '%s ', bits);
  151.         end
  152.        
  153.         % sync the next clock tick to the end of the sync bit, it will be
  154.         % ticked further to the end of the first bit below
  155.         clock = sync_locs(next_tick_index) + fft_samples_per_bit;
  156.         next_tick_index = next_tick_index + 1;
  157.        
  158.         % reset symbol state for the next symbol
  159.         bit_index = 1;
  160.         bits = '';
  161.     end
  162.    
  163.     % tick the clock to the next bit
  164.     clock = clock + fft_samples_per_bit;
  165. end
  166.  
  167. % if there is a symbol that wasn't ended by a sync, write it to disk as well
  168. if bit_index > 0
  169.     fprintf(fd, '%s ', bits);
  170. end
  171.  
  172. fclose(fd);
RAW Paste Data