[PYTHON] Implementation of independent component analysis

Independent component analysis is used for blind separation. In this method, two kinds of sounds are discriminated from a sound in which two different sounds are mixed. This time, we will implement it using one of the methods, projection pursuit. Projection pursuit is the extraction of more non-Gaussian parts from two elements. This is equivalent to extracting the most independent sound from the central limit theorem and finding non-Gauss.

Assuming that there are two devices that first hear the sound. The volume and order of the sounds do not matter.

motivation

Task

algorithm

Projection pursuit Newton approximation

Random number generation

import numpy as np
from scipy.linalg import sqrtm
import math
import matplotlib.pyplot as plt


def generate_data(n=1000):
    global M, b
    x = np.concatenate([np.random.rand(n, 1), np.random.randn(n, 1)], axis=1)
    x[0, 1] = 6   # outlier
    x = (x - np.mean(x, axis=0)) / np.std(x, axis=0)  # Standardization
    M = np.array([[1, 3], [5, 3]])
    x = x.dot(M.T)
    x = np.linalg.inv(sqrtm(np.cov(x, rowvar=False))).dot(x.T).T
    b  = np.random.uniform(0.08, -0.08, size= (2, 1))
    return x

X = generate_data()

for i in range(X.shape[0]):
  plt.scatter(generate_data()[i][0], generate_data()[i][1], color= 'red')
plt.show()

The audio is already mixed in the matrix M shown above.

Definition of function to be introduced

def g_s_3_dif(s):
  return 3 * (s**2)

def g_tan_dif (s):
  return 1 - (math.tanh(s))**2

def g_s_3 (s):
  return s**3

def g_tan (s):
  return math.tanh(s)

Spheroidization

#Is it spherical of x?
def kyujouka ():
  X_= []
  for i in range(X.shape[0]):
    x0, x1 = 0, 0
    n =X.shape[0]
    for j in range(X.shape[0]):
      x0 += X[j][1]
      x1 += X[j][1]
    x_ = (X[i] - [x0/n, x1/n])
    X_.append(x_.tolist())
  X_ = np.array(X_)
#   X_.append(x_)
  global X__
  X__ = []
  for l in range(X_.shape[0]):
    a = np.matmul(X_[l], X_[l].T)/n
    aa = 1/math.sqrt(a)
    x__ = aa*X_[l]
    X__.append(x__)
  return X__
X__ = kyujouka()
print((X__))

When using kurtosis

difference = 1
b  = np.random.uniform(0.08, -0.08, size= (1, 2))
print(b)
X__ = np.array(kyujouka())
while difference > 0.1:
  n = X.shape[0]
  sum_dif, sum_ = 0, 0
  for i in range(X.shape[0]):
    a = np.matmul(b, X__[i].T)
    g = g_s_3(a[0])
    gdif = g_s_3_dif(a[0])
    sum_dif += gdif
    sum_ += X__[i]*g
  b_before = b
  b = (sum_dif/n)*b - (sum_/n)
  b = b/math.sqrt(np.matmul(b, b.T)[0][0])
  difference = abs((b[0][0]) - abs(b_before[0][0]))
y =[]
for i in range(X.shape[0]):
  a = np.matmul(b, X[i].T)
  y.append(a[0])

aa = np.histogram(y, bins = 50)
a_bins = aa[1]
a_hist = aa[0]
X1 = []
for i in range(1, len(a_bins)):
    X1.append((a_bins[i-1]+a_bins[i])/2)
plt.bar(X1,a_hist, width=0.08)

Close to Gauss, projection pursuit is not working well.

When using tanh

difference = 1
b  = np.random.uniform(0.08, -0.08, size= (1, 2))
print(b)
X__ = np.array(kyujouka())
while difference > 0.1:
  n = X.shape[0]
  sum_dif, sum_ = 0, 0
  for i in range(X.shape[0]):
    a = np.matmul(b, X__[i].T)
    g = g_tan(a[0])
    gdif = g_tan_dif(a[0])
    sum_dif += gdif
    sum_ += X__[i]*g
  b_before = b
  b = (sum_dif/n)*b - (sum_/n)
  print(np.matmul(b, b.T)[0][0], b)
  b = b/math.sqrt(np.matmul(b, b.T)[0][0])
  print(b)
  difference = abs((b[0][0]) - abs(b_before[0][0]))

  print(difference)

y =[]
for i in range(X.shape[0]):
  a = np.matmul(b, X[i].T)
  y.append(a[0])
aa = np.histogram(y, bins =50)
print(aa)
a_bins = aa[1]
a_hist = aa[0]
X1 = []
for i in range(1, len(a_bins)):
    X1.append((a_bins[i-1]+a_bins[i])/2)
plt.bar(X1,a_hist, width=0.05)

Well separated