import sys
import time
import math
import datetime
import threading
import numpy as np
import matplotlib.pyplot as plt
def show_data(x,y):
num_list_x = np.arange(len(x))
plt.figure(figsize=(20, 10))
plt.xlim(0, len(num_list_x)+50)
plt.plot(num_list_x, x, color='b',linestyle='-', label='x')
plt.legend()
plt.show()
if all(y)!= None:
plt.figure(figsize=(20, 10))
num_list_y = np.arange(len(y))
plt.xlim(0, len(num_list_y)+50)
plt.plot(num_list_y, y, color='r',linestyle='-', label='y')
plt.legend()
plt.show()
def logistic(t):
return 1.0/(1+np.exp(-t))
def n_multiply(t, x):
return t * x
def trans_inverse(x):
if not isinstance(x, np.ndarray):
x = np.array(x)
return np.linalg.inv(x.transpose())
def ica(x):
m = len(x) # 样本数目,这里是 2 个
n = len(x[0]) # mic 数目,这里是 1000 个
w = [[0.0]*n for t in range(n)] # 建立个 1000*1000 的列表
iw = [[0.0]*n for t in range(n)]# 建立个 1000*1000 的列表
w1 = [[0.0]*n for t in range(n)]# 建立个 1000*1000 的列表
for i in range(n): # 将对角线初始化为1
w[i][i] = 1
alpha = 0.001 # 初始化学习率为
for time in range(200): # 迭代次数最多不超过 200 次
for i in range(m): # 分离出"样本数目"个样本
for j in range(n):
t = np.dot(w[j], x[i])
t = 1 - 2*logistic(t)
w1[j] = n_multiply(t,x[i]) # w1[j] = t*x[i]
iw = trans_inverse(w) # iw = w^T^(-1)
iw = w1 + iw
w1 = np.add(w1, iw) # w1 += iw
w1 = np.dot(w1, alpha)
w = np.add(w, w1)
print(time, ":\t", w)
return w
def mix(A, x,y):
mix = np.dot(A, np.array([x, y]))
for i in range(len(mix[0])):
mix[0][i] = mix[0][i] + mix[1][i]
show_data(mix[0], None)
return mix
def decode(w,x):
ps = np.dot(x, w)
return ps[0], ps[1]
if __name__ =="__main__":
s1 =[math.sin(float(x)/20) for x in range(0, 1000, 1)]
s2 = [float(x)/50 for x in range(0, 50, 1)]* 20
#s1 = [math.sin(float(x)/20) for x inrange(0, 100, 1)]
#s2 = [float(x)/50 for x in range(0, 50,1)] * 2
show_data(s1, s2)
#################
A = [[0.6, 0.4], [0.45, 0.55]] # 混合矩阵
x = mix(A, s1, s2) #s1/s2线性加权得到输入数据x
x_mean = np.mean(x) # 去均值
for i in range(x.shape[0]):
for j in range(x.shape[1]):
x[i][j] = x[i][j] - x_mean
# ica
w = ica(x)
# 解混
[ps1, ps2] = decode(w, x)
show_data(ps1, ps2)