First-year seminar
First-year seminar
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.
Click on the file name to see the listing.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");
▷
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, 0, 0);
.
while (size(G,2)>1 && G(1,1)==0)
G=G(2:end,:);
endwhile
.
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
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
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)
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
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(A, t)
.
for index = 1 : length(A)
.
if A(index) >= t
return
endif
.
endfor
.
endfunction
Removes element n from vector A.
Example: Remove element number 3 from a vector:
function R = remove(D, n);
.
R = [D(1:n-1) D(n+1:end)];
.
endfunction
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)
octave:3> plotx(A, 2.5, 4.5)
octave:4> plotx(A, 2.5, 4.5, 'b')
octave:5> plotx(A, 'k')
▷
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
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
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)
.
set(findall('type','line'), 'linewidth', 2);
set(gca(), 'fontsize', 18);
set(gca(), 'linewidth', 2);
.
endfunction
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)
octave:3> deco
▷
function deco()
.
set(findall('type','line'), 'linewidth', 2);
set(gca(), 'fontsize', 18);
set(gca(), 'linewidth', 2);
.
endfunction
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
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.
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')
▷
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')
▷
function t = isbackdelta(A, p, amin, amax, delta, n)
.
if ((p-n) < 2)
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
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
▷
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.
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')
▷
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
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])
▷
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
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
▷
function F = fixback(T, A, Bt, Bv)
.
F = A - interp1(Bt, Bv, T);
.
endfunction
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
▷
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
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
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
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);
octave:5> a = mavg(z, 5);
octave:6> plotx(t,a);
▷
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
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
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
octave:4> Bxd = difference(Bx);
octave:5> plot(Bt, Bxd);deco
▷
function D = difference(V)
.
sz = size(V);
if (sz(2) > sz(1))
D = [0; diff(V)']';
else
D = [0 diff(V)']';
endif
.
endfunction
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
octave:5> Lx = limits(Bxd, -1.5, 1.5);
octave:6> plot(Bt, Lx);axis([0 25 -1.2 1.2]);deco
▷
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
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])
octave:6> Lxg = grow(Lx);
octave:7> plot(Bt, Lxg);axis([10 15 -1.5 1.5])
octave:8> Lxg = grow(Lxg);
octave:9> plot(Bt, Lxg);axis([10 15 -1.5 1.5])
octave:10> Lxg = grow(Lxg);
octave:11> plot(Bt, Lxg);axis([10 15 -1.5 1.5])
▷
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
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);
octave:5> De = envelope(Da, 0.999);
octave:6> plot(Dt, De);deco
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
▷
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
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
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);
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)
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
▷
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
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
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");
▷
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')
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])
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
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();