0

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 (e.g. this one and minimodem).

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.
Here's my script so far: 

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 > 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()
Evan
  • 1,960
  • 4
  • 26
  • 54

1 Answers1

1

Seems that you use multiple channles and the sound you need is embedded in one of them.

So far I have found few problems in your scripts:

  1. Nyquist rate is not a half rate of your sound. It is the rate which could sample the original sound wave, and should be at least 2 times bigger than the sound sampling rate. Hence,

    nyq = 0.5 * Fs

    is wrong.

  2. If you take advantage of the noiseless sound to demodulate, then the Differentiator can be omitted.

  3. For the low-pass filter:

    h = signal.firwin(numtaps = 10, cutoff = freq_range[1], nyq = nyq)

    the cutoff frequency is your data sample rate, please read this.

  4. filt is the final signal which can extract the specific data you desire.

  5. How to choose points in sampled_signal to recreate the original signal actually depends on the ratio between the original signal rate and the sampling rate. Just like the first link you provided, assuming the data were written in 11025 Hz and the sampling or recording rate is 44100 Hz, then the code you gave:

    sampled_signal = env[Fs/Fbit/2:params[3]+1:]

    should be:

    sampled_signal = filt[Fs/Fbit*2:params[3]:Fs/Fbit*4]

    where Fs/Fbit*2 is the beginning, params[3] is the ending, Fs/Fbit*4 is the step length.

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.

It is normal, because of different sample rates, you can add some special characters acting like a start-sign and end-sign in you text, and try to find them, then you might find the data with correct lenght you need.

Wei Yang
  • 78
  • 5