title: Demodulating an FSK signal in Python
body: <p>I am trying to write a Python script that can demodulate an FSK modulated audio file and return the data encoded in the audio.  The data being transmitted is GPS NMEA strings which are embedded as the audio channel in video files.  Basically, text is encoded with FSK modulation, and I am trying to retrieve the text using Python.  The device I am using to encode the data can also decode it, so I have been able to generate the correct output, but I need to be able to do it using software.  </p>
  
  <p>I have done some background reading to introduce myself to signal processing and FSK, and I have looked at example scripts (e.g. <a href="http://www.allaboutcircuits.com/technical-articles/fsk-explained-with-python/" rel="nofollow">this one</a> and <a href="http://www.whence.com/minimodem/" rel="nofollow">minimodem</a>).  </p>
  
  <p>I managed to write a Python script that runs successfully, although the output is incorrect.  The correct output derived from the encoding/decoding device has 8,280 raw binary (0 and 1) characters, the Python output has 1,344,786.  I think I am missing a symbol synchronizer, but I'm not sure how this works.</p>
  
  <p>My question now is: how can I add symbol synchronization to the script and/or symbol timing?  Are there better examples or explanations of how to do FSK demodulation in Python?  I would appreciate any feedback or direction.  Thank you.<br>
  Here's my script so far:  </p>
  
  <pre><code>from scipy.io.wavfile import read
  import numpy as np
  import wave
  import matplotlib.pyplot as plt
  import scipy.signal as signal
  from scipy.signal import blackman, butter
  from scipy.fftpack import fft, rfft, rfftfreq, irfft
  import scipy.signal.signaltools as sigtool
  import binascii
  
  # Read in data; 'wav' allows getting paramters, 'wav1' is actual signal data
  wavfile = 'Sample4_160224_mono.wav'
  wavfile1 = open(wavfile, 'r')
  wav = wave.open(wavfile, 'r')
  wav_1 = read(wavfile1)
  params = wav.getparams()
  N = params[3]  #Sample size
  wav1 = read(wavfile1)
  wav2 = wav1[1][0:N]
  
  duration = float(params[3] / params[2])
  n_samples = len(wav2)
  Fs = params[2]
  nyq = 0.5 * Fs  #Nyquist rate
  Fbit = (params[2]*params[0]*16)/100
  print "Fbit", Fbit
  
  # Windowing function
  w = blackman(n_samples)
  print "W is", w
  
  # FFT
  wfft = rfft(wav2 * w)
  wfft_norm = wfft/N
  wfft_norm = abs(wfft_norm[range(N/2)])
  
  # Working with frequencies...
  freqs = rfftfreq(len(wfft_norm))
  index = np.argmax(np.abs(wfft))  #Returns the index of the maximum absolute value of the windowed FFT
  freq = freqs[index]  #Returns the frequency from the above index
  freq_range = [freq - 0.01, freq + 0.01]
  freq_in_Hz = abs(freq * params[2])  #Converts the Hz
  freq_range_Hz = [abs(freq_range[0] * params[2]), abs(freq_range[1] * params[2])]
  
  # Differentiator
  diff = np.diff(wav2)
  
  # Envelope detector
  env = np.abs(sigtool.hilbert(diff))
  print "ENV", len(env)
  # Low-pass filter
  h = signal.firwin(numtaps = 10, cutoff = freq_range[1], nyq = nyq)
  filt = signal.lfilter(h, 1, env)
  
  # Signal's mean
  mean = np.mean(filt)
  
  #Do some crazy stuff to get binary **maybe wrong**
  rx_data = []
  sampled_signal = env[Fs/Fbit/2:params[3]+1:]
  
  for bit in sampled_signal:
      if bit &gt; mean:
          rx_data.append(int(1))
      else:
          rx_data.append(int(0))
  
  # Save raw binary output
  rx_data1 = ''.join(map(str, (rx_data)))
  outfile1 = open('FSK_wav6_output_binary.txt', 'w')
  outfile1.write(rx_data1)
  outfile1.close()  
  </code></pre>