Feature request: FSK modulation/demodulation support

[andresv]

It would be nice if liquid-dsp would be able to modulate/demodulate FSK signals. What do you think about it?

[jgaeddert]

I have a development branch to do this, but FSK is a tricky beast. Modulating and demodulating is pretty straightforward. Recovering data when running with channel impairments (timing, carrier/phase recovery, equalization, etc.) is much tricking with non-linear modulation schemes.

I'll try to merge in what I have and push it upstream so you can give it a try. Documentation for this module, though, is lacking.

[andresv]

I would give it a try!

Currently I am using multimon-ng for AX25 packet decoding (FSK9600). It takes input from liquid-dsp FM demodulator and decodes quite nicely if carrier offset is around zero. However because of doppler shift it almost never is. Packet decoding rate would be much improved If liquid-dsp FSK demodulator would be able to do carrier offset correction etc.

[jgaeddert]

Do you have any sample data you could share with me? e.g. over-the-air data captures I could use for testing? AX25 would be a great start.

[andresv]

This is short recording from ESTCube-1: https://github.com/cubehub/samples/blob/master/ax25_fsk9600_1024k_i16.wav It has four AX25 packets. Notice that in this recording packets are shifted about 15 kHz "right" in baseband.

You could use my doppler tool for doing baseband shifting if needed:

sox -t wav ax25_fsk9600_1024k_i16.wav -esigned-integer -b16  -r 126000 -t raw - | doppler const -s 126000 -i i16 --shift 14500 > zero_shift_ax25_fsk9600_1024k_i16.iq

This command produces raw IQ file where packets are almost in the centre. Parameter --shift says how much baseband should be "moved" in Hz.

Currently I am using this long command for doing AX25 decoding:

sox -t wav ax25_fsk9600_1024k_i16.wav -esigned-integer -b16  -r 126000 -t raw - > 1.iq; cat 1.iq | doppler const -s 126000 -i i16 --shift 14500 | demod -s 126000 -r 48000 -i i16 -o i16 --bandwidth 4500 fm --deviation 3500 --squarewave | multimon-ng -t raw -a FSK9600 /dev/stdin

Typical output:

multimon-ng  (C) 1996/1997 by Tom Sailer HB9JNX/AE4WA
             (C) 2012-2014 by Elias Oenal
available demodulators: POCSAG512 POCSAG1200 POCSAG2400 EAS UFSK1200 CLIPFSK FMSFSK AFSK1200 AFSK2400 AFSK2400_2 AFSK2400_3 HAPN4800 FSK9600 DTMF ZVEI1 ZVEI2 ZVEI3 DZVEI PZVEI EEA EIA CCIR MORSE_CW DUMPCSV
Enabled demodulators: FSK9600
demod 0.9.1

doppler 1.0.1

constant shift mode
    IQ samplerate   : 126000
    IQ input type   : i16
    IQ output type  : i16

    frequency shift : 14500 Hz
