% function Wavelet_OMP

clc;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)

提取出来的都不对

这是啥语言