请教：这个压缩感知的水印程序提取不出来水印，有没有大神可以帮我看看怎么修改啊。我盯着两天了，不会改额。。。
% function Wavelet_OMPclc;clear
% 读文件
X=imread('lena256.bmp');
if ndims(X)>2
X=rgb2gray(X);
end
X=imresize(X,[256 256]);
X=double(X);
[a,b]=size(X);
% 小波变换矩阵生成
ww=DWT(a);
% 小波变换让图像稀疏化（注意该步骤会耗费时间，但是会增大稀疏度）
X1=ww*sparse(X)*ww';
X1=full(X1);
% 随机矩阵生成
M=256;
R=randn(M,a);
% 得到测量矩阵，观测域结果
Y=R*X1;
cda0=Y;
cda1=cda0; % cda1 = 256_256
%嵌入水印
I=imread('watermark.bmp');
I=rgb2gray(I); %转换为灰度图像
I=double(I)/255;
I=ceil(I); %返回大于或等于I的最小整数
figure(1);
subplot(1,2,1);imshow(I);title('水印图像')
dimI=size(I);
rm=dimI(1);
cm=dimI(2);
%生成水印信息
mark=I;
alpha=100;
k1=randn(1,8);
k2=randn(1,8);
for i=1:rm % i=1:32
for j=1:cm % j=1:32
x=(i-1)*8;y=(j-1)*8;
if mark(i,j)==1
k=k1;
else
k=k2;
end
cda1(x+1,y+8)=cda0(x+1,y+8)+alpha*k(1);
cda1(x+2,y+7)=cda0(x+2,y+7)+alpha*k(2);
cda1(x+3,y+6)=cda0(x+3,y+6)+alpha*k(3);
cda1(x+4,y+5)=cda0(x+4,y+5)+alpha*k(4);
cda1(x+5,y+4)=cda0(x+5,y+4)+alpha*k(5);
cda1(x+6,y+3)=cda0(x+6,y+3)+alpha*k(6);
cda1(x+7,y+2)=cda0(x+7,y+2)+alpha*k(7);
cda1(x+8,y+1)=cda0(x+8,y+1)+alpha*k(8);
end
end
% OMP算法重构
X2=zeros(a,b); % 恢复矩阵
for i=1:b % 列循环
rec=omp(cda1(:,i),R,a);
X2(:,i)=rec;
end
% 原始图像
subplot(1,2,2);imshow(uint8(X));title('原始图像');
% 压缩传感恢复的图像
figure(2);
X3=ww'*sparse(X2)*ww; % 小波反变换得嵌入后的图像
X3=full(X3);
% 误差(PSNR)
errorx=sum(sum(abs(X3-X).^2)) % MSE误差
psnr=10*log10(255*255/(errorx/a/b)) % PSNR
a1=X3;
a_1=uint8(X3); %读入图像数据格式
imwrite(a_1,'withmark.bmp','bmp');
subplot(2,2,1),imshow(a_1,[]),title('嵌入水印后的图像');
%disp('嵌入水印处理时间');
%embed_time=cputime-start_time,
%攻击实验，测试鲁棒性
disp('对嵌入水印的图像的攻击实验，请输入选择项：');
disp('1—添加白噪声');
disp('2—高斯低通滤波');
disp('3—图像剪切');
disp('4—图像平移');
disp('5—JPEG压缩');
disp('6—直接检测水印');
disp('其他—不攻击');
d=input('请输入选择（1-6）：');
start_time=cputime;
switch d
case 1
Aimage1=a1;
noise0=20*randn(size(Aimage1));
Aimage1=Aimage1+noise0;
subplot(2,2,2);imshow(uint8(Aimage1),[]);title('加入白噪声后的图像');
M1=uint8(Aimage1);
%M_1=uint8(M1);
%imwrite(M_1,'whitenoise.bmp','bmp');
case 2
Aimage2=a1;
H=fspecial('gaussian',[4,4],0.6);
Aimage2=imfilter(Aimage2,H);
subplot(2,2,2);imshow(Aimage2,[]);title('高斯低通滤波后的图像');
M1=Aimage2;
%M_1=uint8(M1);
%imwrite(M_1,'gaussian.bmp','bmp');
case 3
Aimage3=a1;
Aimage3(1:128,1:512)=512;
subplot(2,2,2);imshow(uint8(Aimage3));title('部分剪切后的图像');
M1=Aimage3;
%M_1=uint8(M1);
%imwrite(M_1,'cutpart.bmp','bmp');
case 4
Aimage4=a1;
a=translate(strel(1),[64,64]);
Aimage_4=imdilate(Aimage4,a);
subplot(2,2,2);imshow(Aimage_4,[]);title('平移后的图像');
M1=Aimage_4;
case 5
Aimage5=a1;
Aimage5=im2double(Aimage5);
cnum=10;
dctm=dctmtx(8);
P1=dctm;
P2=dctm.';
imageDCT=blkproc(Aimage5,[8,8],'P1*x*P2',dctm,dctm.');
DCTvar=im2col(imageDCT,[8,8],'distinct').';
n=size(DCTvar,1);
DCTvar=(sum(DCTvar.*DCTvar)-(sum(DCTvar)/n).^2)/n;
[dum,order]=sort(DCTvar);
cnum=64-cnum;
mask=ones(8,8);
mask(order(1:cnum))=zeros(1,cnum);
im88=zeros(9,9);
im88(1:8,1:8)=mask;
im128128=kron(im88(1:8,1:8),ones(16));
dctm=dctmtx(8);
P1=dctm.';
P2=mask(1:8,1:8);
P3=dctm;
Aimage_5=blkproc(imageDCT,[8,8],'P1*(x.*P2)*P3',dctm.',mask(1:8,1:8),dctm);
subplot(2,2,2);imshow(Aimage_5,[]);title('经JPEG压缩后的图像');
M1=Aimage_5;
case 6
subplot(2,2,2);imshow(a_1,[]);title('未受攻击的含水印图像');
M1=a_1;
otherwise
disp('您输入的是无效数字，图像未受攻击，将直接检测水印');
subplot(2,2,2);imshow(a_1,[]);title('未受攻击的含水印图像');
M1=a_1;
end
%提取水印
%subplot(2,2,3);imshow(mark),title('原始水印图像');
%psnr_watermarked=M1;
dca1=M1;
p=zeros(1,8);
for i=1:dimI(1)
for j=1:dimI(2) %j=1:32
x=(i-1)*8;y=(j-1)*8;
p(1)=dca1(x+1,y+8);
p(2)=dca1(x+2,y+7);
p(3)=dca1(x+3,y+6);
p(4)=dca1(x+4,y+5);
p(5)=dca1(x+5,y+4);
p(6)=dca1(x+6,y+3);
p(7)=dca1(x+7,y+2);
p(8)=dca1(x+8,y+1);
%sd1=sum(sum(p.*k1))/sqrt(sum(sum(p.^2)));
%sd2=sum(sum(p.*k2))/sqrt(sum(sum(p.^2)));
%if sd1>sd2
if corr2(p,k1)>corr2(p,k2),warning off MATLAB:divideByZero;
mark1(i,j)=1;
else
mark1(i,j)=0;
end
end
end
subplot(2,2,4);imshow(mark1,[]),title('从受攻击图像中提取的水印图像');
%imwrite(mark1,'getmark.bmp','bmp');
%最后计算
%disp('攻击与提取处理时间')
%attack_recover_time=cputime-start_time
% disp('载体图像与受攻击图像峰值信噪比')
% PSNR=psnr(psnr_cover,psnr_watermarked,c,r)
disp('原水印图像与提取水印图像互相关系数')
NC=nc(mark1,mark)
