{"id":947,"date":"2019-08-11T12:39:35","date_gmt":"2019-08-11T11:39:35","guid":{"rendered":"https:\/\/mightydevices.com\/?p=947"},"modified":"2019-09-24T09:23:35","modified_gmt":"2019-09-24T08:23:35","slug":"decoding-fsk-transmission-recorded-by-rtl-sdr-dongle","status":"publish","type":"post","link":"https:\/\/mightydevices.com\/index.php\/2019\/08\/decoding-fsk-transmission-recorded-by-rtl-sdr-dongle\/","title":{"rendered":"Decoding FSK transmission recorded by RTL-SDR dongle"},"content":{"rendered":"\n

In this article I’m about to present the process of decoding previously recorded FSK transmission with the use of Python, NumPy and SciPy. <\/p>\n\n\n\n

The data was gathered using SDR# software and it is stored as a wave file that contains all the IQ baseband data @ 2.4Msps<\/p>\n\n\n\n

The complete source code for this project is available here: https:\/\/github.com\/MightyDevices\/python-fsk-decoder<\/a> <\/p>\n\n\n\n

Step 0: Load the data<\/h2>\n\n\n\n

The data provided is in form of a ‘stereo’ wave file where the left channel contains the In-Phase data and the right one has the Quadrature data. The wave file is using 16-bit PCM samples. For the sake of further processing let’s convert the data into the list of complex number in form: \\[z[n] = in\\_phase[n] + j * quadrature[n]\\]<\/p>\n\n\n\n

If we don’t wan’t to deal with large numbers along the way it it also wise to scale things down to \\((-1, 1)\\) according to how many bits per sample are there.<\/p>\n\n\n\n

import numpy as np\nimport scipy.signal as sig\nimport scipy.io.wavfile as wf\nimport matplotlib.pyplot as plt\n\n# read the wave file\nfs, rf = wf.read('data.wav')\n# get the scale factor according to the data type\nsf = {\n    np.dtype('int16'): 2**15,\n    np.dtype('int32'): 2**32,\n}[rf.dtype]\n\n# convert to complex number c = in_phase + j*quadrature and scale so that we are\n# in (-1 , 1) range\nrf = (rf[:, 0] + 1j * rf[:, 1]) \/ sf<\/code><\/pre>\n\n\n\n
Step 1: Center The data around the DC<\/h2>\n\n\n\n
First of all we need to tune the system so that it receives only what will become the subject of further analysis. Initial data spectrum looks like this:<\/p>\n\n\n\n
\"\"
Spectrum shows that we have some DC offset issues (noting unusual for cheap dongles) as well as the fact that we are tuned 366kHz below the band of the actual transmission<\/figcaption><\/figure>\n\n\n\n
In order to move the signal of interest to the DC we use the concept of mixing <\/strong>which is no more no less than multiplying by a sine and cosine<\/strong> that are tuned to the offset frequency. Thank’s to the magic of complex numbers we can combine the whole sine\/cosine gig into single operation using the equation \\[e^{i\\theta}=cos(\\theta) + isin(\\theta)\\] If we want to generate appropriate sine and cosine functions that represent 360kHz oscillations in 2.4Msps system then the \\(\\theta\\) has to be as follows: \\[\\theta(n)=n*2\\pi*\\frac{360kHz}{2.4MHz}\\] where \\(n\\) is the sample number, ranging from 0 to the number of samples within the input data.<\/p>\n\n\n\n
# offset frequency in Hz (read from the previous plot)\noffset_frequency = 366.8e3\n# baseband local oscillator\nbb_lo = np.exp(1j * (2 * np.pi * (-offset_frequency \/ fs) *\n                      np.arange(0, len(rf))))\n\n# complex-mix to bring the rf signal to baseband (so that is centered around\n# something around 0Hz. doesn't have to be perfect.\nbb = rf * bb_lo<\/code><\/pre>\n\n\n\n
This is what we end up with. All is shifted, including that naaasty DC spike.<\/p>\n\n\n\n
\"\"
Result of complex mixing. Everything is shifted down by the offset frequency.<\/figcaption><\/figure>\n\n\n\n
Step 2:  Remove unwanted data in the frequency domain<\/h2>\n\n\n\n
It’s easy to see that the signal of interest occupies only a small part of the band so the obvious next step will be to limit the sampling rate. The process is called Decimation<\/strong> and must be accompanied by proper filtering<\/strong>. Without filtering all the out-of-band data would simply alias <\/strong>into the band of interest. Let’s design appropriate filter, apply it and decimate (remove all samples except every n-th)<\/p>\n\n\n\n
In this example I’ve decimated by 4<\/strong> even though one could go as high as 8, but the remaining samples will come in handy when we reach the point of symbol synchronization.<\/p>\n\n\n\n
# limit the sampling rate using decimation, let's use the decimation by 4\nbb_dec_factor = 4\n# get the resulting baseband sampling frequency\nbb_fs = fs \/\/ bb_dec_factor\n# let's prepare the low pass decimation filter that will have a cutoff at the\n# half of the bandwidth after the decimation\ndec_lp_filter = sig.butter(3, 1 \/ (bb_dec_factor * 2))\n# filter the signal\nbb = sig.filtfilt(*dec_lp_filter, bb)\n# decimate\nbb = bb[::bb_dec_factor]<\/code><\/pre>\n\n\n\n
\"\"
Spectrum of decimated baseband signal. Spectrum spans for 600kHz instead of 2.4MHz which is the result of decimating by 4<\/figcaption><\/figure>\n\n\n\n
Step 3: Remove unwanted data in time domain<\/h2>\n\n\n\n
The overall length of the recording is much-much longer than the transmission itself. It would be wise to get rid of the moments of “radio silence” so that we can focus solely on what’s meaningful. Needless to say the computation performance will also benefit greatly from that.<\/p>\n\n\n\n
The easiest way around that is to select the signal that is above certain magnitude threshold.<\/p>\n\n\n\n
\"\"
Whole recording and a short ‘blip’ of actual data transmission. Data transmitted is always greater than 0.01 in terms of magnitude<\/figcaption><\/figure>\n\n\n\n
Here’s the selection code. I’ve selected samples spanning from the 1st one that exceeded 0.01 level to the last one. This is to avoid any potential discontinuities that may occur for when signal strength drops for couple of samples.<\/p>\n\n\n\n
# using the signal magnitude let's determine when the actual transmission took\n# place\nbb_mag = np.abs(bb)\n# mag threshold level (as read from the chart above)\nbb_mag_thrs = 0.01\n\n# indices with magnitude higher than threshold\nbb_indices = np.nonzero(bb_mag > bb_mag_thrs)[0]\n# limit the signal\nbb = bb[np.min(bb_indices) : np.max(bb_indices)]<\/code><\/pre>\n\n\n\n
\"\"
Selected data spans only for 10ms as opposed to complete recording that is over 2 seconds long.<\/figcaption><\/figure>\n\n\n\n
Step 4: Demodulation<\/h2>\n\n\n\n
Now it’s the perfect time to do the demodulation. The transmitter uses FSK <\/strong>which is basically a form of digitally controlled Frequency Modulation<\/strong> where ones and zeros are transmitted as tones below and above the center frequency. If you take a look at the baseband spectrum from Step 2<\/strong> you will clearly notice two peaks in the spectrum that are the result of transmitter sending consecutive zeros and ones.<\/p>\n\n\n\n
In order to know whether a zero or one is being transmitted at the moment we need to know the instantaneous frequency. Since we are dealing with complex numbers this is actually easier to accomplish than one may think. During the transmission the consecutive samples form a circle when plotted on a complex plane:<\/p>\n\n\n\n
\"\"
Data points plotted on complex plane<\/figcaption><\/figure>\n\n\n\n
In order to determine the current frequency value one needs simply to calculate the rate at which the samples rotate around the circle, meaning that we need to know the angle<\/strong> between consecutive samples. Calculating the angle between two complex numbers is quite easy: \\[angle=arg(z_{n-1}*\\bar{z_{n}})\\] where \\( \\bar{z_{n}} \\) denotes the complex conjugate<\/strong> of the n-th sample. Here’s the code for the demodulation:<\/p>\n\n\n\n
# demodulate the fm transmission using the difference between two complex number\n# arguments. multiplying the consecutive complex numbers with their respective\n# conjugate gives a number who's angle is the angle difference of the numbers\n# being multiplied\nbb_angle_diff = np.angle(bb[:-1] * np.conj(bb[1:]))\n# mean output will tell us about the frequency offset in radians per sample\n# time. If the mean is not zero that means that we have some offset. lets' get\n# rid of it\ndem = bb_angle_diff - np.mean(bb_angle_diff)<\/code><\/pre>\n\n\n\n
Second part of the snippet above helps to remove any frequency offset that may come from us not being able to provide the exact value for the frequency shift in Step 1. This operation ensures that both zeros and ones are now evenly spaced around the center frequency. And finally here’s the result of the demodulation:<\/p>\n\n\n\n
\"\"
Demodulated data, zeros and ones can now clearly be seen. Transmitter sync pattern is also apparent.<\/figcaption><\/figure>\n\n\n\n
Step 5: Guess the data rate.<\/h2>\n\n\n\n
This is as easy as looking at the spectral plot of the demodulated data. Since we are dealing with binary data transmission it should have a spectrum with shape similar to \\(H(f) = \\frac{\\sin(f)}{f}\\) Such spectra have nulls (zeros in magnitude-scale or dips in decibel-scale) at the exact multiplies of the bit frequency a.k.a bitrate<\/strong>. <\/p>\n\n\n\n
\"\"
Demodulated data spectrum. Null @ 100kHz is pretty apparent.<\/figcaption><\/figure>\n\n\n\n
From the shape of the spectrum we can make a initial guess about the bit rate to be 100kbit\/s<\/strong><\/p>\n\n\n\n
Step 6: Symbol synchronization – data recovery<\/h2>\n\n\n\n
With all the information gathered we can now proceed with data sampling. In order to determine the correct sampling time a symbol synchronization scheme must be employed.<\/p>\n\n\n\n
Such schemes (or algorithms if you like) are often build around controllable clocks (oscillators) that provide the correct timing for the sampler. The clock ‘tick’ rate is controlled by an error term derived from the algorithm. In this article I’ll discuss the simplest method: Early-Late<\/p>\n\n\n\n
Imagine that you were to sample data three times per symbol, taking the middle sample as the final symbol value. You can imagine that the best moment to sample a bit is right in the middle. That means that samples before and after the middle sample should have similar values:<\/p>\n\n\n\n
\"\"
Sampling right in the middle. First and last sample have similar values.<\/figcaption><\/figure>\n\n\n\n
If the sampling time is not aligned with the incoming stream of bits then the first and last sample are no longer similar in value. If you subtract the value of the last sample from the value of the first you’ll get the positive values when you are sampling too early (the sampling clock needs to be slowed down) and negative when sampling is late (the sampling clock needs to be sped up) with respect to incoming data. The result of the subtraction is nothing else than the error term itself.<\/p>\n\n\n\n
  • \"\"
    Data sampled too early<\/figcaption><\/figure><\/li>
  • \"\"
    Data sampled too late<\/figcaption><\/figure><\/li><\/ul>\n\n\n\n
    In order to support sign change (applying the process for the negative data from the demodulator) we additionally scale the output by the middle sample value. The whole error term is then fed to the sampling clock generator implemented as the Numerically Controlled Oscillator. This oscillator is nothing more than accumulator (phase accumulator<\/strong>, to be exact) to which we add a number (frequency word<\/strong>) once per every iteration of the algorithm and wait for it to reach a certain level (1.0 in my code). If the value is reached then we sample. Simple as that.<\/p>\n\n\n\n
    Frequency word<\/strong> (the pace at which the sampling clock is running) is constantly adjusted by the error term from above, so that the clock will eventually synchronize to the data.<\/p>\n\n\n\n
    Following code implements the early-late synchronizer. Keep in mind that I’ve intentionally negated the error value so that I can add it (in proportion) to the frequency word (nco_step and the phase accumulator (nco_phase_acc).<\/p>\n\n\n\n
    # bitrate assumption, will be corrected for using early-late symbol sync (as\n# read from the spectral content plot from above)\nbit_rate = 100e3\n\n# calculate the nco step based on the initial guess for the bit rate. Early-Late\n# requires sampling 3 times per symbol\nnco_step_initial = bit_rate * 3 \/ bb_fs\n# use the initial guess\nnco_step = nco_step_initial\n# phase accumulator values\nnco_phase_acc = 0\n# samples queue\nel_sample_queue = []\n\n# couple of control values\nnco_steps, el_errors, el_samples = [], [], []\n\n# process all samples\nfor i in range(len(dem)):\n    # current early-late error\n    el_error = 0\n    # time to sample?\n    if nco_phase_acc >= 1:\n        # wrap around\n        nco_phase_acc -= 1\n\n        # alpha tells us how far the current sample is from perfect\n        # sampling time: 0 means that dem[i] matches timing perfectly, 0.5 means\n        # that the real sampling time was between dem[i] and dem[i-1], and so on\n        alpha = nco_phase_acc \/ nco_step\n        # linear approximation between two samples\n        sample_value = (alpha * dem[i - 1] + (1 - alpha) * dem[i])\n        # append the sample value\n        el_sample_queue += [sample_value]\n\n        # got all three samples?\n        if len(el_sample_queue) == 3:\n            # get the early-late error: if this is negative we need to delay the\n            # clock\n            if el_sample_queue:\n                el_error = (el_sample_queue[2] - el_sample_queue[0]) \/ \\\n                           -el_sample_queue[1]\n            # clamp\n            el_error = np.clip(el_error, -10, 10)\n            # clear the queue\n            el_sample_queue = []\n        # store the sample\n        elif len(el_sample_queue) == 2:\n            el_samples += [(i - alpha, sample_value)]\n\n    # integral term\n    nco_step += el_error * 0.01\n    # sanity limits: do not allow for bitrates outside the 30% tolerance\n    nco_step = np.clip(nco_step, nco_step_initial * 0.7, nco_step_initial * 1.3)\n    # proportional term\n    nco_phase_acc += nco_step + el_error * 0.3\n\n    # append\n    nco_steps += [nco_step]\n    el_errors += [el_error]<\/code><\/pre>\n\n\n\n
    As the result of the algorithm we end up with sample times and their values which look like this when plotted over the demodulator output:<\/p>\n\n\n\n
    \"\"<\/figure>\n\n\n\n
    As one can tell the algorithm works nicely as there are no samples taken during the bit transitions, only when the bits are ‘stable’<\/p>\n\n\n\n
    Step 7: Determining the bit value<\/h2>\n\n\n\n
    Probably the easiest step. Just use the sampled data and produce the output value: zero or one depending on it’s sign and your done!<\/p>\n\n\n\n
    \"\"
    Decoded bit stream<\/figcaption><\/figure>\n","protected":false},"excerpt":{"rendered":"
    In this article I’m about to present the process of decoding previously recorded FSK transmission with the use of Python, NumPy and SciPy. The data was gathered using SDR# software and it is stored as a wave file that contains all the IQ baseband data @ 2.4Msps The complete source code for this project is…<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[4],"_links":{"self":[{"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/posts\/947"}],"collection":[{"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/comments?post=947"}],"version-history":[{"count":29,"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/posts\/947\/revisions"}],"predecessor-version":[{"id":990,"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/posts\/947\/revisions\/990"}],"wp:attachment":[{"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/media?parent=947"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/categories?post=947"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/mightydevices.com\/index.php\/wp-json\/wp\/v2\/tags?post=947"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}