We use GNU Octave for data processing, plotting, and simple programming exercises. All student work is done in an online environment using Octave Online. This allows students to access their data and files from any computer on campus or at home.

Various macros and functions are provided here to simplify the data processing tasks and for practicing Octave use.

Function index

Loading data

Click on the file name to see the listing.
csv.m: A = csv(name)

Loads a csv file specified by name.csv and shows a preview plot

Example: To load a file tilt_acc.csv into matrix A

octave:1> A = csv("tilt_acc");
csv file preview

Note that you don't need to add the .csv at the end of the file name.

function G = csv(name) . if (index(name, ".csv")) fn = name; else fn = [name ".csv"]; endif . G = csvread(fn, 3, 0); . rows = size(G,1); columns = size(G, 2); . maxrows = 500; . if ( rows < maxrows ) P = G; else warning("off"); P = G(1:rows/maxrows:rows,:); warning("on"); endif . switch(columns) case(2) plot(P(:,1),P(:,2),'r','linewidth',2); legend('X'); legend('boxoff'); xlabel('Time (s)'); title(strrep(fn,'_','\_')); set(gca(), 'fontsize', 16); set(gca(), 'linewidth', 2); case(3) plot(P(:,1),P(:,2),'r','linewidth',2,P(:,1),P(:,3),'b','linewidth',2); legend('X', 'Y'); legend('boxoff'); xlabel('Time (s)'); title(strrep(fn,'_','\_')); set(gca(), 'linewidth', 2); set(gca(), 'fontsize', 16); case(4) plot(P(:,1),P(:,2),'r','linewidth',2,P(:,1),P(:,3),'b','linewidth',2,P(:,1),P(:,4),'g','linewidth',2); legend('X', 'Y', 'Z'); legend('boxoff'); xlabel('Time (s)'); title(strrep(fn,'_','\_')); set(gca(), 'fontsize', 16); set(gca(), 'linewidth', 2); endswitch . endfunction
wav.m: [t, s] = wav(name)

Loads a wav file specified by name.wav.

Example: To load a file snd1.wav into vectors t and s

octave:1> [t, s] = wav("snd1");

Note that you don't need to add the .wav at the end of the file name.

The time values of the sound samples are stored in t and the sound wave samples are in s. function [t,W] = wav(name) . fp = fopen([name ".wav"], "r"); . H = fread(fp, 22, "uint16"); . S = H(13); N = (H(21) + H(22)*65536) / 2; printf("Sample rate: %d Hz\n", S) printf("Number of samples: %d\n", N) . W = fread(fp, N, "int16", "l"); t = linspace(0, (N-1)/S, N)'; . fclose(fp); . endfunction

disk.m: disk

Shows the amoun of disk space used in Mbytes

Example: To show current disk usage

octave:1> disk
Using 5.3 MBytes of disk space

printf("Using %.1f MBytes of disk space\n",sum([dir().bytes])/1024^2)

Selecting data

point.m: point(A, t)

Find the index of vector or matrix that corresponds to time t.

Example: Find the index of a line in matrix A that corresponds to a time value of 3.5

octave:1> A = csv("move1");
octave:2> point(A, 3.5)
ans = 233
octave:3> A(230:236,:)
ans =
       3.467000   0.509370  -0.254530   8.857050
       3.468000   1.583170  -0.735770   8.853460
       3.485000   1.834560  -1.070960   8.854660
       3.502000   1.165380  -0.392200   8.444050
       3.505000   1.530490  -0.499940   8.960000
       3.520000   1.227630  -0.304810   8.536230
       3.535000   0.990600  -0.550220   9.244910
function index = point(v, t) . for index = 1 : length(A) . if v(index) >= t return endif . endfor . endfunction
remove.m: A = remove(A, n)

Removes element n from vector A.

Example: Remove element number 3 from a vector:

octave:1> a = [8 7 6 5 4]
a =

8 7 6 5 4

octave:1> a = remove(a, 3)
a =

8 7 5 4

function R = remove(D, n); . R = [D(1:n-1) D(n+1:end)]; . endfunction

