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.
I have done some background reading to introduce myself to signal processing and FSK, and I have looked at example scripts, but I feel a little lost in the specific application of building my own demodulator in software.
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.
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.
EDIT: I have changed my original script. Now there is a problem where the mark and space frequencies (1200 Hz and 2200 Hz, respectively) have the same shape when plotted, and the output of the script is all 1s. The script works by bandpass filtering each frequency, comparing them, and deciding a 1 when mark is greater than space (based on this article).
Here's a spectrogram (from Baudline) of the audio file:
EDIT: This is my updated script with as suggested in comments:
import scipy.signal.signaltools as sigtool
from scipy import signal
from scipy.signal import butter, lfilter, blackman
import wave
import struct
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import rfft, rfftfreq
mark_f = 1200
space_f = 2200
Fs = 44100
waveFile = wave.open("Sample4_160224_mono.wav", "r")
length = waveFile.getnframes()
sig = []
for i in range(0, length):
waveData = waveFile.readframes(1)
data = struct.unpack("<h", waveData)
sig.append(data)
# Bandpass filter mark and space frequencies:
low_bpass = signal.remez(5, [.1,.109,.11,.13,.131,.14],[0,1,0])
mark_filtered = signal.lfilter(low_bpass, 1, sig)
high_bpass = signal.remez(5, [.2,.209,.2100,.23,.231,.2400],[0,1,0])
space_filtered = signal.lfilter(high_bpass, 1, sig)
# Envelope detector on mark and space bandpass filtered signals:
mark_env = np.abs(sigtool.hilbert(mark_filtered))
space_env = np.abs(sigtool.hilbert(space_filtered))
# Compare and decision:
seq = range(0, len(mark_env))
rx = []
for i in seq:
if mark_env[i] > space_env[i]:
rx.append(1)
if mark_env[i] < space_env[i]:
rx.append(0)
if mark_env[i] == space_env[i]:
rx.append(0)
