In [1]:
# Imports and boilerplate to make graphs look better
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy
import wave
from IPython.display import Audio
def setup_graph(title='', x_label='', y_label='', fig_size=None):
fig = plt.figure()
if fig_size != None:
fig.set_size_inches(fig_size[0], fig_size[1])
ax = fig.add_subplot(111)
ax.set_title(title)
ax.set_xlabel(x_label)
ax.set_ylabel(y_label)
Now that we've got some ideas for demodulating fsk, let's do some freakin' modulation!
For starters, let's just try doubling the frequency of the carrier wave. What mathematical operation can we use to double a wave's frequency? Well, going back to high school trigonometric identities, sin(2x) = 2*sin(x)*cos(x). So if our starting wave is sin(x) and we want to double it's frequency, we can just multiply it by 2*cos(x).
In [2]:
samp_rate = 1000
len_in_sec = 1
t = np.linspace(0, 1, samp_rate * len_in_sec)
hz_4 = 1*np.sin(4 * 2 * np.pi * t)
hz_8 = hz_4 * (2 * np.cos(4 * 2 * np.pi * t))
plt.plot(t, hz_4)
plt.show()
plt.plot(t, hz_8)
plt.show()
Interesting, so just multiplying by 2*cos(x) doubles the frequency. Or, to think of it another way, just multiplying by a 90-degree phase-shifted wave of the same frequency (cos(x)... or as it could alternatively be called, sin((pi/2)-x)) doubles the frequency of a wave, but when the frequency doubles, the amplitude halves, and we have to do a scalar doubling in order to restore the original amplitude.
So let's now try to use this to modulate a carrier wave...
In [3]:
samp_rate = 1000 # samples/second
len_in_sec = 1
space_freq = 8 # Hz
t = np.linspace(0, 1, samp_rate * len_in_sec)
carrier = np.sin(space_freq * 2 * np.pi * t)
# Note: in FSK venaculare, "space" represents 0 and "mark" 1.
mark_multiplier_array = 2*np.cos(space_freq * 2 * np.pi * t)
modulation_array = np.array([1]*500 + [i for i in mark_multiplier_array[500:]])
setup_graph(title='modulation array for "01"', fig_size=(10,4))
plt.plot(modulation_array)
plt.margins(0.1)
In [4]:
setup_graph(title='modulated "01"', fig_size=(10,4))
plt.plot(carrier * modulation_array)
plt.margins(0.1)
Excellent, so that worked well.
Another simpler way to do FSK modulation could be to just calculate arrays for each of the different frequencies, and then just concatenate those arrays together. In this case, we wouldn't really be modulating a "carrier" wave so much as building the modulated wave, one bit-worth at a time.
Let's try that approach...
In [5]:
samp_rate = 1000 # samples/second
len_in_sec = 1
space_freq = 8 # Hz ("space" = "0")
mark_freq = 16 # Hz ("mark" = "1")
t = np.linspace(0, .5, samp_rate * len_in_sec)
space = np.sin(space_freq * 2 * np.pi * t)
mark = np.sin(mark_freq * 2 * np.pi * t)
modulated_01 = np.append(space, mark)
setup_graph(title='modulation array for "01"', fig_size=(10,4))
plt.plot(modulated_01)
plt.margins(0.1)
In [6]:
def bfsk_modulate(bit_array, space_freq, mark_freq, baud, sample_rate):
seconds_per_bit = 1 / baud
samples_per_bit = sample_rate * seconds_per_bit
t = np.linspace(0, seconds_per_bit, samples_per_bit)
space = np.sin(space_freq * 2 * np.pi * t)
mark = np.sin(mark_freq * 2 * np.pi * t)
signal = np.array([])
for bit in bit_array:
if bit == 0:
signal = np.append(signal, space)
elif bit == 1:
signal = np.append(signal, mark)
return signal
In [7]:
sig_010110 = bfsk_modulate([0,1,0,1,1,0], 8, 16, 1, 1000)
setup_graph(title='bfsk-modulated "010110"', fig_size=(18,5))
plt.plot(sig_010110)
plt.margins(0.05)
In [8]:
sig_010110 = bfsk_modulate([0,1,0,1,1,0], 400, 800, 5, 3200)
In [9]:
from IPython.display import Audio
Audio(sig_010110, rate=3200)
Out[9]:
In [10]:
import scipy
import scipy.io.wavfile
def write_audio_file(filename, filedata, sample_rate):
scipy.io.wavfile.write(filename, sample_rate, filedata)
write_audio_file('raw_data/bfsk_010110.wav', sig_010110, 3200)
In [11]:
def fsk_modulate(bit_str, bit_freq_map, baud, sample_rate):
seconds_per_bit = 1 / baud
samples_per_bit = int(sample_rate * seconds_per_bit)
t = np.linspace(0, seconds_per_bit, samples_per_bit)
# maps from bit sequence (like "10") to the modulated wave representing that "symbol"
symbol_map = {bit_seq: np.sin(freq * 2 * np.pi * t) for bit_seq, freq in bit_freq_map.items()}
signal = np.array([])
bits_per_symbol = len(list(bit_freq_map.keys())[0]) # Assume all keys are the same length
for symbol in [bit_str[i:i+bits_per_symbol] for i in range(0, len(bit_str), bits_per_symbol)]:
symbol_wave = symbol_map[symbol]
signal = np.append(signal, symbol_wave)
return signal
In [12]:
bit_freq_map = {
"00": 5,
"01": 10,
"10": 15,
"11": 20
}
sig_00011011 = fsk_modulate("00011011", bit_freq_map, 1, 5000)
setup_graph(title='fsk-modulated "00011011"', fig_size=(18,5))
plt.plot(sig_00011011)
plt.margins(0.05)
Now let's listen to one of audible frequencies (by just multiplying each of the frequencies by 10), and speed it up (by increasing the baud rate).
In [13]:
bit_freq_map = {
"00": 500,
"01": 1000,
"10": 1500,
"11": 2000
}
sig_00011011 = fsk_modulate("00011011", bit_freq_map, 8, 5000)
setup_graph(title='fsk-modulated "00011011"', fig_size=(18,5))
plt.plot(sig_00011011)
plt.margins(0.05)
In [14]:
Audio(sig_00011011, rate=5000)
Out[14]:
I want to be able to write a software FSK (Frequency Shift Key) Modem. First, I'll try to figure out the Demodulation step...
The simplest approach is zero-crossing detection - just see how often the sign changes to determine the frequency in a particular point.
In [15]:
samp_rate = 1000
len_in_sec = 1
carrier_freq = 20 # Hz
t = np.linspace(0, 1, samp_rate * len_in_sec)
carrier = 1*np.sin(carrier_freq * 2 * np.pi * t)
plt.plot(t, carrier)
Out[15]:
In [16]:
zero_crossings = np.where(np.diff(np.sign(carrier)))[0]
In [17]:
zero_crossings
Out[17]:
In [18]:
len(zero_crossings)
Out[18]:
In [19]:
# Note that in a single wave cycle, there will be 2 zero crossings
frequency_detected = len(zero_crossings) / 2
frequency_detected
Out[19]:
Excellent! So using zero-crossing detection, we accurately determined the frequency to be 20Hz
To try to determine the frequency, we can just take the dot product of a number of prospective frequencies, and see which has the largest dot product, which should correspond to the closest approximation of the frequency.
In [20]:
hz_10 = 1*np.sin(10 * 2 * np.pi * t)
hz_20 = 1*np.sin(20 * 2 * np.pi * t)
hz_30 = 1*np.sin(30 * 2 * np.pi * t)
hz_40 = 1*np.sin(40 * 2 * np.pi * t)
In [21]:
[np.dot(carrier, hz_10), np.dot(carrier, hz_20), np.dot(carrier, hz_30), np.dot(carrier, hz_40)]
Out[21]:
In [22]:
sum([hz_20[i]*carrier[i] for i in range(len(carrier))])
Out[22]:
Success. The 20Hz dot product is clearly much larger than the 10, 20, or 40Hz dot products, indicating the frequency is 20Hz.
It's worth noting, however, that if the frequency is just a bit off, the dot product is still near zero, so there's not a lot of room for error in terms of the frequency with this approach. Only 0.5Hz off has a near-zero dot product.
In [23]:
[np.dot(carrier, 1*np.sin(freq * 2 * np.pi * t)) for freq in [19.5, 19.6, 19.7, 19.8, 19.9, 20]]
Out[23]:
In [24]:
def dot_between_freqs(f1, f2):
t = np.linspace(0, 1, samp_rate * len_in_sec)
f1_samples = 1*np.sin(f1 * 2 * np.pi * t)
f2_samples = 1*np.sin(f2 * 2 * np.pi * t)
return np.dot(f1_samples, f2_samples)
In [25]:
center_freq = 20
frequency_diffs = np.linspace(-10, 10, 500)
dots = [dot_between_freqs(center_freq, center_freq+d) for d in frequency_diffs]
setup_graph(title='frequency deviation vs dot product', x_label='frequency deviation (in Hz)', y_label='dot product', fig_size=(14,7))
plt.plot(frequency_diffs, dots)
Out[25]:
Interesting! So it looks like a reflected damped sine wave. Now, does the variation in the dot product change based on absolute deviation from the center frequency, or is it based on percentage? Let's try a higher frequency with the same deviations...
In [26]:
center_freq = 20000
frequency_diffs = np.linspace(-10, 10, 500)
dots = [dot_between_freqs(center_freq, center_freq+d) for d in frequency_diffs]
setup_graph(title='frequency deviation vs dot product', x_label='frequency deviation (in Hz)', y_label='dot product', fig_size=(14,7))
plt.plot(frequency_diffs, dots)
Out[26]:
So it looks like the dot product varies as you move away from the center frequency purely based on absolute difference in Hz, rather than based on percentage of deviation (relative to the frequency).
The above method is kind of like a manual FFT. Let's see how the normal FFT fares.
In [27]:
fft_output = np.fft.rfft(carrier)
_ = plt.plot(fft_output)
In [28]:
len(carrier)
Out[28]:
In [29]:
[np.abs(fft_output[10]), np.abs(fft_output[20]), np.abs(fft_output[30])]
Out[29]:
Just as with the manual wave dot products above, you can see that the FFT output corresponding to 20Hz is the largest, indicating the wave is 20Hz.
Now, it's worth noting that the fft output length defaults to the num_samples/2. But for efficiency, we can constrain the fft to some lower number of output bins...
In [30]:
fft_output100 = np.fft.rfft(carrier, n=100)
plt.plot(fft_output100)
Out[30]:
And it works, but since the sample size didn't change, we need to divide the frequencies we're looking for by 10 as well (so below, fft_output100[2] corresponds to the 20Hz component.
In [31]:
[np.abs(fft_output100[1]), np.abs(fft_output100[2]), np.abs(fft_output100[3])]
Out[31]:
In [32]:
aprs_msg = "W0HAK>NY5N:>Hello World'"
aprs_msg_bits = ''.join(["{0:b}".format(ord(c)).zfill(8) for c in aprs_msg])
aprs_msg_bits
Out[32]:
In [33]:
bit_freq_map = {
"0": 2200,
"1": 1200
}
baud = 1200
sample_rate = 44100
aprs_msg_signal = fsk_modulate(aprs_msg_bits, {"0": 2200, "1": 1200}, baud=1200, sample_rate=44100)
In [34]:
len(aprs_msg_signal)
Out[34]:
In [35]:
Audio(aprs_msg_signal, rate=44100)
Out[35]:
In [36]:
def fsk_demodulate(raw_signal, bit_freq_map, baud, sample_rate):
seconds_per_bit = 1 / baud
samples_per_bit = int(sample_rate * seconds_per_bit)
t = np.linspace(0, seconds_per_bit, samples_per_bit)
# maps from bit sequence (like "10") to the modulated wave representing that "symbol"
wave_to_symbol_map = {bit_seq: np.sin(freq * 2 * np.pi * t) for bit_seq, freq in bit_freq_map.items()}
bit_str = ""
for index in range(0, len(raw_signal), samples_per_bit):
best_symbol = ""
highest_dot_abs = 0
for symbol, symbol_wave in wave_to_symbol_map.items():
raw_window = raw_signal[index:index+samples_per_bit]
dot_abs = np.abs(np.dot(symbol_wave[0:len(raw_window)], raw_window))
if dot_abs > highest_dot_abs:
best_symbol = symbol
highest_dot_abs = dot_abs
bit_str += best_symbol
return bit_str
In [37]:
demodulated_aprs_msg = fsk_demodulate(aprs_msg_signal, bit_freq_map, baud, sample_rate)
demodulated_aprs_msg
Out[37]:
In [38]:
demodulated_aprs_msg == aprs_msg_bits
Out[38]:
Aw yeah!
Continued in fsk_modem_research.ipynb.