PCA算法的主要步驟是：
讓客戶滿意是我們工作的目標(biāo)，不斷超越客戶的期望值來自于我們對這個行業(yè)的熱愛。我們立志把好的技術(shù)通過有效、簡單的方式提供給客戶，將通過不懈努力成為客戶在信息化領(lǐng)域值得信任、有價值的長期合作伙伴，公司提供的服務(wù)項(xiàng)目有：域名申請、網(wǎng)頁空間、營銷軟件、網(wǎng)站建設(shè)、忻城網(wǎng)站維護(hù)、網(wǎng)站推廣。
（1） 對向量X進(jìn)行去中心化
（2） 計算向量X的協(xié)方差矩陣，自由度可以選擇0或1
（3）計算協(xié)方差矩陣的特征值和特征向量
（4）選取最大的k個特征值及其特征向量
（5）用X與特征向量相乘
python實(shí)現(xiàn)：
from sklearn.datasets import load_iris
import numpy as np
def pca(X, k):
X = X - X.mean(axis=0)
X_cov = np.cov(X.T, ddof = 0)
eigenvalues, eigenvectors = eig(X_cov)
klarge_index = eigenvalues.argsort()[-k:][::-1]
k_eigenvectors = eigenvectors[klarge_index]
return np.dor(X, k_eigenvectors.T)
iris = load_iris()
X = iris.data
k = 2
X_pca = pca(X, k)
程序說明：
y = pca(mixedsig)，程序中mixedsig為 n*T 階混合數(shù)據(jù)矩陣，n為信號個數(shù)，T為采樣點(diǎn)數(shù), y為 m*T 階主分量矩陣。
程序設(shè)計步驟：
1、去均值
2、計算協(xié)方差矩陣及其特征值和特征向量
3、計算協(xié)方差矩陣的特征值大于閾值的個數(shù)
4、降序排列特征值
5、去掉較小的特征值
6、去掉較大的特征值(一般沒有這一步）
7、合并選擇的特征值
8、選擇相應(yīng)的特征值和特征向量
9、計算白化矩陣
10、提取主分量
程序代碼
%程序說明：y = pca(mixedsig)，程序中mixedsig為 n*T 階混合數(shù)據(jù)矩陣，n為信號個數(shù)，T為采樣點(diǎn)數(shù)
% y為 m*T 階主分量矩陣。
function y = pca(mixedsig)
if nargin == 0
error('You must supply the mixed data as input argument.');
end
if length(size(mixedsig))2
error('Input data can not have more than two dimensions. ');
end
if any(any(isnan(mixedsig)))
error('Input data contains NaN''s.');
end
%——————————————去均值————————————
meanValue = mean(mixedsig')';
mixedsig = mixedsig - meanValue * ones(1,size(meanValue,2));
[Dim,NumofSampl] = size(mixedsig);
oldDimension = Dim;
fprintf('Number of signals: %d\n',Dim);
fprintf('Number of samples: %d\n',NumofSampl);
fprintf('Calculate PCA...');
firstEig = 1;
lastEig = Dim;
covarianceMatrix = cov(mixedsig',1); %計算協(xié)方差矩陣
[E,D] = eig(covarianceMatrix); %計算協(xié)方差矩陣的特征值和特征向量
%———計算協(xié)方差矩陣的特征值大于閾值的個數(shù)lastEig———
rankTolerance = 1e-5;
maxLastEig = sum(diag(D)) rankTolerance;
lastEig = maxLastEig;
%——————————降序排列特征值——————————
eigenvalues = flipud(sort(diag(D)));
%—————————去掉較小的特征值——————————
if lastEig oldDimension
lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig + 1))/2;
else
lowerLimitValue = eigenvalues(oldDimension) - 1;
end
lowerColumns = diag(D) lowerLimitValue;
%—————去掉較大的特征值(一般沒有這一步)——————
if firstEig 1
higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig))/2;
else
higherLimitValue = eigenvalues(1) + 1;
end
higherColumns = diag(D) higherLimitValue;
%—————————合并選擇的特征值——————————
selectedColumns =lowerColumns higherColumns;
%—————————輸出處理的結(jié)果信息—————————
fprintf('Selected[ %d ] dimensions.\n',sum(selectedColumns));
fprintf('Smallest remaining (non-zero) eigenvalue[ %g ]\n',eigenvalues(lastEig));
fprintf('Largest remaining (non-zero) eigenvalue[ %g ]\n',eigenvalues(firstEig));
fprintf('Sum of removed eigenvalue[ %g ]\n',sum(diag(D) .* (~selectedColumns)));
%———————選擇相應(yīng)的特征值和特征向量———————
E = selcol(E,selectedColumns);
D = selcol(selcol(D,selectedColumns)',selectedColumns);
%——————————計算白化矩陣———————————
whiteningMatrix = inv(sqrt(D)) * E';
dewhiteningMatrix = E * sqrt(D);
%——————————提取主分量————————————
y = whiteningMatrix * mixedsig;
%——————————行選擇子程序———————————
function newMatrix = selcol(oldMatrix,maskVector)
if size(maskVector,1)~ = size(oldMatrix,2)
error('The mask vector and matrix are of uncompatible size.');
end
numTaken = 0;
for i = 1:size(maskVector,1)
if maskVector(i,1) == 1
takingMask(1,numTaken + 1) == i;
numTaken = numTaken + 1;
end
end
newMatrix = oldMatrix(:,takingMask);
本文來自CSDN博客，轉(zhuǎn)載請標(biāo)明出處：
/* --- 演示 STC 1T 系列單片機(jī) 用PCA功能實(shí)現(xiàn)16位定時器 --------*/
#include "reg51.h"
#include "intrins.h"
#define FOSC 18432000L
#define T100Hz (FOSC / 12 / 100)
typedef unsigned char BYTE;
typedef unsigned int WORD;
/*Declare SFR associated with the PCA */
sfr CCON = 0xD8; //PCA control register
sbit CCF0 = CCON^0; //PCA module-0 interrupt flag
sbit CCF1 = CCON^1; //PCA module-1 interrupt flag
sbit CR = CCON^6; //PCA timer run control bit
sbit CF = CCON^7; //PCA timer overflow flag
sfr CMOD = 0xD9; //PCA mode register
sfr CL = 0xE9; //PCA base timer LOW
sfr CH = 0xF9; //PCA base timer HIGH
sfr CCAPM0 = 0xDA; //PCA module-0 mode register
sfr CCAP0L = 0xEA; //PCA module-0 capture register LOW
sfr CCAP0H = 0xFA; //PCA module-0 capture register HIGH
sfr CCAPM1 = 0xDB; //PCA module-1 mode register
sfr CCAP1L = 0xEB; //PCA module-1 capture register LOW
sfr CCAP1H = 0xFB; //PCA module-1 capture register HIGH
sfr PCAPWM0 = 0xf2;
sfr PCAPWM1 = 0xf3;
sbit PCA_LED = P1^0; //PCA test LED
BYTE cnt;
WORD value;
void PCA_isr() interrupt 7 using 1
{
CCF0 = 0; //Clear interrupt flag
CCAP0L = value;
CCAP0H = value 8; //Update compare value
value += T100Hz;
if (cnt-- == 0)
{
cnt = 100; //Count 100 times
PCA_LED = !PCA_LED; //Flash once per second
}
}
void main()
{
CCON = 0; //Initial PCA control register
//PCA timer stop running
//Clear CF flag
//Clear all module interrupt flag
CL = 0; //Reset PCA base timer
CH = 0;
CMOD = 0x00; //Set PCA timer clock source as Fosc/12
//Disable PCA timer overflow interrupt
value = T100Hz;
CCAP0L = value;
CCAP0H = value 8; //Initial PCA module-0
value += T100Hz;
CCAPM0 = 0x49; //PCA module-0 work in 16-bit timer mode
//and enable PCA interrupt
CR = 1; //PCA timer start run
EA = 1;
cnt = 0;
while (1);
}
? PCA的處理步驟：
1，均值化
2，求協(xié)方差矩陣（我知道的有兩種方法，這是第一種，按部就班的求，第二種是：（A*A‘/(N-1)））
3，求協(xié)方差的特征值和特征向量
4，將特征值按照從大到小的順序排序，選擇其中最大的k個，然后將其對應(yīng)的k個特征向量分別作為列向量組成特征向量矩陣
5，將樣本點(diǎn)投影到選取的特征向量上
matlab實(shí)現(xiàn)源代碼
%PCA算法，matlab實(shí)現(xiàn)
function?F=pcad(A,n)%A是M*N
%測試實(shí)例A=[2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1;2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9]
%結(jié)果F=[0.8280，-1.7776，0.9922，0.2742，1.6758，0.9129，-0.0991，-1.1446，-0.4380，-1.2238]
%PCA第一步：均值化
X=A-repmat(mean(A,2),1,size(A,2))%去均值
%PCA第二步：求特征協(xié)方差矩陣
B=COV(X')%求協(xié)方差
%PCA第三步：求特征協(xié)方差矩陣的特征值和特征向量
[v,d]=eig(B)%求特征值和特征向量
%PCA第四步：將特征值按照從大到小的順序排序
d1=diag(d);%取出對角矩陣，也就是把特征值提出來組成一個新的M*1的d1矩陣
[d2?index]=sort(d1);?%特征值以升序排序?d2是排序后的結(jié)果?index是數(shù)排序以前的排名位置
cols=size(v,2);%?特征向量矩陣的列數(shù)
for?i=1:cols???%對特征向量做相反位置的調(diào)整?是個降序排列。這個過程把特征值和特征向量同時做相應(yīng)的降序排列
vsort(:,i)?=?v(:,index(cols-i+1)?);?%?vsort?是一個M*col(注:col一般等于M)階矩陣，保存的是按降序排列的特征向量,每一列構(gòu)成一個特征向量
%vsort保存的是協(xié)方差矩陣降序后的特征向量，為M*M階
dsort(i)?=?d1(index(cols-i+1));??%?dsort?保存的是按降序排列的特征值，是一維行向量，1*M
end??%完成降序排列
M=vsort(:,1:n)%提取主成分量
%PCA第五步：將樣本點(diǎn)投影到選取的特征向量上
F=(X'*M)'%最終的投影
package?PCA;
import?java.util.ArrayList;
import?java.util.Collections;
import?java.util.List;
import?org.jblas.ComplexDoubleMatrix;
import?org.jblas.DoubleMatrix;
import?org.jblas.Eigen;
public?class?PCA?{
?/**
??*?Reduce?matrix?dimension?????減少矩陣維度
??*?@param?source?? 源矩陣
??*?@param?dimension? 目標(biāo)維度
??*?@return?Target?matrix? 返回目標(biāo)矩陣
??*/?
public?static?void?main(String[]?args){
DoubleMatrix?d?=?new?DoubleMatrix(new?double[][]{{-1,-1,0,2,0},
{-2,0,0,1,1}});
DoubleMatrix?result?=?PCA.dimensionReduction(d,?2);
System.out.println(result);
}
public?static?DoubleMatrix?dimensionReduction(DoubleMatrix?source,?int?dimension)?{
//C=X*X^t/m 矩陣*矩陣^異或/列數(shù)
DoubleMatrix?covMatrix?=?source.mmul(source.transpose()).div(source.columns);
ComplexDoubleMatrix?eigVal?=?Eigen.eigenvalues(covMatrix);
ComplexDoubleMatrix[]?eigVectorsVal?=?Eigen.eigenvectors(covMatrix);
ComplexDoubleMatrix?eigVectors?=?eigVectorsVal[0];
//通過特征值將符號向量從大到小排序
ListPCABean?beans?=?new?ArrayListPCA.PCABean();
for?(int?i?=?0;?i??eigVectors.columns;?i++)?{
beans.add(new?PCABean(eigVal.get(i).real(),?eigVectors.getColumn(i)));
}
Collections.sort(beans);
DoubleMatrix?newVec?=?new?DoubleMatrix(dimension,?beans.get(0).vector.rows);
for?(int?i?=?0;?i??dimension;?i++)?{
ComplexDoubleMatrix?dm?=?beans.get(i).vector;
DoubleMatrix?real?=?dm.getReal();
newVec.putRow(i,?real);
}
return?newVec.mmul(source);
}
static?class?PCABean?implements?ComparablePCABean?{
double?eigenValue;
ComplexDoubleMatrix?vector;
public?PCABean(double?eigenValue,?ComplexDoubleMatrix?vector)?{
super();
this.eigenValue?=?eigenValue;
this.vector?=?vector;
}
@Override
public?int?compareTo(PCABean?o)?{
return?Double.compare(o.eigenValue,?eigenValue);
}
@Override
public?String?toString()?{
return?"PCABean?[eigenValue="?+?eigenValue?+?",?vector="?+?vector?+?"]";
}
}
}
org.jblas的jar包去百度下一個，我不知道怎么上傳文件
private void lbl_CutImage_Paint(object sender, PaintEventArgs e)
{
int imgWidth = this.lbl_CutImage.Width - 4;
int imgHeight = this.lbl_CutImage.Height - 4;
if (imgWidth 1) { imgWidth = 1; }
if (imgHeight 1) { imgHeight = 1; }
// 創(chuàng)建緩存圖像，先將要繪制的內(nèi)容全部繪制到緩存中，最后再一次性繪制到 Label 上，
// 這樣可以提高性能，并且可以防止屏幕閃爍的問題
Bitmap bmp_lbl = new Bitmap(this.lbl_CutImage.Width, this.lbl_CutImage.Height);
Graphics g = Graphics.FromImage(bmp_lbl);
// 將要截取的部分繪制到緩存
Rectangle destRect = new Rectangle(2, 2, imgWidth, imgHeight);
Point srcPoint = this.lbl_CutImage.Location;
srcPoint.Offset(2, 2);
Rectangle srcRect = new Rectangle(srcPoint, new System.Drawing.Size(imgWidth, imgHeight));
g.DrawImage(this.screenImage, destRect, srcRect, GraphicsUnit.Pixel);
SolidBrush brush = new SolidBrush(Color.FromArgb(10, 124, 202));
Pen pen = new Pen(brush, 1.0F);
重慶分公司
重慶分公司
