%4阶统计量方法
filename1=['E:\课堂作业\神经网络\NN\风景图\3.bmp'];
imgdata=double(imread(filename1,'bmp')); %%%读取图像3
[h,w]=size(imgdata); %%%图像的大小为[h,w]
len=10^3; %%%%选取10万个样本
for i=1:len %%%%
x=floor(rand()*(w-16))+1;
y=floor(rand()*(h-16))+1;
TI=imgdata(y:(y+15),x:(x+15)); %%%从图中随机选取，小图的大小 16*16
X(:,i)=reshape(TI,256,1); %%将每个小图整形成[256,1]的数据；
end
X=X-mean(X,2)*ones(1,len);
Cx=X*X'/len; %%样本去均值
[U V]=svd(Cx); %%%
epsilon=1e-2;
n=length(find(diag(V)>epsilon)) %%%将大于>epsilon的奇异值取出
U=U(:,1:n); %%%[256,n]
V=V(1:n,1:n);
V=diag((diag(V).^(-0.5))); %%%[n,n]
M=V*U'; %%%[n,n]*[n,256]->[n,256]
Z=M*X; %%[n,n]*[256,len]=[n,len]
sum=0;
for i=1:len
sum=sum+(norm(Z(:,i)))^2*Z(:,i)*Z(:,i)';
end %%%[n,n]
Rz4=sum/len; %%% 四阶矩大小为[n,n]
[B S]=svd(Rz4); %%B:[n,n]
A=pinv(M)*B; %%[256,n],A的每一列极为一个基；
for i=1:n
A(:,i)=255*(A(:,i)-min(A(:,i)))/(max(A(:,i))-min(A(:,i)));
end
%imgdata=uint8(imgdata);
A=uint8(A);
figure;
for i=1:n
base=reshape(A(:,i),16,16); %%%%将基转换为[16,16]的图
subplot(ceil(n/floor(sqrt(n))),fix(sqrt(n)),i),imshow(base);
end