Plotting

plotx.m:
plotx(A)
plotx(A, t1)
plotx(A, t1, t2)
plotx(A, t1, t2, style)
plotx(A, style)

Plot the x-component of the sensor data in matrix A that was loaded with csv. If t1 is given, plot starts from that time point. Plot end time can be given by t2. An optional style can be given to change the plot color and line style. The style parameter can be given even when the time limits are not used.

Example: plot the x-component of data in file move1.csv

octave:1> A = csv("move1");
octave:2> plotx(A)
x component plot
octave:3> plotx(A, 2.5, 4.5)
x component plot with time limits
octave:4> plotx(A, 2.5, 4.5, 'b')
x component plot with time limits and style
octave:5> plotx(A, 'k')
x component plot with time limits and style

function plotx(A, t1=0, t2=NaN, m='r-') if (strcmp(typeinfo(t1),'string') || strcmp(typeinfo(t1),'sq_string')) m = t1; t1 = 0; endif . if (strcmp(typeinfo(t2),'string') || strcmp(typeinfo(t2),'sq_string')) m = t2; t2 = NaN; endif . n1=point(A,t1); n2=point(A,t2); . plot(A(n1:n2,1), A(n1:n2,2), m) . endfunction
ploty.m:
ploty(A)
ploty(A, t1)
ploty(A, t1, t2)
ploty(A, t1, t2, style)
ploty(A, style)

Plot the y-component of the sensor data in matrix A that was loaded with csv. The parameters are the same as for plotx.

function ploty(A, t1=0, t2=NaN, m='r-') if (strcmp(typeinfo(t1),'string') || strcmp(typeinfo(t1),'sq_string')) m = t1; t1 = 0; endif . if (strcmp(typeinfo(t2),'string') || strcmp(typeinfo(t2),'sq_string')) m = t2; t2 = NaN; endif . n1=point(A, t1); n2=point(A, t2); . plot(A(n1:n2,1), A(n1:n2,3), m) . endfunction
wavplot.m:
wavplot(t, s)
wavplot(t, s, n)

Plot the sound wave data with time in t and the sound signal in s. By default, the function plots the min/max values of data in blocks of 256 points, reducing the plot size and allowing large sound files to be viewed. The block size can be adjusted to keep the total point count below about 1000. In this case, the third parameter n sets the block size.

function wavplot(t,s,n=256) . N = (floor(length(s)/n)-1)*n; . np = 1; for i=1:n:N tt(np)=t(i); ss(np)=min(s(i:i+n-1)); np = np + 1; tt(np) = t(i+n/2); ss(np)=max(s(i:i+n-1)); np = np + 1; endfor . plot(tt,ss) . endfunction
deco.m: deco

Makes the linewidth of the last plot wider and the label text larger.

Example: decorate a plot

octave:1> A = csv("move1");
octave:2> plotx(A)
x component plot
octave:3> deco
plot line width and font size change

function deco() . set(findall('type','line'), 'linewidth', 2); set(gca(), 'fontsize', 18); set(gca(), 'linewidth', 2); . endfunction

Simulating data

damp.m: F = damp(T, tzero, offset, period, amp, decay, phase)

Calculate a damped cosine function shape as a function of time T. Parameters are: starting time tzero, background offset offset, amplitude amp, oscillation period in seconds period, the decay factor decay, and the initial phase of the cosine phase (in radians).

Example: Simulate a decaying cosine function starting at 21.3 second point with an amplitude 5.8, offset 9.6, and period of 3.4 seconds. The decay factor is 0.006, initial phase is zero:

octave:1> Fy = damp(St, 21.3, 9.6, 3.4, 5.8, .006, 0);

function F = damp(T, t0, offset, period, amp, decay, phase=0) . for i = 1 : length(T) . if ( T(i) < t0 ) F(i) = offset; else y = cos(phase + 2*pi/period*(T(i)-t0)); a = amp * exp(-decay*(T(i)-t0)); F(i) = y * a + offset; endif . endfor . endfunction

Background correction

