【图像隐藏】基于小波变换+SURF、RANSAC、LT码、CRC码多种算法实现图像隐藏（抗多种攻击）matlab源码

1 内容介绍

网络技术和多媒体技术迅猛发展,为了更好地保障图像信息传输的安全性和可靠性,解决数字图像的版权保护问题,图像信息隐藏技术已成为图像处理领域的研究热点之一.采用MATLAB开发环境,提出了一种图像信息隐藏算法,实现了图像信息隐藏的两大基本功能,即信息的嵌入和信息的提取.同时,模拟常见的攻击方法,验证了添加隐藏信息后的图像的抗攻击能力.实验表明,该算法具有较强的鲁棒性,添加隐藏信息后的图像能有效地应对多种常见攻击,且隐藏的信息具有良好的不可见性.

2 仿真代码

I2_2=attack1(d,I0_wm,fig);

%%%去除权威0的行和列，使本文算法可以抵抗RT、RST的组合攻击

% [ro2,co2]=size(I2_2);

%     m=0;I2_21=[];    

%     for i=1:ro2

%         if sum(uint8(I2_2(i,:)))>512

%             m=m+1;

%             I2_21(m,:)=uint8(I2_2(i,:));

%         end

%     end

%     n=0;

%     for i=1:co2

%         if sum(I2_21(:,i))>512

%             n=n+1;

%             I2_22(:,n)=I2_21(:,i);

%         end

%     end

%     I2_2=I2_22;

%     figure,imshow(uint8(I2_2));title('去除行列为0后的图像')

% I1=im2double(imread('TestImages/lena.bmp'));

% rect=[20 20 460 460];

% I2_2=imcrop(uint8(I2_0),rect);

% 先scale后rotation

% Scl=614;

% I2_1=imresize(I2_0,[Scl Scl],'bilinear');    %%无论是imrotate还是imresize采用'bicubic'是最好的

% I2_2=imrotate(I2_1,10,'bilinear');

I2=im2double(uint8(I2_2));                     %先必须化为uint8型，然后用im2double归一化

% Get the Key Points

%% 受攻击后的图像与原图像（备注：此处有误，其实应该为嵌入水印的载体图像）进行第一次特征点的匹配

Options.upright=true;

Options.tresh=0.0002;

Ipts1=OpenSurf(I1,Options);

Ipts2=OpenSurf(I2,Options);

% 挑选Scale在一定范围内的点，因尺度太大或太小，受攻击后易丢失，但受攻击后的图像不需要选择

% m=0;

% for i=1:length(Ipts1)

%     if Ipts1(i).scale>=2.0&&Ipts1(i).scale<=6.0

%         m=m+1;

%         Ipts1_impv(m)=Ipts1(i);

%     end

% end

% save('fp_mult.mat','Ipts1_impv');

% fp_tst=load('fp_mult.mat');

% Ipts1_impv=fp_tst.Ipts1_impv;

point_cnt=0;ip_tmp1=[];

for i=1:length(Ipts1)

   ip=Ipts1(i);

   

   S = 2 * fix(2.5 * ip.scale);

%    R = fix(S / 2);

   

   if S>=10.0&&S<=30

       ip_tmp1=[ip_tmp1;Ipts1(i)];

       point_cnt=point_cnt+1;

%        pt =  [(ip.x), (ip.y)];

%        ptR = [(R * cos(ip.orientation)), (R * sin(ip.orientation))];

%        

%        if(ip.laplacian >0), myPen =[0 0 1]; else myPen =[1 0 0]; end

%        

%        rectangle('Curvature', [1 1],'Position', [pt(1)-R, pt(2)-R, S, S],'EdgeColor',myPen);

%        

%        plot([pt(1), pt(1)+ptR(1)]+1,[pt(2), pt(2)+ptR(2)]+1,'g');

   end

   

   

end

% orientation_max=max(ip_tmp1.orientation);

% nc_max=max(Nc_sta(:,1));

%  tmp=Nc_sta(:,1)>=0.86*nc_max;

%  nc_coord_alt=Nc_sta(tmp,:);

[a1,idx1]=sort(abs([ip_tmp1.doh]),'descend');

ip_sort=ip_tmp1(idx1);

nn=0;ip_fp=ip_sort(1);

