Software Defined Radio for the Radio Amateur

2018-Feb-18 - Tutorial 1: Oscillators

The Java files for this tutorial are available on github.com/ac2cz/SDR
The Python files are on github.com/ac2cz/pySdr

Perhaps the first homebrew project I built that was not directly from a kit was an oscillator. They are a key part of the radio architecture, perhaps as a VFO or as a crystal oscillator for a mixer. This is a good place for us to start our adventure into software defined radio.

Like homebrew electronics, we need to be able to measure and see the results of our experiments. On the test bench we have a DVM, an oscilloscope and perhaps a spectrum analyzer. We also use a radio (software defined or otherwise) to hear and see our signals.

Let's create an oscillator and a piece of test equipment that we can use to measure it.

An analog oscillator typically has a tuned circuit and we pick C and L values to set the frequency. In the digital realm we do not produce a continuous wave, but instead deal with samples. Each sample is represented by a number. We will use values between -1 and 1 to represent our signals. You will see other formats and we will have to convert to/from them for storage in files or to pass in and out of hardware sound cards. But more on that later.

Test Tools

We need to be able to see the results of what we do. One approach is to print out values and graph them. This will get us started. Given everything is digital we can cut/paste and graph with a spreadsheet. Later we will add some code to draw graphs on the screen and eventually use that as part of our radio.

Building the Oscillator

So we start with a simple oscillator. Let's say we want a sine wave at 1200Hz.

In software our sine wave generator is called a Numerically Controlled Oscillator or NCO. In Java or Python we can create this with a short class called Oscillator which has one constructor that takes the requested frequency, the sample rate in samples per second and then builds a table of sine values. This will save us from calculating the sine value every time we need a sample of our waveform.

The frequency is relative to the sample rate and is used to set the amount that the phase increases for each sample. One complete sample of our waveform has 2 pi radians. We define the phase increment as 2 pi times frequency / samples per second. That means our 1200Hz wave will make one complete cycle in 48000 / 1200 = 40 samples

Finally, we create a method called nextSample (note that in Python I will follow the typical naming convention using lower case and underscores for functions/methods, so it will be next_sample). We will call nextSample each time we need a new sample of our waveform.

Pulled together in a class it looks like the below.

If you want to work in Python, click the bar below to show the example Python code, otherwise click Java.

import math

class Oscillator:
    TABLE_SIZE = 128

    def __init__(self, frequency, samples):
        self.freq = frequency
        self.sampleRate = samples
        self._phase = 0
        self._sinTable = self._init_table()
        self._phaseIncrement = self._calc_phase_inc(self.freq, self.sampleRate)
        
    def _init_table(self):
        table = []
        for n in range(0, self.TABLE_SIZE):
            table.append(math.sin(n*2.0*math.pi/self.TABLE_SIZE))
        return table
    
    def _calc_phase_inc(self, frequency, samplesPerSecond):
        return 2 * math.pi * frequency / samplesPerSecond

    def next_sample(self):
        self._phase = self._phase + self._phaseIncrement
        if (self._phase >= 2 * math.pi):
            self._phase = 0;
        idx = (int)((self._phase * self.TABLE_SIZE/(2 * math.pi))%self.TABLE_SIZE)
        value = self._sinTable[idx]
        return value

package tutorial1.signal;

public class Oscillator {

	private final int TABLE_SIZE = 128;
	private double[] sinTable = new double[TABLE_SIZE];
	private int samplesPerSecond = 48000;
	private int frequency = 0;
	private double phase = 0;
	private double phaseIncrement = 0;
	
	public Oscillator(int samples, int freq) {
		this.samplesPerSecond = samples;
		this.frequency = freq;
		phaseIncrement = 2 * Math.PI * frequency / samplesPerSecond;
		for (int n=0; n<TABLE_SIZE; n++) {
			sinTable[n] = Math.sin(n*2.0*Math.PI/TABLE_SIZE);
		}
	}
	
	public double nextSample() {
		phase = phase + phaseIncrement;
		if (phase >= 2 * Math.PI)
			phase = 0;
		int idx = (int)((phase * TABLE_SIZE/(2 * Math.PI))%TABLE_SIZE);
		double value = sinTable[idx];
		return value;
	}
}

