rcParams["figure.figsize"] = (10, 6)
First of all, we're going to load the test data and compute the information we want. The .npy file contains the following columns:
Once this information is extracted, we compute a mean mass and the $x, y$ coordinates of the center of mass using the fact that we know the distance between the left-right (215 * 2 mm) and top-down (117.5 * 2 mm) extremities of the board.
raw_data = np.load('files/WiiBoard_data.npy')
t = raw_data[:, 0]*1e-6 #t in seconds
tl = raw_data[:, 1]
tr = raw_data[:, 2]
bl = raw_data[:, 3]
br = raw_data[:, 4]
m = raw_data[:, 5]
reference_mass = m.mean()
R = tr + br
L = tl + bl
T = tr + tl
B = br + bl
x, y = (215 * (R - L) / reference_mass,
117.5 * (T - B) / reference_mass) # in mm
First, let's take a look at the mass curve.
plot(t, m, label='mass')
plot(t, ones(m.shape) * reference_mass, label='mean mass')
legend()
<matplotlib.legend.Legend at 0x8513950>
(m > reference_mass + 1).nonzero()
(array([], dtype=int32),)
(m < reference_mass - 1).nonzero()
(array([], dtype=int32),)
The mass is always comprised between + and - 1 kilograms of the average mass. So here, the measurements seem very coherent.
Let's now take a look at the displacement as a function of time.
plot(t, x, label="x")
plot(t, y, label="y")
xlabel("time (s)")
ylabel("coordinate of center of gravity (mm)")
legend()
<matplotlib.legend.Legend at 0x8e07ef0>
As we can see, the coordinates are not centered. Let's see what it looks like if we center them.
plot(t, x - x.mean(), label="x")
plot(t, y - y.mean(), label="y")
xlabel("time (s)")
ylabel("coordinate of center of gravity (mm)")
legend()
<matplotlib.legend.Legend at 0x86ccbb0>
One of the things we see is that the amplitude of $x$ and $y$ displacement is not the same. In fact, the ratio between the two is almost 4!
print x.max() - x.min()
print y.max() - y.min()
4.46394312903 17.906250965
Let's see what the trajectory looks like in the $(x, y)$ plane.
plot(x - x.mean(), y - y.mean())
xlabel("horizontal direction (mm)")
ylabel("vertical direction (mm)")
title("Center of gravity path in the x-y plane")
axis('equal')
(-3.0, 3.0, -10.0, 15.0)
Let's now work with centered data, so that we can use the classical Fourier transform tools.
x_centered = x - x.mean()
y_centered = y - y.mean()
Let's start with the $y$ axis data.
sample_freq = 1 / (t[1] - t[0])
print sample_freq
y_fft = fft.fft(y_centered)
freqs = linspace(0, 1, num=y.shape[0], endpoint=False) * sample_freq / 2
plot(freqs, abs(y_fft))
90.991810737
[<matplotlib.lines.Line2D at 0x90998b0>]
As expected, there are some low frequency components and no high frequencies. Let's focus on the low frequencies.
plot(freqs, abs(y_fft))
xlim((0, 0.5))
(0, 0.5)
We can then output the frequency content of each band of frequencies.
def print_freq_content(y_fft, freqs):
EX1 = sum(abs(y_fft[(freqs > 0.01) & (freqs <0.5)] ** 2))
EX2 = sum(abs(y_fft[(freqs > 0.5) & (freqs < 2)] ** 2))
EX3 = sum(abs(y_fft[ freqs > 2] ** 2))
print "frequency band 1 = %s" % (EX1/(EX1+EX2+EX3)*100)
print "frequency band 2 = %s" % (EX2/(EX1+EX2+EX3)*100)
print "frequency band 3 = %s" % (EX3/(EX1+EX2+EX3)*100)
print_freq_content(y_fft, freqs)
frequency band 1 = 34.4120533146 frequency band 2 = 0.484694114847 frequency band 3 = 65.1032525706
Now, what about the $x$ axis?
x_fft = fft.fft(x_centered)
plot(freqs, abs(x_fft))
[<matplotlib.lines.Line2D at 0x9e3f0b0>]
plot(freqs, abs(x_fft))
xlim((0, 0.5))
(0, 0.5)
print_freq_content(x_fft, freqs)
frequency band 1 = 41.6957576723 frequency band 2 = 1.86334135102 frequency band 3 = 56.4409009767
dx = x[1:] - x[:-1]
dy = y[1:] - y[:-1]
dx2 = dx * dx
dy2 = dy * dy
L = np.sqrt(dx2 + dy2).sum()
L
935.87185703233376
The center of mass travelled 936 mm during the 51 s of the recorded data.
dt = t[1:] - t[:-1]
vitesse_inst = np.sqrt((dx2 + dy2) / dt)
vitesse_inst.mean()
1.9325976879295166
The mean speed of the $x y$ trajectory of the center of mass was 1.93 mm/s.
points=np.vstack((x,y)).transpose()
cov = np.cov(points, rowvar=False)
def eigsorted(cov):
vals, vecs = np.linalg.eigh(cov)
order = vals.argsort()[::-1]
return vals[order], vecs[:,order]
vals, vecs = eigsorted(cov)
theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))
nstd = 2
width, height = 2 * nstd * np.sqrt(vals)
ellip = mpl.patches.Ellipse(xy=pos, width=width, height=height, angle=theta)
ellip.set_alpha(0.2)
gca().add_artist(ellip)
scatter(x, y)
axis('equal')
xlabel('x position (mm)')
ylabel('y position (mm)')
title('Covariance ellipse of point cloud during balance state')
annotate("width = %.2f, height = %.2f, theta = %.2f" % (width, height, theta), (x.min(), y.max()))
<matplotlib.text.Annotation at 0xb42ba50>