IQ Constelation Diagram of FSK signal

Question

I'm want to generate the IQ components of a FSK signal using numpy and then try to demodulate the resulting IQ samples by looking at the derivative of its phase in time.

This is the script I have:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = int(1.8e6) # Hz
rb = int(1e3) # bps

message = b"AAAABBBBCCCCDDDD"

f_one = 30e3
f_zero = 10e3
fc = (f_one + f_zero) / 2

t = np.arange(0, 0.2, 1 / fs)
x = np.zeros(np.size(t))

start_t = 0.05
start_i = int(fs * start_t)
samples_per_bit = fs // rb
print(f"Each bit takes: {samples_per_bit} samples, {samples_per_bit / fs * 1e3:.2f} ms")

bit_index = 0
for i in range(10):
    tempslice = slice(
        start_i + bit_index * samples_per_bit,
        start_i + (bit_index + 1) * samples_per_bit
    )
    if i % 2 == 0:
        x[tempslice] = np.cos(2*np.pi*f_one*t[tempslice])
    else:
        x[tempslice] = np.cos(2*np.pi*f_zero*t[tempslice])
    bit_index += 1

for byte in message:
    for i in range(8):
        one = (byte >> (7 - i)) & 1
        tempslice = slice(
            start_i + bit_index * samples_per_bit,
            start_i + (bit_index + 1) * samples_per_bit
        )
        if one:
            x[tempslice] = np.cos(2*np.pi*f_one*t[tempslice])
        else:
            x[tempslice] = np.cos(2*np.pi*f_zero*t[tempslice])
        bit_index += 1

plt.figure()
plt.title("Signal vs time")
plt.plot(t * 1e3, x)
plt.ylabel("Amplitude")
plt.xlabel("Time [ms]")

### Plot spectogram ###
fft_size = 256

num_rows = len(x) // fft_size
spectrogram = np.zeros((num_rows, fft_size))
for i in range(num_rows):
    w = x[i*fft_size:(i+1)*fft_size] * np.hamming(fft_size)
    spectrogram[i, :] = 20*np.log10(np.abs(np.fft.fftshift(np.fft.fft(w))))
