Discussion Forums

Using GNU Octave

19 replies [Last post]
Anders Feder
Offline
Joined: 2010-04-22
Posts: 618

As you may have read, SETI Institute is opening a sandbox enviroment in the Amazon cloud to let people test out signal processing algoritms on real SETI data. The service will be based on GNU Octave (a Matlab clone) which can be downloaded here. While we wait for cloud service to open, can we get some examples of how to use GNU Octave for signal processing? If we also download the currently available data files, it should be possible to work on them, right? The syntax should be pretty much the same as Matlab, but since I have no Matlab experience, I have trouble just doing a basic Fourier transform. Any recommendations?

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
Here's a textbook example I

Here's a textbook example I found, though not working on actual data:

Define a basic sine wave; x(k)=3sin(2πk/10) as k=1,2,3...256; Type

 octave:51>for k=1:256;s(k)=3*sin(2*pi*k/10);end

Generate some Gaussian noise

 octave:52>n=2*randn(1,256);

Add it to the signal

 octave:54>x=s+n;

Take the magnitude of the Fourier transform of the time
series

 octave:55>p=abs(fft(x));

Plot the result

 octave:56>plot(p)

This generates a sine wave and adds it to a series of Gaussian noise, much like you do when you transmit a radio signal in a noisy medium like the atmosphere or indeed space.

I don't know how to interpret the exact figures involved, though; the plot shows spikes at around 27 and 231 on the lower axis. What could that signify? Apparently, the axises are time and amplitude. Can someone please explain how t=27 and t=231 relates to the example? Also, why does the plot render over the time domain as opposed to the frequency domain?

For convenience, the complete code as it is typed in:

 for k=1:256;s(k)=3*sin(2*pi*k/10);end
 n=2*randn(1,256);
 x=s+n;
 p=abs(fft(x));
 plot(p)

robackrman
Offline
Joined: 2010-04-15
Posts: 235
I have annotated your plot to

I have annotated your plot to help with its interpretation. Because the input was real (not complex), the negative frequencies mirror the positive frequencies. There is no time axis. Time is collapsed.

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
Thanks! What is the unit on

Thanks! What is the unit on the horizontal axis then? Because it doesn't look like hertz. It is related somehow to the number of samples since the graph cuts off at 256.

robackrman
Offline
Joined: 2010-04-15
Posts: 235
The horizontal axis in this

The horizontal axis in this case is "Cycles per 256 Samples."

To calibrate with Hertz a sample rate must be defined. For example, if these samples (x=s+n) were output through a DAC at the rate of 1000 samples-per-second the embedded sinusoid (s) would have frequency 100 Hertz (cycles-per-second).

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
I am trying to make a script

I am trying to make a script that reads a number of values from a data file and plots their probability distribution. So far, I've done this by reading the data as unsigned 8-bit integers and incrementing a matrix k for each appearance of a given value. For instance, if the value 32 appears 98 times in the data, then k(32) = 98 when all the data has been processed.

However, the values in the data files are actually signed 8-bit integers. I can't use the above approach because Octave won't let me use negative indices (i.e. k(-32) = 49 is an illegal operation).

What is a better way to plot the distribution of a sequence of signed 8-bit integers?

robackrman
Offline
Joined: 2010-04-15
Posts: 235
Here is an Octave example

Here is an Octave example that plots histograms of the real and imaginary coefficient values from the first 1000000 complex samples of a setiQuest observation raw data file:

rackrman@dev:~/temp/anders$ octave -q
octave:1> N = 1000000;
octave:2> fd = fopen("2010-03-19-exo060-8bit-1-of-10.dat");
octave:3> signed_bytes = fread(fd, N*2, "schar");
octave:4> unsigned_real = signed_bytes(1:2:N*2) + 128;
octave:5> unsigned_imag = signed_bytes(2:2:N*2) + 128;
octave:6> real_hist = zeros(256, 1);
octave:7> imag_hist = zeros(256, 1);
octave:8> for k = 1:N
> real_hist(unsigned_real(k))++;
> imag_hist(unsigned_imag(k))++;
> endfor
octave:9> subplot(2,1,1);
octave:10> plot(-128:127, real_hist, "g")
octave:11> subplot(2,1,2);
octave:12> plot(-128:127, imag_hist, "r")
octave:13>

and, the resulting plot which shows the expected Gaussian distributions.

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
Thanks!

Thanks! At first I thought one curve looked wider than the other on your plot, but I guess that must be optics.

robackrman
Offline
Joined: 2010-04-15
Posts: 235
Here is an image of the

Here is an image of the histograms plotted over each other made using the previous script and these commands:

octave:14> plot(-128:127, real_hist, "g")
octave:15> hold on
octave:16> plot(-128:127, imag_hist, "r")
octave:17> set(gca(), "xlim", [-50,50])
octave:18>