isback.m: b = isback(V, p, min, max, n)

Determine if a point in vector V is background or not. Returns 1 if point at index p in vector V is a background point or 0 if it isn't. Parameters min and max set the minimum and maximum thresholds for separating background from motion data. The parameter n can be used to set the number of neighboring points to analyze on either side of index p.

Example: Check if points 203 and 206 in the x-component of the acceleration data represent the stationary background or not.

octave:1> A = csv("move1");
octave:2> At=A(:,1);Ax=A(:,2);Ay=A(:,3);
octave:3> isback(Ax, 203, .2, .6, 2)
ans = 1
octave:4> isback(Ax, 205, .2, .6, 2)
ans = 0

function t = isback(A, p, amin, amax, n) if ((p-n) < 1) t = 1; return endif . if ((p+n) > length(A)) t = 1; return endif . for i = p-n : p+n if (A(i) < amin) t = 0; return endif . if (A(i) > amax) t = 0; return endif endfor . t = 1; . endfunction

You can trace the operation of the isback function in the panel below. Use the step, stop, and play buttons at the bottom of the source listing and observe how the code executes and how variables change. The trace starts with the point index p=1 and is repeated for all data points in sequence. The function return value is shown for each function execution in the 'function return value' table.

isbackdelta.m: b = isbackdelta(V, p, min, max, delta, n)

The operation is identical to the isback function, except for the addition of a delta parameter. Data points are only considered to be background points if all neighbors are within the amin to amax limits and the difference between each point pair in the search neighborhood is less than delta.

The benefit of this function over isback is that it can distinguish regions where a small signal falls into the background search limits and might therefore be incorrectly assumed to be a part of the background.

Example: Calculating the background with the back function produces and incorrect background.

octave:1> A = csv("move1");
octave:2> t=A(:,1);ax=A(:,2);
octave:3> [Bt,Bv]=back(t,ax,-1,1,2);
octave:4> plot(t,ax,'r',Bt,Bv,'b')
Incorrect background extraction

Repeating the calculation with the delta parameter used in the backdelta function produces a much better result with the delta parameter set to 0.1.

octave:5> [Bt,Bv]=backdelta(t,ax,-1,1,.1,2);
octave:6> plot(t,ax,'r',Bt,Bv,'b')
Correct background extraction

function t = isbackdelta(A, p, amin, amax, delta, n) if ((p-n) < 1) t = 1; return endif . if ((p+n) > length(A)) t = 1; return endif . for i = p-n : p+n if (A(i) < amin) t = 0; return endif . if (A(i) > amax) t = 0; return endif . if (abs(A(i)-A(i-1)) > delta) t = 0; return endif endfor . t = 1; . endfunction
back.m: [Bt,Bv] = back(T, V, min, max, n)

Calculate the background data for acceleration or velocity data for a single axis in vector V. The time values for the data are given in vector T. The threshold parameters, min, max, and n are the same as for isback function.

The function returns two arrays with background time values in Bt and background values in Bv

Example: Calculate the background profile for the x-axis acceleration data and plot the result.

octave:1> A = csv("move1");
octave:2> At=A(:,1);Ax=A(:,2);Ay=A(:,3);
octave:3> [Baxt, Bax] = back(At, Ax, .2, .6, 2);
octave:4> plotx(A)
octave:5> hold
octave:6> plot(Baxt, Bax, 'b')
octave:7> hold
octave:8> deco
calculated background profile

function [Btime, Bvalue] = back(T, A, amin, amax, n) Bn = 1; Bt = -1; . for i = 1 : length(A) . if (T(i) > Bt && isback(A, i, amin, amax, n)) Btime(Bn) = T(i); Bvalue(Bn) = A(i); Bn = Bn + 1; Bt = T(i); endif . endfor . endfunction

You can trace the operation of the back function in the panel below. Use the step, stop, and play buttons at the bottom of the source listing and observe how the code executes and how variables change.

backdelta.m: [Bt,Bv] = backdelta(T, V, min, max, delta, n)

