Question on signal detection here http://stats.stackexchange.com/questions/135262/statistical-criteria-for-signal-registration
import numpy as np
import matplotlib.pylab as plt
from scipy import ndimage
import pandas as pd
%matplotlib inline
data = np.genfromtxt('data/signal.csv', delimiter=',')
X1 = data[0:13]
X2 = data[13:26]
X3 = data[26:39]
X4 = data[39:52]
plt.figure(figsize=(8, 3))
plt.imshow(np.hstack([X1, X2, X3, X4]), interpolation='none')
plt.show()
th = 0.2
plt.figure(figsize=(8, 3))
plt.imshow(np.hstack([X1, X2, X3, X4]) >= th, interpolation='none')
plt.show()
We can convert these images into a $(x,y)$ data set, and then find the principal direction. To do this, we first pivot the matrix, and then filter out all values that are below the threshold
img = X4
df = pd.DataFrame(img)
stacked = pd.DataFrame(df.stack())
stacked.reset_index(inplace=True)
stacked.columns = ('y', 'x', 'v')
filtered = stacked[stacked.v >= th]
X = filtered[['x', 'y']].values
X_cent = X - X.mean(axis=0)
plt.scatter(X_cent[:, 0], X_cent[:, 1], c='white')
plt.show()
After pivoting, filtering and centering, we can do PCA and find the principal direction
_, _, V = np.linalg.svd(X_cent)
v = V[0, :]
plt.scatter(X_cent[:, 0], X_cent[:, 1], c='white')
plt.quiver(v[0], v[1], angles='xy', scale_units='xy', scale=0.2)
plt.show()
plt.subplot(1, 2, 1)
plt.title('original')
plt.imshow(img)
plt.quiver([6], [6], v[0], v[1],
angles='xy',scale_units='xy',scale=0.4)
plt.subplot(1, 2, 2)
deg = 180 * np.arccos(v[0]) / np.pi
X_rot = ndimage.rotate(img, deg, reshape=False)
plt.title('rotated by %d C' % deg)
plt.imshow(X_rot)
plt.show()
Now, when the image is rotated, we can sum $y$ out and get this
plt.plot(X_rot.sum(axis=0))
plt.show()
This is how no signal looks like:
plt.plot(X1.sum(axis=0))
plt.show()
So maybe we can subtract this from our processed data
base = X1.sum(axis=0)
plt.plot(X_rot.sum(axis=0) - base)
plt.show()
s = X_rot.sum(axis=0) - base
f = np.fft.fft(s)
fabs = np.abs(np.fft.fftshift(f))
fabs = fabs / fabs.max()
plt.plot(fabs[7:])
plt.show()
Let's see what happens with other images
for idx, img in enumerate([X1, X2, X3, X4]):
print 'data set #', idx
df = pd.DataFrame(img)
stacked = pd.DataFrame(df.stack())
stacked.reset_index(inplace=True)
stacked.columns = ('y', 'x', 'v')
filtered = stacked[stacked.v >= th]
X = filtered[['x', 'y']].values
X_cent = X - X.mean(axis=0)
_, _, V = np.linalg.svd(X_cent)
v = V[0, :]
plt.subplot(1, 2, 1)
plt.imshow(img)
plt.quiver([6], [6], v[0], v[1],
angles='xy',scale_units='xy',scale=0.4)
plt.subplot(1, 2, 2)
deg = 180 * np.arccos(v[0]) / np.pi
if abs(deg) > 1e-6:
X_rot = ndimage.rotate(img, deg, reshape=False)
else:
X_rot = img
plt.title('rotated by %d C' % deg)
plt.imshow(X_rot)
plt.show()
plt.subplot(2, 1, 1)
s = X_rot.sum(axis=0) - base
plt.title('1d signal, time')
plt.plot(s)
plt.xticks([])
plt.subplot(2, 1, 2)
f = np.fft.fft(s)
fabs = np.abs(np.fft.fftshift(f))
fabs = fabs / fabs.max()
plt.title('1d signal, frequency')
plt.plot(fabs)
plt.xticks([])
plt.show()
data set # 0
data set # 1
data set # 2
data set # 3
Looks like spectrum of the 1d signal might be a good feature. Since we don't have much data, we may use just 1st and 3rd amplitudes
def extract_features(img, th, base):
df = pd.DataFrame(img)
stacked = pd.DataFrame(df.stack())
stacked.reset_index(inplace=True)
stacked.columns = ('y', 'x', 'v')
filtered = stacked[stacked.v >= th]
X = filtered[['x', 'y']].values
X_cent = X - X.mean(axis=0)
_, _, V = np.linalg.svd(X_cent)
v = V[0, :]
deg = 180 * np.arccos(v[0]) / np.pi
if abs(deg) > 1e-6:
X_rot = ndimage.rotate(img, deg, reshape=False)
else:
X_rot = img
s = X_rot.sum(axis=0) - base
f = np.fft.fft(s)
fabs = np.abs(np.fft.fftshift(f))
fabs = fabs / fabs.max()
return fabs[[6, 8]]
f2 = extract_features(X2, th, base)
f3 = extract_features(X3, th, base)
f4 = extract_features(X4, th, base)
F = np.array([f2, f3, f4])
print F
plt.scatter(F[:, 0], F[:, 1], c=['white', 'black', 'black'])
plt.show()
[[ 0.81180072 0.45749437] [ 0.54543749 0.26524667] [ 0.47091471 0.25209055]]