for i=2:length(ip_sort)

       if nn>=17

           break;

       end

       S_sort = 2 * fix(2.5 * ip_sort(i).scale);

       R_sort = fix(S_sort/2);

%        mm=0;

       kk=0;

       for j=1:length(ip_fp)

            S_fp = 2 * fix(2.5 * ip_fp(j).scale);

            R_fp = fix(S_fp/2);

%             mm=mm+1;

            if sqrt((ip_fp(j).x-ip_sort(i).x)^2+(ip_fp(j).y-ip_sort(i).y)^2)>=R_sort+R_fp

                kk=kk+1;

               

            end

            

       end

       if kk==length(ip_fp)

          nn=nn+1;

          ip_fp=[ip_fp;ip_sort(i)];

       end

end

save('fp_invT_480_lna1.mat','ip_fp');

figure,imshow(I1,[]);hold on;

for i=1:length(ip_fp)

    ip=ip_fp(i);

    

    S = 2 * fix(2.5 * ip.scale);

    R = fix(S / 2);

    

    %         if S>=12&&S<=30

    %             ip_tmp1=[ip_tmp1;ipts(i)];

    %             point_cnt=point_cnt+1;

    pt =  [(ip.x), (ip.y)];

    ptR = [(R * cos(ip.orientation)), (R * sin(ip.orientation))];

    

    if(ip.laplacian >0), myPen =[0 0 1]; else myPen =[1 0 0]; end

    

    rectangle('Curvature', [1 1],'Position', [pt(1)-R, pt(2)-R, S, S],'EdgeColor',myPen);

    

    plot([pt(1), pt(1)+ptR(1)]+1,[pt(2), pt(2)+ptR(2)]+1,'g');

    

end

fp1=load('fp_invT_480_lna1.mat');

Ipts1_impv=fp1.ip_fp;

% Ipts2_impv=[];

n=0;

for i=1:length(Ipts2)

    if Ipts2(i).scale>=2.0&&Ipts2(i).scale<=6.5

        n=n+1;

        Ipts2_impv(n)=Ipts2(i);

    end

end

% Ipts2_impv=Ipts2;

% Put the landmark descriptors in a matrix

D1 = reshape([Ipts1_impv.descriptor],64,[]);

D2 = reshape([Ipts2_impv.descriptor],64,[]);

% Find the best

% matches，最近邻次近邻比值的阈值设定法挑选，Ransac算法再进行一次挑选，是否要用到最小二乘法拟合那个矩阵？？？

err=zeros(1,length(Ipts1_impv));

% cor1=1:length(Ipts1_impv);

cor2=zeros(1,length(Ipts1_impv));

ratio=zeros(length(Ipts1_impv),1);

dis_store=zeros(length(Ipts1_impv),2);

dis_ratio=zeros(length(Ipts1_impv),5);

for i=1:length(Ipts1_impv)

    distance=sum((D2-repmat(D1(:,i),[1 length(Ipts2_impv)])).^2,1);

    dis_sort=sort(distance,'ascend');

    dis_store(i,:)=[dis_sort(1),dis_sort(2),];

    ratio(i)=dis_sort(1)/dis_sort(2);        

    [err(i),cor2(i)]=min(distance);

    dis_ratio(i,:)=[dis_store(i,:) ratio(i) i cor2(i)];

end

dis_ratio1=sortrows(dis_ratio,3);

mask1=dis_ratio1(:,3)<=0.9;

dis_ratio2=zeros(sum(mask1),5);

nn=0;

for i=1:length(dis_ratio1)

    if mask1(i)

        nn=nn+1;

        dis_ratio2(nn,:)=dis_ratio1(i,:);

    end

end

% 利用Ransac算法对匹配的点再进行一次挑选

max_itera=50;          %设置最大迭代次数

sigma=1.8;        %设置拟合矩阵[m11,m12,m13;m21,m22,m23;0,0,1]还原的偏差

pretotal=0;     %符合拟合模型的数据的个数

k=0;

% sample=round(18*rand(50,3));

% save('Index.mat','sample');

sam_tmp=load('Index.mat');

Samp=sam_tmp.sample;

while pretotal <= size(dis_ratio2,1)*0.75 &&  k

    SampIndex=Samp(k+1,:);