Calculate the background data for acceleration or velocity data for a single axis in vector V. The operation is similar to the back function except that a delta parameter has been added and the function uses the isbackdelta function to extract background points.

The function returns two arrays with background time values in Bt and background values in Bv

Example: Calculate the background profile for the x-axis acceleration data and plot the result.

octave:1> A = csv("move1");
octave:2> t=A(:,1);ax=A(:,2);
octave:3> [Bt,Bv]=backdelta(t,ax,-1,1,2);
octave:4> plot(t,ax,'r',Bt,Bv,'b')
Correct background extraction

function [Btime, Bvalue] = backdelta(T, A, amin, amax, da, n) Bn = 1; Bt = -1; . for i = 1 : length(A) . if (T(i) > Bt && isbackdelta(A, i, amin, amax, da, n)) Btime(Bn) = T(i); Bvalue(Bn) = A(i); Bn = Bn + 1; Bt = T(i); endif . endfor . endfunction
backave.m: [Bt,Bv] = backave(T, V, min, max, delta, n)

Calculate the background data for acceleration or velocity data for a single axis in vector V. The operation is similar to the back function but instead of copying data points to the bacground vector, the function finds each stop period in the data and returns just two data points for each stop section with the values set to the background average value for each stop section. The functon parameters are the same as those for backdelta.

The function returns two arrays with background time values in Bt and background values in Bv

Example: Calculate the background profile for the x-axis acceleration data and plot the result.

octave:1> A = csv("move1");
octave:2> t=A(:,1);ax=A(:,2);
octave:3> [Bt,Bv]=backave(t,ax,-1,1,.1,2);
octave:4> plot(t,ax,'r',Bt,Bv,'b.-','markersize',8)
octave:5> axis([0 10 .2 .8])
Average background extraction

function [Btime, Bvalue] = backave(T, A, amin, amax, delta, n) Bn = 1; Bfound = 0; . for i = 1 : length(A) . if (isbackdelta(A, i, amin, amax, delta, n)) if (Bfound == 0) Bfound = i; endif else if (Bfound > 0) Bave = mean(A(Bfound:i)); Btime(Bn) = T(Bfound); Bvalue(Bn) = Bave; Btime(Bn+1) = T(i); Bvalue(Bn+1) = Bave; Bn = Bn + 2; Bfound = 0; endif endif endfor . if (Bfound > 0) Bave = mean(A(Bfound:i)); Btime(Bn) = T(Bfound); Bvalue(Bn) = Bave; Btime(Bn+1) = T(i); Bvalue(Bn+1) = Bave; endif . endfunction
fixback.m: F = fixback(T, V, Bt, Bv)

Returns the background-corrected acceleration or velocity data F. The time vector is given in T with the corresponding acceleration or velocity data in V. The background data calculated with the back function should be in Bt and Bv.

Example: Calculate the background-corrected x-axis acceleration

octave:1> A = csv("move1");
octave:2> At=A(:,1);Ax=A(:,2);Ay=A(:,3);
octave:3> [Baxt, Bax] = back(At, Ax, .2, .6, 2);
octave:4> Fx = fixback(At, Ax, Baxt, Bax);
octave:5> plot(At, Fx);deco
background-corrected x component plot

function F = fixback(T, A, Bt, Bv) . F = A - interp1(Bt, Bv, T); . endfunction

Signal processing

integrate.m: V = integrate(T, F)

Integrates the acceleration or velocity data in F, returning either velocity data (integrating acceleration) or position data (integrating velocity). The time values should be in vector T.

Example: Calculate the x-axis velocity from the background-corrected acceleration data

octave:1> A = csv("move1");
octave:2> At=A(:,1);Ax=A(:,2);Ay=A(:,3);
octave:3> B = back(At, Ax, .2, .6, 2);
octave:4> Fx = fixback(At, Ax, B);
octave:5> Vx = integrate(At, Fx);
octave:6> plot(At, Vx);deco
x-axis velocity plot

function V = integrate(t, A) V(1) = 0; . for i = 2 : length(A) dt = t(i) - t(i-1); dv = dt * A(i); V(i) = V(i-1) + dv; endfor . V = V'; . endfunction
ave.m:
a = ave(V)
a = ave(T, V, t1, t2)

