分块压缩感知图像处理仿真BCS.rar

分块压缩感知图像处理仿真BCS.rar

clc;
clear;
tic;
im=phantom（‘Modified Shepp-Logan‘256）;
G=abs（fftshift（fft2（im）））;
figure;
imshow（G[]）;
colorbar
title（‘二维Fourier变换频谱图‘）;
im=im2uint8（mat2gray（im））;
figure;
imshow（im[]）;
title（‘256*256原图像‘）;

[m n]=size（im）;
im2=zeros（mn）;
bm=12;
bn=12;	%分块大小
p=0.7;	%采样率
d=bm*bn;	%信号长度
N=round（d*p）;	%测量值个数
bx=ceil（m/bm）;
by=ceil（n/bn）;	%分块数
imc=zeros（bx*bmby*bn）;
for i=1:bx*bm	%分块没有覆盖到的地方补0
	for j=1:by*bn
		if i<=m&&j<=n
			imc（ij）=im（ij）;
		else imc（ij）=0;
		end
	end
end
for x=1:bx	%对每个小块处理
	for y=1:by
		b=imc（（1+（x-1）*bm）:（x*bm）（1+（y-1）*bn）:（y*bn））;
		xin=reshape（bd1）;
		xin=double（xin）;	%输入信号
		A=dctmtx（d）;		%DCT矩阵
		Phi=randn（Nd）;	%高斯随机矩阵
		T=Phi*inv（A）;		%传感矩阵
		s=Phi*xin;			%测量值
		hat_y=zeros（1d）;	%OMP算法
		Aug_t=[];
		aug_y=[];
		pos_array=[];
		r_n=s;
		times=1;
		while norm（r_n）>0.01
			for col=1:d;
				product（col）=abs（T（:col）‘*r_n）;
			end
			[valpos]=max（product）;
			Aug_t=[Aug_tT（:pos）];
			T（:pos）=zeros（N1）;
			aug_y=pinv（Aug_t‘*Aug_t）*Aug_t‘*s;
			r_n=s-Aug_t*aug_y;
			pos_array（times）=pos;
			times=times+1;
		end
		hat_y（pos_array）=aug_y;
		hat_x=hat_y;
		hat_x2=inv（A）*hat_x‘;	%复原信号
		b2=reshape（hat_x2bmbn）;
		%将计算好的分块组合
		imc2（（1+（x-1）*bm）:（x*bm）（1+（y-1）*bn）:（y*bn））=b2;
	end
end
imc2=uint8（imc2）;
im2=imc2（1:m1:n）;
mse=sum（sum（abs（im-im2）.^2））/（m*n）;	%计算均方误差
toc;		%计算时间
figure;	%显示图像
imshow（im2[0 255]）;
title（strcat（‘OMP重构结果:‘‘采样率=‘num2str（p）‘分块=‘...
num2str（bm）‘בnum2str（bn）‘计算时间=‘...
num2str（round（toc））‘s MSE=‘ ...
num2str（round（mse*10）/10）））;
colormap gray;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       1730  2014-12-11 19:12  BCS.m

----------- ---------  ---------- -----  ----

                 1730                    1