主分量分析(PCA).doc

实验1:线性PCA算法
一． 算法原理：
寻找一组正交基组成的矩阵，有，使得是对角阵。则P的行向量（也就是一组正交基），就是数据的主元向量。
对进行推导：

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

则是一个对角阵而则是对称阵的特征向量排成的矩阵。
这里要提出的一点是，是一个的矩阵，而它将有个特征向量。其中是矩阵的秩。如果，则即为退化阵。此时分解出的特征向量不能覆盖整个空间。此时只需要在保证基的正交性的前提下，在剩余的空间中任意取得维正交向量填充的空格即可。它们将不对结果造成影响。因为此时对应于这些特征向量的特征值，也就是方差值为零。
求出特征向量矩阵后我们取，则，由线形代数可知矩阵有性质，从而进行如下计算：

可知此时的就是我们需要求得变换基。至此我们可以得到PCA的结果：
l 的主元即是的特征向量，也就是矩阵的行向量。
l 矩阵对角线上第i个元素是数据在方向的方差。

计算PCA求解的一般步骤：
（1）采集数据形成的矩阵。为观测变量个数，为采样点个数。
（2）在每个观测变量（矩阵行向量）上减去该观测变量的平均值得到矩阵。
（3）对的协方差阵进行特征分解：

式中：是对角阵，，各为的特征根，为特征矩阵，它的各列为特征矢量。求出特征根后和特征矩阵后，对特征根进行重新排列，使得，特征矢量进行相应的交换。
（4）把前乘到数据阵上，得：
的各行即为的主分量，他们在中是依能量大小排列的。

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

其中是一个的矩阵，是一个的矩阵，而是的对角阵。形式如下：

其中，是原矩阵的奇异值。由简单推导可知，如果对奇异值分解加以约束：的向量必须正交，则矩阵即为PCA的特征值分解中的，则说明PCA并不一定需要求取，也可以直接对原数据矩阵进行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中的观察信号相比，有了很大的改进，能够看出信号中含有方波，锯齿波和随机噪声。