主分量分析(PCA)

一． 算法原理：

寻找一组正交基组成的矩阵P，有YPX，使得CY1YYT是对角阵。

n1则P的行向量（也就是一组正交基），就是数据X的主元向量。 对CY进行推导：

1CYYYTn11(PX)(PX)Tn11PXXTPTn11P(XXT)PTn11CYPAPTn1

T 定义AXX，则A是一个对称阵。对A进行对角化求取特征向量得：

AEDET

则D是一个对角阵而E则是对称阵A的特征向量排成的矩阵。

A是一个mm的矩阵， 这里要提出的一点是，而它将有r(rm)个特征向量。

其中r是矩阵A的秩。如果rm，则A即为退化阵。此时分解出的特征向量不能覆盖整个m空间。此时只需要在保证基的正交性的前提下，在剩余的空间中任意取得mr维正交向量填充E的空格即可。它们将不对结果造成影响。因为此时对应于这些特征向量的特征值，也就是方差值为零。

TT 求出特征向量矩阵后我们取PE，则APDP，由线形代数可知矩阵P有

1T性质PP，从而进行如下计算：

1CYPAPTn11P(PTDP)PTn11(PPT)D(PPT)n11(PP1)D(PP1)n11CYDn1

可知此时的P就是我们需要求得变换基。至此我们可以得到PCA的结果：

T X的主元即是XX的特征向量，也就是矩阵P的行向量。  矩阵CY对角线上第i个元素是数据X在方向pi的方差。

计算PCA求解的一般步骤：

（1）采集数据形成MN的矩阵。M为观测变量个数，N为采样点个数。 （2）在每个观测变量（矩阵行向量）上减去该观测变量的平均值得到矩阵X。 （3）对X的协方差阵CXX'进行特征分解：

CXX'UU'

MMMNMNMMMMMM式中：是对角阵，Diag[1,2,...,M]，各i为C的特征根，U为特征矩阵，它的各列u1,u2,....,uM为特征矢量。求出特征根后和特征矩阵后，对特征根进行重新排列，使得12...M，特征矢量u1,u2,....,uM进行相应的交换。

（4）把U'前乘到数据阵X上，得：PU'X

MNMMMNP的各行即为X的主分量，他们在P中是依能量大小排列的。

PCA方法和线形代数中的奇异值分解(SVD)方法有内在的联系，一定意义上来说，PCA的解法是SVD的一种变形和弱化。对于mn的矩阵X，通过奇异值分解可以直接得到如下形式：

XUVT

其中U是一个mm的矩阵，V是一个nn的矩阵，而是mn的对角阵。形式如下：

1r00 r，是原矩阵的奇异值。由简单推导可知，如果对奇

00其中12异值分解加以约束：U的向量必须正交，则矩阵U即为PCA的特征值分解

T中的E，则说明PCA并不一定需要求取XX，也可以直接对原数据矩阵X进行SVD奇异值分解即可得到特征向量矩阵，也就是主元向量。

二． MATLAB程序：

%-------PCA程序------

% y = PCA(x)，x为 M*T 阶混合数据矩阵，M为信号个数，T为采样点数 % y为 M*T 阶主分量矩阵 function y = PCA(x)

[M,T] = size(x); %获取输入矩阵的行数和列数（信号个数，采样点数） average= mean(x')'; %均值 for i=1:M

x(i,:)=x(i,:)-average(i)*ones(1,T); end

%------------特征根和特征向量分解-------- % CovMatrix = cov(x',1); %计算协方差矩阵

% [eigvector,eigvalue] = eig(CovMatrix); %计算协方差矩阵的特征值和特征向量

% %降序排列特征值,sortedvalue为特征根，order为特征根对应的次序 % [sortedvalue,order] = sort(diag(eigvalue),'descend'); % srotedeigvector=zeros(M,M); %计算重新排列的特征向量 % for i=1:M

% srotedeigvector(:,i)=eigvector(:,order(i)); % end

%------------奇异值分解（svd）

[U S V]=svd(x); y=U'*x;

%------------主程序------ clear all;clc;

N=200;n=1:N;%N为采样点数

s1=2*sin(0.02*pi*n);%正弦信号

t=1:N;s2=2*square(100*t,50);%方波信号

a=linspace(-1,1,25);s3=2*[a,a,a,a,a,a,a,a];%¾锯齿信号 s4=rand(1,N);%噪声

s=[s1;s2;s3;s4];%信号组成分量 4*N A=rand(4,4);

x=A*s;%观察信号 %源信号波形

figure(1);subplot(4,1,1);plot(s1);axis([0 N -5,5]);title('源信号'); subplot(4,1,2);plot(s2);axis([0 N -5,5]); subplot(4,1,3);plot(s3);axis([0 N -5,5]); subplot(4,1,4);plot(s4);xlabel('Time/ms');

%混合信号波形

figure(2);subplot(4,1,1);plot(x(1,:));title('观察信号'); subplot(4,1,2);plot(x(2,:));

subplot(4,1,3);plot(x(3,:));subplot(4,1,4);plot(x(4,:));

y = PCA(x);

%主分量分析后的信号成分

figure(3);subplot(4,1,1);plot(y(1,:));title('主分量分析后的信号'); subplot(4,1,2);plot(y(2,:)); subplot(4,1,3);plot(y(3,:));

subplot(4,1,4);plot(y(4,:));xlabel('Time/ms');

[U1,U2]=princomp(x'); w=U2';

figure(4);subplot(4,1,1);plot(w(1,:));

title('主分量分析后的信号（调用库函数）'); %用于比较 subplot(4,1,2);plot(w(2,:)); subplot(4,1,3);plot(w(3,:));

subplot(4,1,4);plot(w(4,:));xlabel('Time/ms');

三． 实验结果：

图1-源信号波形图

图2-观察信号波形图

图3-PCA分析后的信号波形图

图4-调用Matlab库函数得到的主分量分析结果

从图3和图4可以看出，使用本人所编的程序与使用Matlab中的库函数得到的结果基本相同，区别仅在于分离出的信号的符号可能相反，从而说明了程序的正确性。图3中PCA分解后的得到的信号与图1中源信号相比仍有较大的区别，但与图2中的观察信号相比，有了很大的改进，能够看出信号中含有方波，锯齿波和随机噪声。