用Python实现FSK信号的解调

Question

我正在尝试编写一个Python脚本，它可以解调FSK调制的音频文件并返回音频中编码的数据。正在传输的数据是GPS-NMEA字符串，这些字符串作为音频通道嵌入到视频文件中。基本上，文本是用FSK调制编码的，我尝试使用Python检索文本。我用来编码数据的设备也可以解码，所以我已经能够生成正确的输出，但是我需要能够使用软件来完成。在

我做了一些背景阅读来介绍自己信号处理和FSK，并查看了示例脚本（例如this one和minimodem）。在

我设法编写了一个成功运行的Python脚本，尽管输出不正确。编码/解码设备的正确输出有8280个原始二进制（0和1）字符，Python输出有1344786个。我想我错过了一个符号同步器，但我不确定这是如何工作的。在

我现在的问题是：如何将符号同步添加到脚本和/或符号计时？有没有更好的例子或解释如何在Python中进行FSK解调？如有任何反馈或指导，我将不胜感激。
以下是我目前为止的剧本：

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()  

Answer 1

似乎你使用多个通道，你需要的声音嵌入其中一个。在

到目前为止，我在你的脚本中发现了一些问题：

1. 奈奎斯特速率不是你声音的一半。它是能够对原始声波进行采样的速率，至少应大于声音采样率的2倍。因此

   nyq = 0.5 * Fs
   

   是错误的。

2. 如果利用无噪声的声音进行解调，则可以省略微分器。

3. 对于低通滤波器：

   ^{pr2}$

   截止频率是您的数据采样率，请阅读this。

4. filt是最终的信号，它可以提取您想要的特定数据。

5. 如何选择采样信号中的点来重建原始信号实际上取决于原始信号速率与采样率的比值。就像您提供的第一个链接一样，假设数据以11025赫兹写入，采样或记录速率为44100赫兹，则您给出的代码：

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

   应该是：

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

   其中，Fs/Fbit*2是开始，params[3]是结束，Fs/Fbit*4是步长。

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.

这是正常的，因为不同的采样率，你可以在你的文本中添加一些类似于开始符号和结束符号的特殊字符，并尝试找到它们，这样就可以找到所需长度正确的数据。在