GNU Octave Introduction

GNU Octave is a math package which has many capabilities similar to MatLAB. It is not a replacement for MatLAB, but Octave can run many MatLAB scripts and is intended to be compatible.

Getting GNU Octave

If you’re on Windows, go to the octave-forge project and get the latest .7z binary distribution of GNU Octave.

Alternative to Binaries

If you don’t want to bother downl0ading and installing GNU Octave you can run it in a browser instead. The NCLab project is entirely open source and free and allows on-line computing. It can run octave in a browser window – saving you the hassle of installing binaries on your desktop. Of course it’s a bit slower (for me anyway) and it requires a couple of different things to get plots displayed, etc. But overall it’s an excellent effort and well worth checking out. It’s very handy if all you have to do is design a quick filter to get some coefficients.

Creating Digital Signals

Lets jump straight in to it and create some signals. Firstly, let’s create a signal with two super-imposed sine-waves at different frequencies:

[code]
fs=8000;
t=0:1/fs:1;
f1=250;
sine1 = sin(2*pi*f1*t);
f2=1000;
sine2 = sin(2*pi*f2*t);
y = sine1 + sine2;
plot(t(1:100), y(1:100));
[/code]

Here, we’ve created two sine waves (sine1, and sine2) at 250Hz and 1000Hz respectively. Then we generate the sum of the waves into y (because it becomes the plots y axis). Then we plot the first 100 samples so we can see what the signal looks like in the time domain. The X axis is time.

Let’s look more closely:
Create a variable, fs which sets our sampling frequency (essential for a digital signal!)

[code]fs=8000[/code]

Create a series (an array) from 0 to 1 with an incremental change of 1/fs. This is our quantised time (or our sample points).

[code]t=0:1/fs:1;[/code]

Generate a sine-wave. sine1 is set to an array because t is a series of values. So for each point in time (sample point) we get a value of 2 * pi * f * (current time in s). This is a 250Hz sine wave.

[code]f1=250;
sine1 = sin(2*pi*f1*t);[/code]

This does the same for a 1000Hz sine wave.

[code]f2=1000;
sine2 = sin(2*pi*f2*t);[/code]

Sum both waves into a new series, y.

[code]y = sine1 + sine2;[/code]

Plot the series 1:100 of each input, i.e. time and magnitude.

[code]plot(t(1:100), y(1:100));[/code]

After the plot command you should see the sum plot, like below:

octave-test-signals1

Digital Filter Implementation

There are generally two types of Digital filter, Finite Impulse Response (FIR) and Infinite Impulse Response (IIR). Here, we’re going to go through the quick implementation of a first order (single pole) IIR. The analog equivalent of the circuit is a simple RC circuit. We’re going to implement a digital version of this simple low pass filter.

As engineers, we generally talk in terms of frequency and magnitudes of order when designing filters.

Let’s say our sampling frequency is 8kHz, and we want to have a filter that has a -3dB point of 1kHz. To start with, we need to generate some coefficients for the IIR filter we’re about to implement and analyse. We can use the most basic filter design from the DSP Book – Recursive Filters section. There, we get Equation 19-2 which describes how to obtain coefficients for a single-order low pass filter, which is:
a_o = 1 - xb_1 = x

x is defined in Equation 19-3 on the same page:
x = \epsilon^{-2 \pi f_c}

In our case where fs = 8000kHz and our cutoff frequency fc = 1000kHz = 0.125 * 8000 we get:

-2 \pi f_c = -2 \pi * 0.125
This gives us a value of
x = \epsilon^{-0.785398}
Which ends up being:
x = 0.455938

This makes our coefficients look like this:

a_0 = 1 - 0.455938 = 0.54406b_1 = 0.455938

GNU Octave Implementation of IIR

Now that we have the coefficients we require, we should eb able to put that into GNU Octave and get some results so that we can analyse the filter. IIR Filters can be unstable if they are designed wrong, and so we must always check IIR filter performance to make sure they are stable. IIR Filters have the ability to runaway with themselves, or oscilate because they include a feedback, unlike FIR Filters. It’ll also give us a chance to check the phase response of the filter too.

