数字信号处理之-预处理-数据白化

数据白化

1、白化的原理

随机向量的“零均值化”和“空间解相关”（也叫白化）是最常用的两个预处理过程，其中“零均值化”比较简单，而“空间解相关”涉及一些矩阵的知识。

设有均值为零的随机信号向量 X ，其自相关矩阵为
$$R_x=E[x{x}^T]\ne \text{I}$$

很明显，$R_x$是对称矩阵，且是非负定的（所有特征值都大于或等于0）。

现在，寻找一个线性变换B对X进行变换，即$y=Bx$，使得
$$R_y=B\cdot E[x{x}^T]\cdot B^T=\text{I}$$
上式的含义是：y的各分量是不相关的，即$E[Y_i{Y}_i]={\delta }_{ij}$ 。通常将这个过程称为“空间解相关”、“空间白化”或“球化”。 B称为空间解相关矩阵（空间白化矩阵、球化矩阵）。

由$R_x$的性质可知，其存在特征值分解：
$$R_x=Q\sum Q^T$$
Q是正交矩阵，$\sum$是对角矩阵，其对角元素是$R_x$特征值。

令$B={\sum }^{-1/2}{Q}^T$,则有
$$R_y=(\sum ^{-1/2}{Q}^T)Q\sum Q^T(\sum ^{-1/2}{Q}^T)T=\text{I}$$
因此，通过矩阵B线性变换后，y的各个分量变得不相关了。

对于$R_x$来说，特征值分解和奇异值分解是等价的，而奇异值分解的数值算法比特征值分解的数值算法具有更好的稳定性，因此一般都用奇异值分解来构造白化矩阵B。

术语“白”来自于白噪声的能谱在所有频率上是一个常数这一事实，就像含有各种颜色的白光谱一样。白化的本质就是去相关加缩放。从上述原理可以看出，白化的过程跟PCA是一样的，经过PCA处理过之后的数据是不相关的，且具有压缩数据的作用，但只能保证各分量信号之间“不相关”，不能保证的“独立性”，但它能够简化盲分离算法或改善分离算法的性能。

此外，白化矩阵B肯定不是唯一的。容易看到，任何矩阵UB（ U为正交矩阵）也是白化矩阵。这是因为对 $y=UBx$，下式成立：
$$E[y{y}^T]=UB\cdot E[x{x}^T]\cdot B^T{U}^T=UI{U}^T=\text{I}$$
一个重要的例子是矩阵$Q{\sum }^{-1/2}{Q}^T$ 。这也是一个白化矩阵，因为它是用正交矩阵Q左乘式$B={\sum }^{-1/2}{Q}^T$的B得到的。这个矩阵称为$Cx$的逆均方根，并用$C{x}^{-1/2}$表示，因为它来自于均方根概念向矩阵的标准推广。

2、图像数据的白化代码（matlab）：

该pca函数接口形式为：
$$[Y,V,E,D]=pca(X)$$
其中，X为输入数据，X的每一列是一个输入样本。返回值Y是对X进行PCA分析后的投影矩阵。V是与X有关的协方差矩阵特征向量的白化矩阵，E是对应的特征向量（列）构成的矩阵，D是对应的特征值构成的对角矩阵（特征值处于对角线上）。返回值中的白化矩阵，特征向量和特征值都是按照对应特征值大小进行排序后了的。

function [Y,V,E,D] = pca(X)
    % do PCA on image patches
    % INPUT variables:
    % X         matrix with image patches as columns
    % OUTPUT variables:
    % Y         the project matrix of the input data X without whiting
    % V         whitening matrix
    % E         principal component transformation (orthogonal)
    % D         variances of the principal components

    % 去除直流成分
    X = X-ones(size(X,1),1)*mean(X);
    % 其中r=size(A,1)该语句返回的时矩阵A的行数，c=size(A,2) 该语句返回的时矩阵A的列%数

    % Calculate the eigenvalues and eigenvectors of the new covariance matrix.
    covarianceMatrix = X*X'/size(X,2); % 求出其协方差矩阵

    % E是特征向量构成，它的每一列是特征向量，D是特征值构成的对角矩阵
    % 这些特征值和特征向量都没有经过排序
    [E, D] = eig(covarianceMatrix); 

    % Sort the eigenvalues and recompute matrices
    % 因为sort函数是升序排列，而需要的是降序排列，所以先取负号,diag(a)是取出a的对角元素构成
    % 一个列向量，这里的dummy是降序排列后的向量，order是其排列顺序
    [dummy,order] = sort(diag(-D));
    E = E(:,order);%将特征向量按照特征值大小进行降序排列，每一列是一个特征向量
    Y = E'*X;
    d = diag(D); %d是一个列向量

    % dsqrtinv是列向量，特征值开根号后取倒,仍然是与特征值有关的列向量
    % 其实就是求开根号后的逆矩阵
    dsqrtinv = real(d.^(-0.5));
    Dsqrtinv = diag(dsqrtinv(order));% 是一个对角矩阵，矩阵中的元素时按降序排列好了的特征值（经过取根号倒后）
    D = diag(d(order));				 % D是一个对角矩阵，其对角元素由特征值从大到小构成
    V = Dsqrtinv*E';				 % 特征值矩阵乘以特征向量矩阵
end

在进行 ICA或者ISA处理之前，首先要对训练样本图像集进行预处理，包括去均值和白化，分为两步进行：首先要对观测数据阵X进行主成分分析，其次对X进行白化。

1. 对X进行主成分分析

   求X的协方差矩阵，如下式：
   $$C=X\ast X^{{}^{\prime }}$$
   对C进行奇异值分解，如下式：
   $$C=U\ast D\ast U{}^{\prime }$$
   其中，D为特征值的对角阵，U为C的奇异值分解中的左奇异阵，U的各个分量$[u_1,u_2,...,u_m]$为C的特征向量。

2. 对X进行白化

   观测数据阵X的白化过程如下：
   $$Z=M\ast X$$
   其中，M为白化矩阵，$M=\text{inv}(\sqrt{D})\ast U^{\prime }$，其中D为:X的协方差矩阵C的奇异值分解中的特征根对角阵。Z为白化后的数据阵。