Note: there is a bug in this code which we discover and fix in Tutorial 5. It is fine for this tutorial, but don't cut/paste this code for other purposes without reading about the bug.

Testing the Oscillator

Having created the oscillator, let's write a very short program to call it.

In our test we will simply create an instance of the Oscillator with the sample rate set to 48000 and the frequency set to 1200. We will then call it in a loop to output the values that are generated. We can then graph the values in a spreadsheet to see if they look right.

import signal

def main():
    buffer = []
    sample_len = 100
    osc = signal.Oscillator(1200, 48000)
    for n in range(0, sample_len):
        buffer.append(osc.next_sample())
        print(buffer[n])

main()

package tutorial1;
import tutorial1.signal.Oscillator;

public class SdrTutorial1 {

	public static void main(String[] args) {

		Oscillator osc = new Oscillator(48000, 1200);
		
		int len = 100;
		for (int n=0; n< len; n++) {
			double value = osc.nextSample();
			System.out.println(value);
		}
	}
}

When we run the program it prints out 100 values of the waveform. I pasted them into excel and while the values were still highlighted I clicked Insert > Scatter Plot Chart, which plotted the values. The commands may be different in another spreadsheet, but you should be able to plot them.

What we get looks great! We see two complete cycles of the waveform in the first 80 samples. But is that what we expect?

At 48000 samples per second, a 1200Hz tone would take 48000 / 1200 = 40 samples for a single cycle. So this looks good! It represents a 1200Hz tone at 48000 samples per second. This is the experimental method applied to SDR. We have created our oscillator, run it in a test situation and measured the results. The results agree with our theory. This is often not the case of course. Typically when I create some new code it does not do what I expect. So this experiment and measure approach is critical to make sure it is working.

One final note. If we arbitrarily decide this set of numbers are samples at 96000 samples per second, then this is a 2400Hz tone because it completes one cycle in 40 samples. 96000/40 = 2400. There is no change to the numbers or the file. Everything in DSP is relative to the sample rate. This can be confusing if you are new. There is nothing in the numbers that "stores" the sample rate in some way.

Analyzing the signal in the frequency domain

Let's finish by running some analysis on the signal using one more test tool. We will often want to look at the frequencies in our signal both to understand if our radio is working and to display the results to the user. On the electronics test bench we would use a spectrum analyzer. In software we use a Fourier Transform.

There are lots of explanations of the Fourier Transform. I like the detailed explanation in The Engineers and Scientists Guide to DSP. In fact the whole book is recommended and is available online for free.

We will use the Fast Fourier Transform (FFT) which is an algorithm that efficiently implements the Discrete Fourier Transform or DFT. Don't worry too much about the names. The DFT is just the Digital equivalent of the Fourier Transform. The FFT is not another different transform, it is just a way to implement it.

The Fourier Transform decomposes the signal into a set of component sine and cosine waves. Specifically the numbers that it returns are the magnitudes of those sine and cosine waves. The sine and cosine waves span the frequencies from 0 to half our sampling rate. The magnitude of each sine and cosine represents the magnitude of our signal at those frequencies. In fact the sines and cosines collectively are our signal. We can literally add them back together and get our signal again.

This decomposition allows us to analyze the signal as its component frequencies and to visualize them. We get to see the spectrum. We have all seen this either on a spectrum analyzer, on a pan adapter or on the waterfall of an SDR.

Before we jump in, we need to understand the output that we will get. We will supply 128 samples and the real FFT will return 65 results. The results are in pairs made up of Real (Re) and Imaginary (Im) values. The real parts are the amplitudes of the cosine waves that make up our signal and the imaginary parts are the amplitudes of the sine waves. They are called Real and Imaginary to match the algorithm for the Complex DFT, which returns Complex numbers, even though the results from the Real DFT are real numbers. For the real DFT this is just a convenient way of returning the results, but it is confusing. Don't worry about it for now, just think of it as a way to pack and return the data from the algorithm.

We don't typically look at the Cosine (Real) and Sine (Imaginary) magnitudes when we are doing frequency analysis. We are instead interested in how we can use them to calculate the actual Magnitude of our signal at each frequency and the corresponding Phase. We can calculate those from the Cosine (Real) and Sine (Imaginary) magnitudes as follows:

	mag[k] = sqrt( Re[k]^2 + Im[k]^2)
	phase[k] = arctan (Im[k]/Re[k])