Calculates the average of a vector V. The second form of the function calculates an average over a time interval from t1 to t2. In this case, the first argument is the time vector T.

Example: Calculate the x-axis acceleration average over a time interval from 3 s to 4 s.

octave:1> A = csv("move1");
octave:2> At=A(:,1);Ax=A(:,2);Ay=A(:,3);
octave:3> aa = ave(At, Ax, 3, 4);
aa = 0.51478

function value = ave(T, V=T, t1=0, t2=NaN) . if ( T == V ) value = sum(T) / length(T); return endif . n1 = point(T, t1); n2 = point(T, t2); sum = 0; . for i = n1 : n2 sum = sum + V(i); endfor . value = sum / (n2 - n1 + 1); . endfunction
avex.m:
a = avex(A)
a = avex(A, t1, t2)

Calculates the average of the x-axis component of the data in matrix A. The second form of the function calculates an average over a time interval from t1 to t2.

Example: Calculate the x-axis acceleration average over a time interval from 3 s to 4 s.

octave:1> A = csv("move1");
octave:2> aa = avex(A, 3, 4);
aa = 0.51478

function value = avex(A, t1=0, t2=NaN) . n1 = point(A, t1); n2 = point(A, t2); sum = 0; . for row = n1 : n2 sum = sum + A(row,2); endfor . value = sum / (n2 - n1 + 1); . endfunction
mavg.m:
a = mavg(v, w)

Calculates a moving average of the vector v with a window of +/- w points.

Example: Calculate the moving average of noisy accelereometer data with a window of +/- 5 data points.

octave:1> N = csv("noise_acc");
octave:2> t = N(753:end,1);
octave:3> z = N(753:end,4);
octave:4> plotx(t,z);
z-axis noisy acceleration
octave:5> a = mavg(z, 5);
octave:6> plotx(t,a);
z-axis smoothed acceleration

Note that this type of averaging changes the original data. It can be a useful technique for handling noisy data, but care should be taken that the averaging process doesn't create systematic errors in final calculation results.

function A = mavg(D, w=1) . for i = 1 : length(D) n1 = max(1, i - w); n2 = min(length(D), i + w); s = 0; for j = n1:n2 s = s + D(j); endfor A(i) = s / (n2 - n1 + 1); endfor . endfunction
boxave.m:
a = boxave(T, M)
a = boxave(T, M, points, limit)

Calculates the average over blocks of data M. The default block length is 50 data points, but can be set with the optional points parameter. If data variation (max - min difference) is less than a critical value, the block is considered to contain stable data. The default limit is 1, but can be set with the optional limit parameter.

the function returns two vectors, the time and field values for all stable regions found in the data.

Example: Load mgnetic field mapping data and find all stable points in the x-axis data:

octave:1> M = csv("map12");
octave:2> T=M(:,1);Mx=M(:,2);My=M(:,3);Mz=M(:,4);
octave:3> [Xt, Xv] = boxave(T, Mx);
octave:4> plot(T,Mx,'b');hold
octave:5> plot(Xt,Xv,'r.');hold

function [At, Av] = boxave(T, V, points=50, limit=1) . An = 1; Rold = limit; . for i = points : points: length(V) . R = range(V(i-points+1 : i)); if ( R < limit ) if ( R < Rold ) At(An) = T(fix(i-points/2)); Av(An) = mean(V(i-points+1 : i)); Rold = R; endif else if (Rold < limit) An = An + 1; endif Rold = limit; endif . endfor . endfunction

Signal conditioning

difference.m: Dx = difference(Bx)

Calculates the point-to-point differences of the vector Bx. Dx[0] is always zero. For other elements Dx[i] = Bx[i] - Bx[i-1]

Example: Remove slow background variations from time-series data:

octave:1> B = csv("mag");
octave:2> Bt=A(:,1);Bx=A(:,2);
octave:3> plot(Bt,Bx);deco
Raw time series data
octave:4> Bxd = difference(Bx);
octave:5> plot(Bt, Bxd);deco
Difference data

