FastICA算法的实现过程及其python实现

ICA算法的数学原理

参考我的这篇文章：ICA算法的数学原理

算法实现过程

python实现

(1)从硬盘中读取两段参数相同的语音信号S1和S2，大小都是（1，m）
(2)人为指定一个混合矩阵A，将S1和S2进行混合得到矩阵D,D的大小是(2,m)
(3)将混合的结果存入硬盘中：mix1.wav和mix2.wav
(4)对D作用FastICA算法，得到估算的源信号Sr,Sr的大小是(2,m)
(6)将估算的源信号存入硬盘中:reconsruct_sound1.wav和reconsruct_sound2.wav

完整代码：

#FastICA by Leo
import numpy as np
import math
import random
import matplotlib.pyplot as plt



'''
函数名称：center_data
函数功能：对数据中心化：对输入矩阵的每个元素，都减去该元素所在行(每一行一共有m个元素)的均值
输入参数：X          要处理的矩阵，大小为(n,m)
返回参数：X_center   进行中心化处理之后的矩阵，大小为(n,m)
作者：Leo Ma
时间：2020.01.04
'''
def center_data(X):
    #沿着行的方向取均值，即计算n个麦克风在m个时刻中的均值，X_means的shape是(n,)
    X_means = np.mean(X,axis = 1)
    #将X_means增加一个新行，shape变为(n,1),X的每一列都与之对应相减
    return X-X_means[:, np.newaxis]
    


'''
函数名称：whiten_data
函数功能：对数据白化处理
输入参数：X          要处理的矩阵，大小为(n,m)
返回参数：Z          白化处理之后的矩阵，大小为(n,m)
         V          白化变换矩阵
作者：Leo Ma
时间：2020.01.04
'''
def whiten_data(X):
    #计算X的协方差矩阵，cov_X = E{XX^T)}
    cov_X = np.cov(X)
    #计算协方差矩阵的特征值和特征向量
    eigenValue,eigenVector = np.linalg.eig(cov_X)
    #将特征值向量对角化，变成对角阵，然后取逆
    eigenValue_inv = np.linalg.inv(np.diag(eigenValue))
    #计算白化变换矩阵V
    V = np.dot(np.sqrt(eigenValue_inv), np.transpose(eigenVector))
    #计算白化处理后得矩阵Z,Z=VX
    Z = np.dot(V,X)
    return Z,V
    
    
    
'''
函数名称：gx
函数功能：定义s的CDF，这里选择tanh(),对输入矩阵的每个元素都作用tanh()
输入参数：x          要处理的矩阵，大小为(n,m)
        alpha       常量，值域[1,2]，通常取alpha=1
返回参数：gx         tanh()的计算结果
作者：Leo Ma
时间：2020.01.04
'''
def gx(x,alpha=1):
    return np.tanh(alpha*x)
    
    
    
'''
函数名称：div_gx
函数功能：定义s的pdf，tanh()的导数是1-tanh()**2
输入参数：x          要处理的矩阵，大小为(n,m)
        alpha       常量，值域[1,2]，通常取alpha=1
返回参数：div_gx         tanh()导数的计算结果
作者：Leo Ma
时间：2020.01.04
'''    
def div_gx(x,alpha=1):
    return alpha*(1-gx(x)**2)

    

'''
函数名称：decorrelation_data
函数功能：对数据(W)进行去相关
输入参数：W                       要处理的矩阵，大小为(n,n)
返回参数：W_decorrelation         去相关之后的W
作者：Leo Ma
时间：2020.01.04
'''  
def decorrelation_data(W):
    #对WW.T进行特征值分解，D是特征值，P是特征向量
    D, P = np.linalg.eigh(np.dot(W, np.transpose(W