The real and imaginary distributions appear to coincide - neither "wider" than the other.

Which observation data did you use in your analysis which showed an unexpected distribution in the difference between the real and imaginary coefficients? We should investigate and understand that result.

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
Yes, you're right. It was the

Yes, you're right. It was the red color that made the lower curve appear more 'compressed' than the bright green one :)

I used the first second of the Kepler-4b data (2010-01-22-kepler-exo4-1420mhz-one-second.dat). It's very likely that it was due to a bug in my code of course, but I can't immediately spot it.

maxs-pooper-scooper
Offline
Joined: 2010-07-28
Posts: 58
or you can use: hist()

FYI: octave has a histogram command

<code>
octave:1> N= 1000000;
octave:2> data= randn(N,1)+j*randn(N,1);
octave:3> [xi,yi]= hist( real(data), 100, 1.0 );
octave:4> [xq,yq]= hist( imag(data), 100, 1.0 );
octave:5> plot( yi,xi,'g', yq,xq,'r' )
</code>

You may or may not have run across it yet!

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
I sometimes have a feeling

I sometimes have a feeling that Octave keeps some kind of state between command entries that is not immediately apparent to me. Is there some way to clear that state? I.e. a way to reset all variables and so on (short of actually restarting the program, which get a little tedious).

maxs-pooper-scooper
Offline
Joined: 2010-07-28
Posts: 58
clear -- might be what you are looking for

If you run a function, all the variables have the same scope as the function, i.e. once the function exits, the variables gone.

The workspace keeps the variables in scope by default.
The command 'who' lets you know what variables are declared.
The command 'whos' lets you know what variables are declared and how much memory they consume.
The command 'clear' will remove all variables.

You might have run across these commands by now.

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
I've created a

I've created a general-purpose function for creating Gaussian white noise from a uniform entropy source (such as data downloaded from QRBG) based on the Box-Muller transform. I figured it would have applicability for other setiQuesters, so I am posting it below. Comments and corrections are welcome.

# Generates Gaussian white noise from a uniform entropy source file.
function Z0, Z1 = gwn ( blocksize, filename )
  persistent fd;
  if (exist("filename"))
    if (isscalar(fd))
    fclose(fd);
  endif
  fd = fopen(filename);
  endif
  if (!exist("blocksize"))
    blocksize = 1;
  endif
  b1 = fread(fd, blocksize, "uint8");
  b2 = fread(fd, blocksize, "uint8");
  for n = 1:blocksize
    u1 = (b1(n)+1)/256;
    u2 = (b2(n)+1)/256;
    z0 = sqrt(-2*log(u1))*cos(2*pi*u2);
    z1 = sqrt(-2*log(u1))*sin(2*pi*u2);
    Z0(n) = z0;
    Z1(n) = z1;
  endfor
endfunction

If you save the function in a file of its own called gwn.m, Octave will be able to autoload it when you invoke the function.

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
When using the "getdata"

When using the "getdata" example in setiCloud, and try to use the built-in kurtosis function on the "data" variable, I get an error saying that it is not a matrix. What should I do to convert this variable to a matrix?

robackrman
Offline
Joined: 2010-04-15
Posts: 235
Is this the error?: "error:

Is this the error?:

"error: kurtosis: x has to be a matrix or a vector"

I am not familiar with the example; however, I suspect getdata is returning a string of 8-bit values, not a vector.  Try the following:

vctr = x(1:length(x))*1.0;
kurtosis(vctr)

substitute your variable (returned by getdata) for x.

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
Yes, that did the trick.

Yes, that did the trick. Thanks. I wonder what exactly your code does though? It looks like it is just applying the identity operator to every element in the string?

robackrman
Offline
Joined: 2010-04-15
Posts: 235
Yeah, it's a trick I use for

Yeah, it's a trick I use for converting a string to an array.  See below for more detail.  There may be a different and better way to do it.
rackrman@dev:~$ octave -q octave:1> str="abcdefghijklmnopqrstuvwxyz"; octave:2> array=str(1:length(str))*1.0; octave:3> disp(str) abcdefghijklmnopqrstuvwxyz octave:4> disp(array) Columns 1 through 18: 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 Columns 19 through 26: 115 116 117 118 119 120 121 122 octave:5> whos -long -variables *** local user variables: Prot Name Size Bytes Class ==== ==== ==== ===== ===== rw- __nargin__ 1x1 8 double rwd array 1x26 208 double rwd str 1x26 26 char

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
Is there a way to determine

Is there a way to determine whether a field (e.g. "extract{1}.SetiMetaData.Name") exists in a data structure in Octave? Trying to access a non-existing field directly in any way results in an error. Using the exist() function works for variables, but apparently not fields.

Any suggestions?

Anders Feder
Offline
Joined: 2010-04-22
Posts: 618
To answer my own question,

To answer my own question, isfield() does this.