function D = difference(V) . sz = size(V); if (sz(2) > sz(1)) D = [0; diff(V)']'; else D = [0 diff(V)']'; endif . endfunction
limits.m: L = limits(V, vmin, vmax)

Compares data points to minimum and maximum thresholds. Points below the minimum threshold vmin are set to -1, points above the maximum threshold vmax are set to +1 and points between thresholds are set to 0.

Example: Threshold a time sequence for pulse feature extraction:

octave:1> B = csv("mag");
octave:2> Bt=A(:,1);Bx=A(:,2);
octave:3> Bxd = difference(Bx);
octave:4> plot(Bt, Bxd);deco
Difference data
octave:5> Lx = limits(Bxd, -1.5, 1.5);
octave:6> plot(Bt, Lx);axis([0 25 -1.2 1.2]);deco
Thresholded data

function L = limits(V, vmin, vmax) . for i = 1 : length(V) L(i) = 0; if (V(i) < vmin) L(i) = -1; elseif (V(i) > vmax) L(i) = 1; endif endfor . endfunction
grow.m: L = grow(V)

Increases the +1 and -1 domain sizes by one data point in in either direction for data sets that have been thresholded to -1/0/+1 domains. Repeated application of the grow function can be used to eliminate noise from a data set.

Example: Expand the detected pulse features in thresholded data:

octave:1> B = csv("mag");
octave:2> Bt=A(:,1);Bx=A(:,2);
octave:3> Bxd = difference(Bx);
octave:4> Lx = limits(Bxd, -1.5, 1.5);
octave:5> plot(Bt, Lx);axis([10 15 -1.5 1.5])
Thresholded data
octave:6> Lxg = grow(Lx);
octave:7> plot(Bt, Lxg);axis([10 15 -1.5 1.5])
Thresholded data
octave:8> Lxg = grow(Lxg);
octave:9> plot(Bt, Lxg);axis([10 15 -1.5 1.5])
Thresholded data
octave:10> Lxg = grow(Lxg);
octave:11> plot(Bt, Lxg);axis([10 15 -1.5 1.5])
Thresholded data

function G = grow(V) . G = V; . for i = 2 : length(V)-1 if (V(i) != 0) if (V(i+1) == 0) G(i+1) = V(i); elseif (V(i-1) == 0) G(i-1) = V(i); endif endif endfor . endfunction
envelope.m:
E = envelope(V)
E = envelope(V, factor)
E = envelope(V, factor, flat)

Calculates the envelope profile of a positive signal V. After every local peak in V, the decay of the envelope towards the baseline occurs at a rate defined by factor. The value of factor should be between 0 and 1. Values closer to 1 result in a slower return to baseline. A reasonable default value is 0.9999. It is possible to enforce a constant region after each local maximum during which there is no envelope decay. The length of the constant region is given by the optional parameter flat that defines the length as a number of data points in V. For microphone data, values of 100 to 1000 may be reasonable.

Example: Calculate the envelope profile for a dart shot sound:

octave:1> D = csv("dart1");
octave:2> Dt=D(:,1);Ds=D(:,2);
octave:3> Da = abs(Ds);
octave:4> plot(Dt, Da);
Microphone data
octave:5> De = envelope(Da, 0.999);
octave:6> plot(Dt, De);deco
Fast envelope

The decay in this example is set by factor 0.999, which is quite fast. For a slower factor of 0.9995, we get:

octave:7> De = envelope(Da, 0.9995);
octave:8> plot(Dt, De);deco
Slower envelope

function E = envelope(S, f=.9999, d=0) . v = 0; fn = 1; . for i = 1 : length(S) if ( v < S(i) ) E(i) = S(i); v = S(i); fn = 1; else if ( fn > d ) v = v * f; endif fn = fn + 1; E(i) = v; endif endfor . endfunction

Finding patterns

message.m:
Text = message(L)
Text = message(L, charspace, bitspace)

Finds a text message from a thresholded -1/0/1 bit pattern L. The optional charspace parameter sets the number of zero points that should mark a space between characters. Default for charspace is 50 zeros. The bitspace sets the number of zeros that mark a space between bits in a single character. Deafult is 10.

function Txt = message(L, Charspace=50, Bitspace=10) . # Number of zeros found in series Nzero = 0; # Number of 1 to -1 or -1 to 1 changes found Nflips = 0; # Bit pattern found so far for one character Bits = ""; # Extracted text Txt = " "; . for i = 2 : length(L) . # Count +1/-1 sign changes if (L(i-1) * L(i) < 0) Nflips = Nflips + 1; endif . if (L(i) != 0) Nzero = 0; else # Count zeros Nzero = Nzero + 1; . # Detect the start of a new character if (Nzero > Charspace) if (length(Bits)) Txt = [Txt char(bin2dec(Bits))]; disp(Txt); endif Nzero = 0; Nflips = 0; Bits = ""; endif . # Detect the end of a single bit if (Nflips != 0 && Nzero > Bitspace) if (Nflips < 3) disp ('1'); Bits = [Bits "1"]; elseif (Nflips >= 3) disp ('0'); Bits = [Bits "0"]; else disp ('Error at') disp (i) disp (Nflips) endif Nflips = 0; endif . endif . endfor . endfunction
radiusx.m:
R = radiusx(A, G)
[R, Rw] = radiusx(A, G)
R = radiusx(A, G, start, end)
[R, Rw] = radiusx(A, G, start, end)

Takes the acceleration matrix A and rotation rate G (loaded with the csv function) and calculates the average x-axis component of the rotation radius vector. The averaging can be limited to a certain time section of the data by specifying the desired start and end time limits. The function returns the average radius for the given time interval and optionally the point-by-point radius data that can be used for plotting.

Example: Calculate the average rx for a given set of acceleration data (a30.csv) and gyroscope data (g30.csv):

octave:1> A = csv("a30");
octave:2> G = csv("g30");
octave:3> plotx(A);
x-axis acceleration

The background level is calculated from the first 1-second interval and the radius value can be averaged over a time interval from 2 to 10 seconds:

octave:4> R = radiusx(A, G, 2, 10)
R = 0.41940

To view the calculated radius data, use the optional return value Rw:

octave:5> [R, Rw] = radiusx(A, G, 2, 10);
octave:6> plot(A(:,1), Rw)
x-axis radius

This data contains divide-by-zero points where the rotation rate is nearly zero. Set reasonable axis values to see to interesting data and add the average value plot:

octave:7> axis([0 15 -1 1]);
octave:8> hold
octave:9> plot([2 10],[R R],'r-')
octave:8> hold
x-axis radius

function [R, Rw] = radiusx(A, G, t1=0, t2=NaN) . n1 = lookup(A(:,1), t1); n2 = lookup(A(:,1), t2); . nz = lookup(A(:,1), 1); Back = mean(A(1:nz,2)); . Rw = (Back - A(:,2)) ./ G(:,4) .^ 2; R = mean(Rw(n1 : n2)); . endfunction
radiusy.m:
R = radiusy(A, G)
[R, Rw] = radiusy(A, G)
R = radiusy(A, G, start, end)
[R, Rw] = radiusy(A, G, start, end)

Usage is identical to radiusx

function [R, Rw] = radiusy(A, G, t1=0, t2=NaN) . n1 = lookup(A(:,1), t1); n2 = lookup(A(:,1), t2); . nz = lookup(A(:,1), 1); Back = mean(A(1:nz,3)); . Rw = (Back - A(:,3)) ./ G(:,4) .^ 2; R = mean(Rw(n1 : n2)); . endfunction
accposx.m: source("accposx.m");

This is a script that can be executed by clicking the run button in the Octave-online editor window or by giving the source command as shown. The script is only valid for the demonstration data and should be adapted for actual measured data (file names and averaging time intervals).

octave:1> source("accposx.m");
x-axis linear fit

The x-axis radius values are calculated for 10 different measurements taken at different tablet positions on the turntable. The data are plotted as a function of the measured position of the tablet mounting stage and fitted with a linear function. The zero crossing (where the x-axis acceleration is zero) occurs at a position of -111.27 mm, which means that the accelerometer sensor is about 111.3 mm from the left edge of the tablet mounting stage.

A = csv("a30");G = csv("g30"); Rm(1) = 300; Rx(1) = radiusx(A, G, 2, 10); . A = csv("a25");G = csv("g25"); Rm(2) = 250; Rx(2) = radiusx(A, G, 4, 12); . A = csv("a20");G = csv("g20"); Rm(3) = 200; Rx(3) = radiusx(A, G, 4, 12); . A = csv("a15");G = csv("g15"); Rm(4) = 150; Rx(4) = radiusx(A, G, 5, 12); . A = csv("a10");G = csv("g10"); Rm(5) = 100; Rx(5) = radiusx(A, G, 4, 12); . A = csv("a05");G = csv("g05"); Rm(6) = 50; Rx(6) = radiusx(A, G, 3, 11); . A = csv("a00");G = csv("g00"); Rm(7) = 0; Rx(7) = radiusx(A, G, 4, 10); . A = csv("a-05");G = csv("g-05"); Rm(8) = -50; Rx(8) = radiusx(A, G, 3, 10); . A = csv("a-10");G = csv("g-10"); Rm(9) = -100; Rx(9) = radiusx(A, G, 3, 8); . A = csv("a-15");G = csv("g-15"); Rm(10) = -150; Rx(10) = radiusx(A, G, 3, 10); . Fx = polyfit(Rm, Rx, 1); plot(Rm, Rx, 'b.') hold plot(Rm,Fx(1)*Rm+Fx(2),'r-') hold legend('location', 'northwest') legend('',['x_0 = ' num2str(-Fx(2)/Fx(1))]) legend('boxoff')
edges.m: E = edges(T, D)

Takes as arguments a time vector T and thresholded data D. The thresholding can be done with the limits function. A list of positive edge positions (time values) is returned, where 0 to 1 transitions occur in D.

Example: Assume that we have the time vector T and the thresholded data in Dl:

octave:1> plot(T,Dl);axis([6 11 -.5 1.5])
thresholded data
octave:2> E = edges(T, Dl)
E =

Columns 1 through 8:

6.3870 6.8090 7.3750 7.8700 8.3650 8.9300 9.3570 9.8490

Column 9:

10.4150

Usually you would want time differences, which you can get with diff

octave:3> D = diff(E)
D =

0.42200 0.56600 0.49500 0.49500 0.56500 0.42700 0.49200 0.56600

function E = edges(T, D) . n = 1; E = []; . for i = 2 : length(T) if (D(i-1) == 0 && D(i) == 1) E(n) = T(i); n = n + 1; endif endfor . endfunction

DOS scripts

csvcut.js:

A script for cutting a section out of a larger CSV file in a command window on a Windows PC. Since Octave-online cannot handle very large files, it may be necessary to cut out only the necessary sections for analysis.

Example: Cut a section extending from 2.8 to 3.5 seconds in a shot1.csv file and save the section in dart1.csv

cscript //Nologo //E:jscript csvcut.js shot1.csv 2.8 3.5 > dart1.csv

// Cut a range of time values from a CSV file . t1 = parseFloat(WScript.Arguments(1)); t2 = parseFloat(WScript.Arguments(2)); . fs = new ActiveXObject("Scripting.FileSystemObject"); f = fs.GetFile(WScript.Arguments(0)); is = f.OpenAsTextStream( 1, 0 ); . // Copy header WScript.Echo(is.ReadLine()); . // skip lines until starting time do { ln = is.ReadLine(); lw = ln.split(","); t = parseFloat(lw[0]); } while( !is.AtEndOfStream && t < t1); . // copy lines until end time while( !is.AtEndOfStream && t <<= t2){ WScript.Echo(ln); ln = is.ReadLine(); lw = ln.split(","); t = parseFloat(lw[0]); } . is.Close();