%     SampIndex=round(1+(size(dis_ratio2,1)-1)*rand(3,1));  %产生三个随机索引，找样本用，floor向下取整

    

    samp1_tmp=dis_ratio2(SampIndex,4);     %对原数据随机抽样两个样本

    samp2_tmp=dis_ratio2(SampIndex,5);

    

    samp1=[Ipts1_impv(samp1_tmp).y;Ipts1_impv(samp1_tmp).x];    % 此处特别注意，其实在利用SURF算法时，其中有一步

    samp2=[Ipts2_impv(samp2_tmp).y;Ipts2_impv(samp2_tmp).x];    % 通过旋转Haar小波模板求响应值，将X，Y的坐标对调了

    

    att_matr=Matrix_solv(samp1,samp2);      %对两组数据拟合出矩阵，或其他变种拟合方法

    total=0;dis_ratio3=[];

    for j=1:size(dis_ratio2,1)

        x_ini=Ipts1_impv(dis_ratio2(j,4)).y;y_ini=Ipts1_impv(dis_ratio2(j,4)).x;

        x_atk=Ipts2_impv(dis_ratio2(j,5)).y;y_atk=Ipts2_impv(dis_ratio2(j,5)).x;

        tmp=att_matr*[x_atk;y_atk;1];

        x_red=tmp(1);y_red=tmp(2);

       if sqrt((x_ini-x_red)^2+(y_ini-y_red)^2)

           total=total+1;

           dis_ratio3=[dis_ratio3;dis_ratio2(j,:)];

       end

    end