spectrogram = spectrogram[:, fft_size//2:]

plt.figure()
plt.imshow(spectrogram[::-1,:], aspect='auto', extent = [0, fs/2/1e6, 0, len(x)/fs * 1e3])
plt.title("Original spectrogram")
plt.xlabel("Frequency [MHz]")
plt.ylabel("Time [ms]")

### IQ SAMPLING ###
I = x * np.cos(2*np.pi*fc*t)
Q = x * np.sin(2*np.pi*fc*t)
s = I + 1j * Q

plt.figure()
plt.plot(t * 1e3, np.real(s), label="I")
plt.plot(t * 1e3, np.imag(s), label="Q")
plt.title("After sampling")
plt.ylabel("Amplitude")
plt.xlabel("Time [ms]")
plt.legend()

plt.figure()
plt.plot(I, Q,'.')
plt.title("Constelation")

# Recover signal
plt.figure()
plt.plot(t, I*np.cos(2*np.pi*fc*t) + Q*np.sin(2*np.pi*fc*t))
plt.title("Signal from IQ")
plt.ylabel("Amplitude")
plt.xlabel("Time [ms]")

### FSK DEMODULATION ###
real = np.real(s[:-1]) * np.real(s[1:]) - np.imag(s[:-1]) * np.imag(s[1:])
imag = np.real(s[:-1]) * np.imag(s[1:]) - np.imag(s[:-1]) * np.real(s[1:])
phase_der = np.arctan2(imag, real)

plt.figure()
plt.plot(t[:-1] * 1e3, real, label="I")
plt.plot(t[:-1] * 1e3, imag, label="Q")
plt.title("After demodulation")
plt.ylabel("Amplitude")
plt.xlabel("Time [ms]")
plt.legend()

plt.figure()
plt.plot(t[:-1] * 1e3, np.sign(phase_der))
plt.title("Phase derivative")
plt.ylabel("Amplitude")
plt.xlabel("Time [ms]")
plt.show()

This is a closeup on the time signal I generate from the bitstream and its spectrogram. As you can see I chose very low frequencies for FSK so that I can see easily the effect of this modulation:

Then I try to get the IQ samples and convert it to complex numbers by:

I = x * np.cos(2*np.pi*fc*t)
Q = x * np.sin(2*np.pi*fc*t)
s = I + 1j * Q

From this I plot the constellation diagram. This is the part I don't understand. I was expecting the constellation to show a perfect circle of unit radius, but this is what I get:

Can anyone shine a light here? I reckon this has something to do with the way I calculate the IQ components, but I can't figure out what would be the correct way to do it. I also tried to recover the time signal from the IQ samples and I get exactly the same signal, so I don't know:

Answer 1

The OP is correct with the intuition that the IQ diagram for the FSK signal should remain on the unit circle. All forms of FSK are constant envelope signals, as we see clearly in the time domain plot the OP has given. A complex plot of the signal with the processing done properly will indeed be samples on the unit circle. Viewing an IQ diagram is helpful in understanding different FSK implementations, as I have demonstrated at this post. Basically, in the most general case of FSK it is viewing single tones as spinning phasors; when the signal is at baseband, the upper or positive frequency spins counter-clockwise and the lower or negative frequency spins clockwise.

The issue is the OP is working with the passband signal, not the baseband equivalent signal and has incompletely created I and Q signals such as we would do in implementation; necessary low pass filtering was not yet done. Alternatively, and more simply for analysis, we could use the hilbert command to convert to the analytic signal, move this signal to baseband by multiplying the analytic signal by $e^{-j2\pi f_ct}$ where $f_c$ represents the "carrier" which would be midway between the two FSK tones used. A plot of the result of this would be the IQ diagram with all samples on the unit circle. Note that this process, after using hilbert, is similar to what the OP has subsequently done, given Euler's formula as follows, with some critical differences:

$$e^{-j2\pi f_ct} = cos(2\pi f_ct) - j sin(2\pi f_ct)$$

The analytic signal $x_a(t)$ which is determined using hilbert(x) is itself complex, so the product with $$e^{-j2\pi f_ct}$ would involve four multiplications and two sums; in Python this can be done directly as a product of complex numbers as hilbert(x) * np.exp(-1j*2*np.pi*fc*t).

To see this graphically and why the OP's constellation is distorted consider the processing as currently done by the OP as follows:

The upper and lower diagrams are mathematically equivalent; in the lower diagram I use wider arrows to represent complex (I+jQ) datapaths. Reviewing the spectrums, we can see how there will be a higher frequency spectrum that if not filtered out will result in significant distortion:

What we are seeing in the OP's IQ diagram is the higher frequency components twirling on the constant envelope locations of the baseband signal (resulting in a "Spirograph"). If you look closely, you can make out two distinct patterns, associated with which of the two frequencies is active (if the simulation is redone with just one of each frequency, you would see each of the patterns separately).

The solution is to either filter out the higher frequency term (and when the carrier is significantly higher than the modulation, this is actually very easy to do), OR convert the real passband signal to a one-sided positive only complex signal (the analytic signal) using the Hilbert Transform given:

$$x_a(t) = x(t) + j \hat{x}(t)$$

Where $x_a(t)$ is the analytic (one-sided) signal, $x(t)$ is the original real passband signal, and $ \hat{x}(t)$ is the Hilbert Transform. (Note the hilbert function in Matlab, Octave and Python scipy.signal returns the analytic signal, the imaginary part of which is the actual Hilbert Transform).

In doing that, the resulting processing of the spectrums would appear as follows (no filtering needed!):

As a side note, I recommend NOT simulating passband signals (waste of samples), and do everything as complex baseband signals directly (so don't torture yourself with the use cosines and sines, but use $e^{j\omega t}$ to cleanly represent individual frequency tones.) This maps directly to what we would see at any carrier frequency (so in this case our carrier is $f=0$), so there is no real reason to simulate all those sinusoids at any other carrier.

Answer 2

The basic problem is that -- except in cases like MSK where a signal can be taken as either FSK or PSK -- a constellation diagram of an FSK is not done in the I/Q space.

I/Q constellation diagrams make sense for PSK or variants (such as QPSK, n-APSK, etc.) where there's some sampling instant at which the phase of a filtered I/Q signal maps to a bit decision. Then you capture the phase of the signal at that instant and plot a bunch of points to see if they cluster where they ought to.

For FSK, the demodulation is different. It appears that you've chosen a wide-ratio FSK (i.e., the difference in frequency is significantly more than the symbol rate). In that case, an optimal receiver in most cases will have a pair of filters, one at each of your bit frequencies. Then you'll find the amplitude coming out of the filters, subtract them, and sample them at the correct bit time. Your constellation diagram (which will be a valid one) will be that difference, sampled at the end of each bit.