噪声估计程序

这一段时间，想着怎么修改噪声估计程序，以使得获得更快的速度，而且多编写了一些，使得对3维和1维的数据也能使用。

不多说以下是程序

%% 该函数估计图像(信号)的噪声水平
% 这里假设噪声为高斯白噪声
% 由于各种问题，函数只估计1、2、3维数据的噪声，一般而言，这是足够的.....
% An Efficient Statistical Method for Image Noise Level Estimation
% Guangyong Chen1, Fengyuan Zhu1, and Pheng Ann Heng1,2
function [delta]=Noisele(M,d)
M=squeeze(M);
SizeM=size(M);
%% 一维情形
if SizeM(2)==1
if nargin<2
d=65;
end
X=zeros(SizeM(1)-d,d);
for ii=1:SizeM(1)-d
F=M(ii:SizeM(1)-d+ii-1);
X(:,ii)=F(:);
end
else
if nargin<2
d=9;
end
%% 二维情形
if length(SizeM)==2
%% 将每一小斑块列化
X=zeros((SizeM(1)-d)*(SizeM(2)-d),d^2);
for ii=1:d
for jj=1:d
F=M(ii:SizeM(1)-d+ii-1,jj:SizeM(2)-d+jj-1);
X(:,jj+(ii-1)*d)=F(:);
end
end
else
%% 三维情形
X=zeros((SizeM(1)-d)*(SizeM(2)-d)*(SizeM(3)-d),d^3);
index=1;
for ii=1:d
for jj=1:d
for kk=1:d
F=M(ii:SizeM(1)-d+ii-1,jj:SizeM(2)-d+jj-1,kk:SizeM(3)-d+kk-1);
X(:,index)=F(:);
index=index+1;
end
end
end
end
end
[delta]=NoiseLevel(X);
end

function [delta]=NoiseLevel(X)
[~,mm]=size(X);
%% 求协方差矩阵
F=cov(X);
%% 求协方差矩阵的特征值并降序排列
[~,D]=eig(F);
D=real(diag(D));
D=sort(D,'descend');
%% 估计噪声大小
for ii=1:mm
t=sum(D(ii:mm))/(mm+1-ii);
F=floor((mm+ii)/2);
F1=F-1;
F2=min(F+1,mm);
if (t<=D(F1))&&(t>=D(F2))
delta=sqrt(t);
break;
end
end
end

程序的改进在于两个方面，一是使用了协方差函数cov

二是在数据的操作上，结合MATLAB的特性，减小了循环次数。

 

以下是一个多维通用版本，但速度相比之前的要慢一半

可以处理高维数据，但是噪声估计算法估计高维数据的速度会很慢

function [dt]=Noisle(M)
SizeM=size(M);
d=9;
if SizeM(2)==1
SizeM=SizeM(1);d=65;end
Length_sizem=length(SizeM);
Lengthm=length(M(:));

Numberplus=ones(Length_sizem,1);
Numberdim=(SizeM(1)-d)*ones(Length_sizem,1);
index=(1:Lengthm)';
Label=index;
for ii=1:Length_sizem
if ii>1
Numberdim(ii)=(SizeM(ii)-d)*Numberdim(ii-1);
Numberplus(ii)=(SizeM(ii-1))*Numberplus(ii-1);
end
indexi=mod(ceil(index/Numberplus(ii)),SizeM(ii));
indexi(indexi>SizeM(ii)-d)=0;
Label(indexi==0)=0;
end

Label(Label==0)=[];
dlength=d^Length_sizem;
X=zeros(Numberdim(Length_sizem),dlength);
for ii=1:dlength
Labelplus=0;fj=ii;
for jj=1:Length_sizem
mj=mod(fj,d);
if mj==0
mj=d;
end
Labelplus=Labelplus+mj*Numberplus(jj);
fj=floor(fj/d);
end
F=M(Label+Labelplus-1);
X(:,ii)=F(:);
end
W=cov(X);
[~,D]=eig(W);
D=real(diag(D));
D=sort(D,'descend');
for ii=1:dlength
t=sum(D(ii:dlength))/(dlength+1-ii);
F=floor((dlength+ii)/2);
F1=F-1;
F2=min(F+1,dlength);
if (t<=D(F1))&&(t>=D(F2))
dt=sqrt(t);
break;
end
end

@版权所有，转载请说明出处