If that is also confusing, then don't worry. You have to understand the results of this and not necessarily all of the mechanics. If you want to understand more, then it is well explained in Chapter 8 of the Engineers and Scientists guide to DSP.

Continuing our experimental method we expand our main routine so that we can calculate and print the magnitude of the FFT. We store the calculated NCO values in a buffer and then run our FFT with the values in the buffer. We change the length to 128 because the FFT works best with lengths that are a power of 2. You can also see that we calculate the magnitude of each frequency using two results, rather than display the actual returned values.

import matplotlib.pyplot as plt
import signal
import numpy as np
import math

def main():
    buffer = []
    sample_len = 128
    osc = signal.Oscillator(1200, 48000)
    for n in range(0, sample_len):
        buffer.append(osc.next_sample())
        print(buffer[n])
    #plt.plot(buffer, label='1200Hz')
    out = np.fft.rfft(buffer)
    psd = []
    psd.append(out[0].real)
    half = int(len(out))
    for k in range(1, half):
        psd.append(math.sqrt(out[k].real*out[k].real + out[k].imag*out[k].imag))
    psd[half-1] = math.sqrt(out[0].imag*out[0].imag + out[half-1].imag*out[half-1].imag)
    plt.plot(psd, label='FFT')
    plt.show()
main()

	public void realForward(double[] a)
	Computes 1D forward DFT of real data leaving the result in a. The physical layout of the output data is 
	as follows:

	if n is even then:

	 a[2*k] = Re[k], 0<=k<n/2
	 a[2*k+1] = Im[k], 0<k<n/2
	 a[1] = Re[n/2]

	This method computes only half of the elements of the real transform. The other half satisfies the 
	symmetry condition. If you want the full real forward transform, use realForwardFull. To get back the 
	original data, use realInverse on the output of this method.

package tutorial1;
import org.jtransforms.fft.DoubleFFT_1D;
import tutorial1.signal.Oscillator;

public class SdrTutorial1 {

	public static void main(String[] args) {

		Oscillator osc = new Oscillator(48000, 1200);
		
		int len = 128;
		double[] buffer = new double[len];
		for (int n=0; n< len; n++) {
			double value = osc.nextSample();
			buffer[n] = value;
		}

		DoubleFFT_1D fft;
		fft = new DoubleFFT_1D(len);
		fft.realForward(buffer);
		System.out.println(buffer[0]); // bin 0
		for (int k=1; k<len/2; k++)
			System.out.println(Math.sqrt(buffer[2*k]*buffer[2*k] + buffer[2*k+1]*buffer[2*k+1]));
		System.out.println(Math.sqrt(buffer[1]*buffer[1] + buffer[len/2]*buffer[len/2])); // bin n/2

	}
}

I ran the program, then cut/pasted the calculated values into a spreadsheet and plotted them as a scatter plot. It is shown below:

Wow! That is pretty exciting. We have a set of 65 result values, which we call "bins", from 0 to 64. We ran an FFT with 128 points, which returned 65 cosine (real) and 65 sine (imaginary) values. We have calculated the magnitude for each bin from the combination of the cosine (real) and sine (imaginary) values. We then see that bin 3 (counting from bin 0) is dramatically raised.

What do the bins mean? Bin 0 holds the DC component of the signal. Bin 1 represents a cosine and a sine wave that each make one cycle in the FFT length. So bin 1 represents a cosine and sine wave that are 128 samples long at 48000 samples per second, so they have 48000/128 = 375 cycles in a second, which is 375Hz. That means bin 3 represents the magnitude of our signal at 3 * 375 = 1125Hz and bin 4 shows the magnitude at 1500Hz. That seems pretty convincing for our 1200Hz tone. We have plotted the spectrum of our oscillator!

With that said, why are many other bins raised up? This might surprise us, given we generated a very pure single frequency in software. Shouldn't it be infinitely thin? Instead we see numerous bins around our target frequency raised. We won't be able to understand why without beginning our tumble down the rabbit hole of DSP.

Note: Don't use this oscillator, because it has a bug, which we will discover and fix in Tutorial 5. No need to post about it in the comments below :)

Index | Next Tutorial

Enter Comments Here:

Comments on this post