Tue, 2010-06-29 14:11

Tue, 2010-06-29 15:26

#1
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)

Tue, 2010-06-29 16:04

#2
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.

Tue, 2010-06-29 16:23

#3
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.

Tue, 2010-06-29 18:08

#4
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).

Thu, 2010-07-08 01:44

#5
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?

Thu, 2010-07-08 05:01

#6
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.

Thu, 2010-07-08 05:38

#7
Thanks!

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

Thu, 2010-07-08 05:47

#8
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.

Thu, 2010-07-08 06:03

#9
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.

Sat, 2010-08-14 17:24

#10
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!

Mon, 2010-07-12 23:33

#11
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).

Sat, 2010-08-14 17:29

#12
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.

Tue, 2010-07-20 04:32

#13
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.

Sun, 2010-07-25 00:38

#14
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?

Sun, 2010-07-25 01:46

#15
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.

Sun, 2010-07-25 03:31

#16
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?

Sun, 2010-07-25 05:46

#17
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 `

Sat, 2010-08-14 04:06

#18
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?