Here is some code which quickly evaluates the

[code]
%Start from nothing!
clear;

% Set the sampling frequency used by our digital filtering system:
fs=8000;

% Set the coefficients up (As already worked out!) for a single-pole low pass
% filter. This should give us -3dB @ 1kHz with -6dB/Octave roll off
a = [ 1 -0.455938 ];
b = [ 0.544062 ];

% Determine the frequency response of the filter design above. Get the output
% in frequency rather than rad/s. Use 64 plot points.
[H,f] = freqz(b, a, 64, fs);

% Show our cutoff frequency
cutoff = -3 * ones(64);
octaveup = -3 * ones(64);

% Plot the result so that we can see if it is correct:
figure(1);
plot(f, 20*log10(abs(H)), f, cutoff, f, octaveup);
xlabel(‘Frequency (Hz)’);
ylabel(‘Magnitude (dB)’);
[/code]

Here, we’re manually setting the filter coefficients which we’ve previously worked out to design our filter. The frequency response (with the frequency in Hz) is worked out using freqz. We then plot the output using some of Octaves Two Dimensional Plot Functions.

The output should look like this:

octave-dspbook-lowpass-response

Confirming that we do indeed have a -3dB cutoff of 1kHz and the curve does indeed imply a low pass filtering effect. -6dB per octave doesn’t look good from here though. The performance doesn’t look great for this filter though, we should be close to -9dB at 2kHz, but instead we’re closer to -6dB. It’s not a great filter!

I had some initial problems with the coefficients (as you can see by the matrix loaded with the coefficients in the above code!) and so I asked a quick question as to why things didn’t work on DSP Stack Exchange.

Hilmar’s response indicates that indeed this is a poor filter design. Let’s have a look at his low pass filter instead, where:

x = tan( f_c * \pi )

Seems simple enough. So reworking the above with that in mind gives us:

x = 0.414213

Using the matix from Hilmar too we end up with the GNU Octave coefficients of:

a[1 -x] (That’s a value of 1 for a0 and -x for a1!)

b[(1-x) (1-x)]/2

So the new Octave program we can try out is Hilmars filter, shown below:

[code]
%Start from nothing!
clear;

% Set the sampling frequency used by our digital filtering system:
fs=8000;

% Set the coefficients up (As already worked out!) for a single-pole low pass
% filter. This should give us -3dB @ 1kHz with -6dB/Octave roll off
a = [ 1 -0.414213 ];
b = [ -0.585786 -0.585786] / 2;

% Determine the frequency response of the filter design above. Get the output
% in frequency rather than rad/s. Use 64 plot points.
[H,f] = freqz(b, a, 64, fs);

% Show our cutoff frequency
cutoff = -3 * ones(64);
octaveup = -9 * ones(64);

% Plot the result so that we can see if it is correct:
figure(1);
plot(f, 20*log10(abs(H)), f, cutoff, f, octaveup);
xlabel(‘Frequency (Hz)’);
ylabel(‘Magnitude (dB)’);
[/code]

We run this filter and get the following frequency response plot:

octave-hilmar-lowpass-frequency-response

Much better! I think everyone should buy Hilmar a beer for that one! It seemed like it just popped out of his head! This filter will be much better with near infinite attenuation near the nyquist frequency.

As can be seen above, the coefficients of a good system can seem very simple and only require minimal computing effort. It is the design of the filter that is important in the first place. At least now we can test IIR filter designs (that we know the coefficients of!) in a quick and simple way. For some other filter designs, check out Robert Bristow-Johnson’s Audio Filter Cookbook.

 

One thought on “GNU Octave Introduction

  1. David

    There is different definition of a and b coefficients in Octave and in the DSP Book you refer to.
    When you compare the equation on p.572 of Octave doc (signal processing chapter) with equation 19-1 in the DSP book you will find that
    a’ = -b
    b’ = a
    with first coefficient of a = 1

    so in coefficients from referred book take form
    a0 = 1-x
    b1 = x

    you conver them to Octave definition by
    a = [ 1, -b1]
    b = [ a0 ]

Leave a Reply