FSK9600: fm ES5E-11 to ES5EC-1 UI  pid=F0
....., ....0.......0.s%.3]...5[3......w-....-..../Xb..1..1.....Z*.B.r-W..i-.3..1..).2.C..1..../.".a,....,...I(...P(."..0..../.2.m*T..X+.s).0.b..-...n+V..q,...`*...G'Sby?'Sri?(..yL(T2.c+...b(...^(TRqL'.b]5&S.iC&T.};'.B95
FSK9600: fm ES5E-11 to ES5EC-1 UI  pid=F0
....., ............%S.I:%."5# O...#."Y3%RbA(%SR=$%..Q.".r..".b.
#..-%#Qr5."P.%.".b.$#..!.!.... O... P...#.RE.%.r}@'.ra5%Q.9;%Q.aZ)...z.Z...3...y+...9&..eP&..Q'$RR!(#..1.!.r!."Q...".!...O... .....OQ......."..56'T..8$.R92
FSK9600: fm ES5E-11 to ES5EC-1 UI  pid=F0
....., ............"P........ .R.. .R.. ....!Q2..!.r1.$.bM6$..M!#..1."P...!.!....q........
.........
.........M.....r!%.?.'..?....?.....c......~*W..>y.....C.......zI.........i...~.+..>....}.....9K..:yC...(...?.+....7.
FSK9600: fm ES5E-11 to ES5EC-1 UI  pid=F0
....., .............?:.....;....;......}....z.O..z.;...yG..68':}0w..y...\.5.;R.wh..t.W;}....y.7.>.s..9~2'..p....h....g..~.pn.?.sp72.q....mk&r.hj...nn.N.s....q0g>.}t.;0zr...y0.~.q..:.n+&..b..&fZ...sa....gk...p,...eh..h

For doing this additional tools must be installed:

• demod
• multimon-ng fork that uses 48000 input sample rate

Disclaimer: sometimes this long command stops with error:

doppler stdout.write error: Error { repr: Os(32) }, insize 8192

If it happens just try again.

[evilpete]

I find the simplest way to demodulate FSK I/Q data is to take the atan2() of I,Q and for each sample and compare it with the previous sample, if the product is > 0 then the signal is in state 1, product is < 0 then the signal is in state 0. This easily determines if I is leading Q or if Q leads I.

it is simple to visualize if you plot atan2(i, q) - atan2(prev_i, preq_q), the resulting graph with have a positive or negative slope.

This can be optimized by by knowing the symbol/baud rate and sampling data very X samples where X = is the number of samples in a transmitted bit ( X = sample_rate/ baud_rate) and re-syncing every once in while to compensate for clock drift and rounding error.

the only requirements to demodulate are that you have an idea of the symbol/baud rate and the sample freq is centered near the carrier between the high an low freq deviation

I can share some example code,

[andresv]

Disclaimer: I have not yet tried liquid-dsp FSK demodulation that was added couple of weeks ago.

evilpete, what about center frequency mismatch? Do you have any example C code snippet that also does offset correction?

[jgaeddert]

This is certainly a simple method but it is not optimum as you're opening up your receiver to a lot of noise. This is what I have used in the past to demodulate GMSK but with a non-negligible performance penalty. I'm working on an implementation of the maximum likelihood detector and it's working well although I still haven't finalized the APIs yet, nor have I figured out a great way to estimate timing and carrier frequency offsets.

On Thu, Jul 30, 2015, at 02:17 AM, Peter Shipley wrote:

I find the simplest way to demodulate FSK I/Q data is to take the atan2() of I,Q and for each sample and compare it with the previous sample, if the product is > 0 then the signal is in state 1, product is < 0 then the signal is in state 0.

it is simple to visualize if you plot atan2(i, q) - atan2(prev_i, preq_q), the resulting graph with have a positive or negative slope. This can be optimized by sampling very X samples where X = is the number of samples in a transmitted bit ( sample_rate/ baud_rate) and re-syncing every once in while to compensate for clock drift and rounding error

---

Reply to this email directly or view it on GitHub: https://github.com/jgaeddert/liquid-dsp/issues/9#issuecomment-126199571

[andresv]

Finally I started to look again into liquid FSK demodulation. I can see that there are 2 implementations: cpfskmodem_example.c and fskmodem_example.c.

I guess for demodulating 2FSK signals that are used in AX25 protocol I have to use the fskmodem object?

For curious here is a sample file that contains just one AX25 packet (2FSK, 9600 bps). singlesplit.iq is 48000sps and uses i16 format for representing IQ data (int16_t I, int16_t Q).

https://dl.dropboxusercontent.com/u/1593479/singlesplit.iq

[jgaeddert]

cpfskmodem is for continuous phase modulation and demodulation and is typically used for minimum-shift keying. Think GMSK but includes many different filter types, modulation indices, and bits/symbol. I actually demodulate these using an equalizer and treat them like linear modulation because it's near optimal. Spectrally you can't see the individual tones.

The fskmodem objects are for FSK where the modulation index is very, very large. This is typical of cheap FSK transmitters and spectrally you can actually see each individual tone. I demodulate these using FFTs as this IS the optimum way to demodulate them. Using a frequency discriminator here would have very poor performance indeed.

Thanks for the file. I'll download it and try it out soon.

[andresv]

I am also slamming together a test application to try out fskmodem and I would like to know if I am using all parameters correctly.

So the signal is 9600 bps 2FSK +-3.5 kHz deviation and IQ data sample rate is 48000 sps, signal part should be around 10 kHz or less.

What about these parameters? For 2FSK m should be 1, however I am little bit confused about other parameters. For example what k is in that case?

unsigned int m           =   1;     // number of bits/symbol
unsigned int k           =   1;     // filter samples/symbol
unsigned int num_symbols = 8192;    // number of data symbols
float        bandwidth   = 10000.0 / 48000.0;    // frequency spacing

Tip: sample file (i16 format) begins with preamble which consists of 0xAA bytes.

[jgaeddert]

The value for k (samples per symbol) is precisely related to your sampling rate relative to your symbol rate. For your case, you've sampled your signal at 48 kHz and your baud rate is 9.6 kHz. This means

k = 48000 / 9.6 = 5 samples/symbol

I'm a bit confused about AX.25, though. From the file you sent it doesn't appear that the sample file you attached is appropriate for the fskmodem. Here's the spectrum of what you sent:

Normally one would expect to see two distinct tones for the type of FSK you'd use for fskmodem objects. I think you want to use the cpfskmodem objects for this.

[andresv]

Actually I got correct data using fskmodem:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "liquid.h"

int main(int argc, char*argv[]) {
    // options
    unsigned int m           = 1;                 // number of bits/symbol
    unsigned int k           = 48000 / 9600;      // filter samples/symbol http://www.dsprelated.com/showthread/comp.dsp/132393-1.php
    float        bandwidth   = 7000.0 / 48000.0;  // frequency spacing

    // create demodulator
    fskdem dem = fskdem_create(m, k, bandwidth);
    //fskdem_print(dem);

    uint8_t iq_buffer[k*4];
    float complex c_input[k];
    uint32_t j, i, bytes_read = 0;

    FILE* iqfile;
    iqfile = fopen ("singlesplit.iq", "r");

    if (iqfile != NULL) {
        unsigned int sym_out;

        while (1) {
            bytes_read = fread(iq_buffer, 1, k*4, iqfile);
            if (bytes_read != 0) {
                // convert i16 to f32
                for (j=0, i=0; j<bytes_read/2; j+=2, i++) {
                    c_input[i] = ((int16_t*)iq_buffer)[j] / 32768.0 + ((int16_t*)iq_buffer)[j+1] / 32768.0 * I;
                }

                // fwrite((uint8_t*)c_input, 8, bytes_read/4, stdout);

                sym_out = fskdem_demodulate(dem, c_input);
                printf("%01X", sym_out);
            }
            else {
                break;
            }
        }

        fclose (iqfile);
    }
    else {
        printf("cannot open file!\n");
    }

    return 0;
}

Output bits:

11111111110101010101010101010101010101010101010101010101100111111111100000111011111111011110001011011101001000010100000100010010100...

I noticed that there is 010101 pattern which in hex is 0xAA. So in python we can do:

In [32]: hex(0b11111111110101010101010101010101010101010101010101010101100111111111100000111011111111011110001011011101001000010100000100010010100)
Out[32]: '0x7feaaaaaaaaaaacffc1dfef16e90a0894L'

From there we can see preamble bytes 0xAA and sync word 0xCFFC1D (I know that because I wrote ESTCube-1 COM firmware), after that comes G3RUH scrambled and NRZI encoded data.

So fskmodem works very nicely and I guess we could close this ticket?

[jgaeddert]

Ah, good to hear! Don't forget to destroy the fskdem object at the end when you're done, otherwise you'll leak that memory out.

fskdem_destroy(dem);

I'm going to close this issue, but if you (or anyone else) is having trouble, feel free to open a new one.

[andresv]

Some more information about AX25 decoding. It seems that fskmodem is able to demodulate some of the bits quite well from this sample file. However I am not able to get all bits correct. I tried to modify bandwidth from 2000 to 8000 and found out that 5000 gave least errors. After that I also tried to modify datarate from 9000 to 10000, however error bit count stayed the same. Test script here.

Now I will try to do the same thing using cpfskmodem.

PS: why fskmodem does not have deviation parameter?

[andresv]

Could you tell me if I am on right track defining parameters for cpfskmodem. First set of parameters are taken from fskmodem and below are cpfskmodem parameters that should work eventually with 2FSK 9600 bps 3500 Hz deviation.

Actually the first problem is with k (filter samples/symbol). It cannot be odd number in case of cpfskmodem. Therefore I guess I cannot use IQ file with 48000 sps straight away, because 48000/9600 = 5 which is odd number. One of the options for me is to resample it to some other value like 38400 sps, however maybe there are some other ways?

// fskmodem options
unsigned int m           = 1;                 // number of bits/symbol
unsigned int k           = 48000 / 9600;      // filter samples/symbol 
float        bandwidth   = 5000 / 48000.0;   // frequency spacing

// cpfskmodem options
unsigned int    bps         = 1;                 // number of bits/symbol
float           h           = 0.73f;             // modulation index https://en.wikipedia.org/wiki/Frequency_modulation#Modulation_index 1/9600.*2*3500
unsigned int    k           = 48000 / 9600;      // filter samples/symbol
unsigned int    m           = 3;                 // filter delay (symbols)
float           beta        = 5000 / 48000.0;    // GMSK bandwidth-time factor

[jgaeddert]

Ah, I remember. The cpfskdem object is incomplete. I have some working code, but this hasn't been ported into the main modem processing code yet. Eventually, though, you would want to resample down to 2 or 4 samples per symbol because this would give you the best processing efficiency (lower sample rate is always better). The arbitrary resampler object (resamp_crcf) helps with this.

I'll try to have a working version for you soon.

[rodrigo455]

@jgaeddert Do you have any work done on this non-coherent cpfsk receiver? I was able to receive a cpfsk signal using a freqdem followed by a symsync, but I was not happy with its performance. I was not getting all the frame syncs. Some adjustments on the symsync loop filter coefs might be required.

[jgaeddert]

@rodrigo455 For non-coherent cpfsk that's probably the right choice; however with bursts you probably don't want to be running all the received samples through the symsync object. This will cause it to drift way off time as you feed noise through. You're better off looking for a special symbol sequence/pattern to correlate against to detect a frame and recover timing information.

[rodrigo455]

first of all thanks for the reply. Could qdetector object be used for that? There's a function that creates it for gmsk which is also non-linear. I would just need to generate a symbol sequence and then create such object.

[jgaeddert]

It certainly can be, but I haven't tried it. I don't know too much about the signal format to know what the right parameters are, though. There are a few things I should mention:

1. The GMSK parameters of the qdetector create method might not match the CP-FSK signal @andresv referenced (GMSK has h=0.5 IIRC, and it wouldn't work with h=0.73). To be fair, you can give the qdetector object any signal to correlate against and it will do so; it's just a bit of a pain to generate the template in a buffer rather than having it do the interpolation etc. for you
2. non-coherent demod of CP-FSK is easier than coherent (qdetector would help with coherent demod), but you incur a few dB performance penalty

A non-coherent demod might look something like this:

• Receive complex baseband I/Q samples in buffer, oversampled by about 4
• Perform coarse carrier frequency offset correction (if possible) and apply low-pass filter
• Run signal through FM discriminator (e.g. fmdemod object) to get instantaneous real-valued frequency representation
• Filter and re-sample to 2 samples/symbol
• correlate result against known pattern for detection and timing recovery
• get bits out

[jgaeddert]

In 6428ef62d73910f I added a method to qdetector to create an object based on a CP-FSK input. It seems to correlate pretty well to the example signal that @andresv provided. Getting the recovered signal back isn't particularly difficult, but it isn't elegant. I'm trying to think of a good way to use the qdetector to not only detect, but also provide a stream of samples to include a coarse frequency/phase and timing correction.

[Mas313]

I find the simplest way to demodulate FSK I/Q data is to take the atan2() of I,Q and for each sample and compare it with the previous sample, if the product is > 0 then the signal is in state 1, product is < 0 then the signal is in state 0. This easily determines if I is leading Q or if Q leads I.

it is simple to visualize if you plot atan2(i, q) - atan2(prev_i, preq_q), the resulting graph with have a positive or negative slope.

This can be optimized by by knowing the symbol/baud rate and sampling data very X samples where X = is the number of samples in a transmitted bit ( X = sample_rate/ baud_rate) and re-syncing every once in while to compensate for clock drift and rounding error.

the only requirements to demodulate are that you have an idea of the symbol/baud rate and the sample freq is centered near the carrier between the high an low freq deviation

I can share some example code,

Although, its to late to ask, can you share the code . Also tell if it is robust enough for practical applications in presence of noise or phase miss-alignments. thanks