%     mask=abs(line*[data ones(size(data,1),1)]');    %求每个数据到拟合直线的距离

%     total=sum(mask

    

    if total>pretotal              %找到符合拟合矩阵数据最多的拟合矩阵

        pretotal=total;

        bestmatr=att_matr;          %找到最好的拟合矩阵

        best_disratio=dis_ratio3;   %找到最符合条件的那些坐标

        best_matr_r=att_matr;

    end  

    k=k+1;

end

best_matr_r1=inv(best_matr_r);

    

% Sort matches on vector distance

% [err, ind]=sort(err);

cor1_1=best_disratio(:,4);

cor2_1=best_disratio(:,5);

%% 法1 用左除命令（m=A\b）来寻求矩阵m的最小二乘解

% coord_ini_matr=[];coord_x_att=[];

% for i=1:length(best_disratio)

%     tmp=

%     coord_ini_matr=[coord_ini_matr;]

% tmp=ones(length(best_disratio),1);

tmp=ones(size(best_disratio,1),1);

coord_ini_matr=[[Ipts1_impv(cor1_1).y]',[Ipts1_impv(cor1_1).x]',tmp];

coord_x_att=[Ipts2_impv(cor2_1).y]';

att_matr_13=coord_ini_matr\coord_x_att;

coord_y_att=[Ipts2_impv(cor2_1).x]';

att_matr_46=coord_ini_matr\coord_y_att;

att_matr_ls=[att_matr_13';att_matr_46';0,0,1];

% T_m{num,1}=att_matr_ls;

%%% 真实的放射攻击矩阵

%%%% 法2：先用已知X坐标拟合m11,m12,m13,再用已知Y坐标拟合m21,m22,m23,并和法1的结果相比较，验证矩阵左除是否是LS解

% att_matr_13_1=lsqnonneg(coord_ini_matr,coord_x_att);

% att_matr_46_1=lsqnonneg(coord_ini_matr,coord_y_att);

% att_matr_ls1=[att_matr_13_1';att_matr_46_1';0,0,1];   %%这不是最小二乘法吧？？？和att_matr_ls差别非常大

% Make vectors with the coordinates of the best matches

Pos1=[[Ipts1_impv(cor1_1).y]',[Ipts1_impv(cor1_1).x]'];

Pos2=[[Ipts2_impv(cor2_1).y]',[Ipts2_impv(cor2_1).x]'];

% Pos1=Pos1(1:15,:);

% Pos2=Pos2(1:15,:);

% Show both images

% I = zeros([size(I1,1) size(I1,2)*2 size(I1,3)]);

% I(:,1:size(I1,2),:)=I1; I(:,size(I1,2)+1:size(I1,2)+size(I2,2),:)=I2;

% figure, imshow(I); hold on;

%%gray image

I = zeros([round(1.6*size(I1,1)) round(size(I1,2)*2.6) ]);

I(1:size(I1,1),1:size(I1,2))=I1; I(1:size(I2,1),size(I1,2)+1:size(I1,2)+size(I2,2))=I2;

figure, imshow(I,[]); hold on;

% Show the best matches

plot([Pos1(:,2) Pos2(:,2)+size(I1,2)]',[Pos1(:,1) Pos2(:,1)]','-');

plot([Pos1(:,2) Pos2(:,2)+size(I1,2)]',[Pos1(:,1) Pos2(:,1)]','o');

%% 等比例缩放与旋转、不等比例缩放与旋转、仅有不等比例缩放的判定，若为不等比例缩放与旋转的攻击，则判定先R还是先S

m11=att_matr_ls(1,1);m12=att_matr_ls(1,2);

m21=att_matr_ls(2,1);m22=att_matr_ls(2,2);

% Sca_flag_un=0;   %%不等比例缩放标志，对于先R在不等比例S，必须先还原S，然后再还原角度

% Sca_flag=0;

% RS_dist_condi1=abs(abs(m12/m11)-abs(m21/m22));  %% 到底是先不等比S还先R的区别条件

% % RS_dist_condi2=abs(abs(m11/m22)-abs(m21/m12));  %% 备注：此处在多想一下，有没有更精准的判定准则？？？

% % RS_dist_condi2 判定条件没多大意义，不等比例S、R的组合攻击都是如此

% if abs(m11-m22)>=0.05        %% 存在不等比例缩放或不等比例缩放与旋转的攻击

%     if RS_dist_condi1<=0.015 %% 判定为先R再S（不等比）Plane这幅图比较特殊，这个值需大一些

% %         [Sx,theta1]=solve('Sx*cosd(theta1)=m11','-Sx*sind(theta1)=m12','Sx','theta1');

% %         [Sy,theta2]=solve('Sy*sind(theta2)=m21','Sy*sind(theta2)=m22','Sy','theta2');

%         theta1=atand(-m12/m11); theta2=atand(m21/m22);

%         Sx_nu=sqrt(m11^2+m12^2);Sy_nu=sqrt(m21^2+m22^2);

%         Sx_nu

%         Sy_nu

%         theta=(theta1+theta2)/2;

%         theta

%         if abs(theta)<0.25

%             disp('仅存在不等比例缩放');

%             Sx_nu=m11;Sy_nu=m22;

%             Sx_nu

%             Sy_nu

%             Rot_flag=0;

%             Sca_flag=1;   %%此处还是用Sca_flag吧，为了和水印判断那段程序统一

%             [ro2,co2]=size(I2);

%             x_s_rev=round(ro2/Sx_nu);y_s_rev=round(co2/Sy_nu); 

%             I4=imresize(im2uint8(I2),[x_s_rev y_s_rev],'bicubic');  %%将不等比例缩放还原

%             Sx_rv1(num,1)=Sx_nu; Sy_rv1(num,1)=Sy_nu;

%             

%            

%         else

%             disp('存在不等比例缩放与旋转的攻击(先R再S)');

%             [ro2,co2]=size(I2);

%             x_s_rev=round(ro2/Sx_nu);y_s_rev=round(co2/Sy_nu);  %%对于先R再S的不等比例攻击必须先还原S，若先还原角度则会造成类似于Shearing的攻击

%             I2=imresize(im2uint8(I2),[x_s_rev y_s_rev],'bicubic');  %%备注：论文中如何描述、阐明？？？

%             Sca_flag=1;                               %%此处还是先用Sca_flag吧，为了和水印判断那段程序统一

%             I3_tmp=imrotate(im2uint8(I2),-theta,'bicubic','crop');

%             I3=im2double(I3_tmp);

%             figure,imshow(I3,[]),title('旋转theta角度后校正的图像');

%             Rot_flag=1;                                       %imrotate攻击标志

%             I4=im2uint8(I3);

%             Sx_rv1(num,1)=Sx_nu; Sy_rv1(num,1)=Sy_nu;

%             thet_rv1(num,1)=theta;

%         end

%     else     %% 判定为先S再R（不等比）

% %         [Sx,theta1]=solve('Sx*cosd(theta1)=m11','Sx*sind(theta1)=m21','Sx','theta1');

% %         [Sy,theta2]=solve('-Sy*sind(theta2)=m12','Sy*cosd(theta2)=m22','Sy','theta2');

%         theta1=atand(m21/m11); theta2=atand(-m12/m22);

%         Sx_nu=sqrt(m11^2+m21^2);Sy_nu=sqrt(m12^2+m22^2);

%         Sx_nu

%         Sy_nu

%         theta=(theta1+theta2)/2;

%         theta

%         if abs(theta)<0.25

%             disp('仅存在不等比例缩放');

%             Sx_nu=m11;Sy_nu=m22;

%             Sx_nu

%             Sy_nu

%             Rot_flag=0;

%             Sca_flag=1;   %%此处还是用Sca_flag吧，为了和水印判断那段程序统一

%             [ro2,co2]=size(I2);

%             x_s_rev=round(ro2/Sx_nu);y_s_rev=round(co2/Sy_nu); 

%             I4=imresize(im2uint8(I2),[x_s_rev y_s_rev],'bicubic');  %%将不等比例缩放还原 

%             Sx_rv1(num,1)=Sx_nu; Sy_rv1(num,1)=Sy_nu;

%         else

%             disp('存在不等比例缩放与旋转的攻击(先S再R)')

%             I3_tmp=imrotate(im2uint8(I2),-theta,'bicubic','crop');

%             I3=im2double(I3_tmp);

%             figure,imshow(I3,[]),title('旋转theta角度后校正的图像'); %%%先还原角度

%             Rot_flag=1;                                       %imrotate攻击标志

%             [ro3,co3]=size(I3);

%             x_s_rev=round(ro3/Sx_nu);y_s_rev=round(co3/Sy_nu); 

%             I3=imresize(im2uint8(I3),[x_s_rev y_s_rev],'bicubic');   %%%再还原缩放尺度

%             I4=I3;

%             Sca_flag=1;     %%此处还是先用Sca_flag吧，为了和水印判断那段程序统一            

%             Sx_rv1(num,1)=Sx_nu; Sy_rv1(num,1)=Sy_nu;

%             thet_rv1(num,1)=theta;

%         end

%     end

% else  %%等比例缩放或旋转或等比例缩放与旋转攻击同时有

%    

%     

%     theta1=atand(-m12/m11); theta2=atand(m21/m22);theta3=atand(m21/m11);theta4=atand(-m12/m22);

%     theta=(theta1+theta2+theta3+theta4)/4;

%     theta

% 

%     S_un1=sqrt(m11^2+m12^2);S_un2=sqrt(m21^2+m22^2);S_un3=sqrt(m11^2+m21^2);S_un4=sqrt(m12^2+m22^2);

%     S_un=(S_un1+S_un2+S_un3+S_un4)/4;

%     S_un

% 

%     if abs(theta)<=0.25

% %         figure,imshow(I2,[]);title('没受到旋转攻击的图片');

%         I3=I2;

%         Rot_flag=0;  

%     else

%         I3_tmp=imrotate(im2uint8(I2),-theta,'bicubic','crop');

%         I3=im2double(I3_tmp);

%         figure,imshow(I3,[]),title('旋转theta角度后校正的图像');

%         Rot_flag=1;                                       %imrotate攻击标志

%         thet_rv1(num,1)=theta;

%     end

%     %%判定有无等比例缩放攻击，若有则进行相应的还原

%     if abs(S_un-1.0)<=0.015

%     figure,imshow(I3);title('没受到等比例缩放攻击的图片');

%     I4=im2uint8(I3);

% %     Sca_flag=0;  

%     else    

%     [ro3,co3]=size(I3);

%     x_s_rev=round(ro3/S_un);y_s_rev=round(co3/S_un);

%     I4=imresize(im2uint8(I3),[x_s_rev y_s_rev],'bicubic');

%     figure,imshow(uint8(I4),[ ]),title('缩放或旋转缩放攻击校正后的图像');

%     Sca_flag=1;             %受到缩放攻击的标志

%     Sca_flag_un=1;

%     S_un_rv1(num,1)=S_un;

%    end

%     

% end

%% 求出仿射矩阵T，利用求逆来还原图像

att_matr_ls1=[m22,m12,0;m21,m11,0;0,0,1];

T_m{num,1}=att_matr_ls1;

T_inv=inv(att_matr_ls1);

transa1=maketform('affine',T_inv);

I4=imtransform(im2uint8(I2),transa1);%WIa为遭遇RST几何攻击后的图像

figure;imshow(uint8(I4),[]);

[ro4,co4]=size(I4);

%     I4=im2uint8(I4);

    m=0;I4_1=[];    

    for i=1:ro4

        if sum(uint8(I4(i,:)))>512*40

            m=m+1;

            I4_1(m,:)=I4(i,:);

        end

    end

    n=0;

    for i=1:co4

        if sum(I4_1(:,i))>512*40

            n=n+1;

            I4_2(:,n)=I4_1(:,i);

        end

    end

    I4=I4_2;

    

    

%% 首先在空域寻找嵌入特殊信息的位置块，以便定位到边缘被破坏的图片的边界以及再一次精准校正

    Nc_total4=[];nc_coord_tal4=[];coord_avg4=[];

    coordx=25;coordy=25;

    [Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_66_1(I4,coordx,coordy,flg_grp(1,:));

    Nc_total4=[Nc_total4;Nc_sta];nc_coord_tal4=[nc_coord_tal4;nc_coord_alt];coord_avg4=[coord_avg4;coord_x_y];

    coordx=25;coordy=463;

    [Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_lena12(I4,coordx,coordy,flg_grp(2,:));

    Nc_total4=[Nc_total4;Nc_sta];nc_coord_tal4=[nc_coord_tal4;nc_coord_alt];coord_avg4=[coord_avg4;coord_x_y];

%     coordx=247;coordy=247;

%     [Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_66(I5,coordx,coordy,flg_grp(3,:));

%     Nc_total=[Nc_total;Nc_sta];nc_coord_tal=[nc_coord_tal;nc_coord_alt];coord_avg=[coord_avg;coord_x_y];

    coordx=463;coordy=25;

    [Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_lena21(I4,coordx,coordy,flg_grp(4,:));

    Nc_total4=[Nc_total4;Nc_sta];nc_coord_tal4=[nc_coord_tal4;nc_coord_alt];coord_avg4=[coord_avg4;coord_x_y];

    coordx=463;coordy=463;

    [Nc_sta,nc_coord_alt,coord_x_y]=odd_even_ext24_66_1(I4,coordx,coordy,flg_grp(5,:));

    Nc_total4=[Nc_total4;Nc_sta];nc_coord_tal4=[nc_coord_tal4;nc_coord_alt];coord_avg4=[coord_avg4;coord_x_y];

    

    %%利用特殊位置的点相对距离求出像素的偏差，并且再一次精准校正

    x1x4_dis=coord_avg4(3,1)-coord_avg4(1,1);

    x2x5_dis=coord_avg4(4,1)-coord_avg4(2,1);

    x_devi=(438-(x1x4_dis+x2x5_dis)*0.5)/(438/512);   %若x_devi为正则需要放大，为负则需要缩小

    y1y4_dis=coord_avg4(2,2)-coord_avg4(1,2);

    y2y5_dis=coord_avg4(4,2)-coord_avg4(3,2);

    y_devi=(438-(y1y4_dis+y2y5_dis)*0.5)/(438/512);         %若y_devi为正则需要放大，为负则需要缩小

    [ro4,co4]=size(I4);

    if abs(x_devi)>=2.0&&abs(y_devi)>=2.0

        ro5=round(x_devi)+ro4;co5=round(y_devi)+co4;

        I5=imresize(I4,[ro5 co5],'bicubic');

        Scl_fc(count,1)=1;

    elseif abs(x_devi)>=2.0&&abs(y_devi)<2.0

3 运行结果

4 参考文献

[1]尤鸿霞, 郭江峰. 基于Matlab的二值图象数据隐藏算法的实现[J]. 南通纺织职业技术学院学报, 2008.

[2]庞镇明. 基于MATLAB小波变换在图像处理中算法的实现[D]. 北京理工大学.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。