精华内容
下载资源
问答
  • matlab正交匹配追踪算法
  • 基于MATLAB算法,基追踪匹配追踪算法,稀疏分解或者压缩感知算法,应用广泛,优化求解算法等,可利于初学者有效学习
  • 压缩感知的正交匹配追踪算法——matlab程序 压缩感知的正交匹配追踪算法——matlab程序
  • 匹配追踪算法OMP matlab代码

    热门讨论 2011-08-11 16:09:27
    匹配追踪算法OMP的 matlab代码,实验中常用到的子程序。很好用的哦
  • omp(正交匹配追踪算法MATLAB编写)
  • 匹配追踪Matlab算法(MP)
  • ROMP正则化的正交匹配追踪算法 ROMP算法是一种贪婪的压缩感知重建算法,在OMP的基础上加入了正则化过程。此外,相比于OMP算法,ROMP在每次迭代过程中选择与残差最相关的m个列向量。下面是Deanna Needell and Roman ...

    ROMP正则化的正交匹配追踪算法

    ROMPROMP算法是一种贪婪的压缩感知重建算法,在OMPOMP的基础上加入了正则化过程。此外,相比于OMPOMP算法,ROMPROMP在每次迭代过程中选择与残差最相关的mm个列向量。下面是DeannaNeedellDeanna NeedellRomanVershyninRoman Vershynin的文章中给出的ROMPROMP的算法。

    在正则化过程,需要寻找索引集合的所有子集,在这里可以使用MatlabMatlab中的ff2nff2n函数实现。对于正则化条件,只用满足uu中的最大元素绝对值小于最小元素绝对值的22倍即可。
    本文中用高斯随机矩阵作为测量矩阵,待观测目标信号长度NN512512,且元素值只有1,01-1,0,1,观测值个数MM120120,目标信号稀疏度KK1010MatlabMatlab实验结果如下:
    在这里插入图片描述

    clear;
    clc;
    N=512;%待观测目标信号长度
    K=10;%目标信号稀疏度
    M=120;%观测值个数
    maxnum=100;
    x=zeros(N,1);
    q=randperm(N);%随机打乱排序
    t=1:N;
    x(q(1:K))=sign(randn(K,1));%生成幅值为1-1的K个非0元素的目标信号y
    A=randn(M,N);%生成M*N的服从高斯分布的测量矩阵
    A=orth(A')';%把测量矩阵按行正交化
    y=A*x;%获得观测向量
    x2=A'*y;%基于L2范数最小化估算的初值
    S=ff2n(K);%20个元素所有可能的组合,(2^K)*K的矩阵
    I=[];
    A1=[];
    r=y;%初始化时,残差r=y
    for i=1:K
        u=A'*r;
        u1=abs(u);
        [temp,index]=sort(u1);%对向量u从小到大排序,index储存索引
        len1=length(u);
        J=index(len1-K+1:end);%储存K个的最大系数索引,此处存在问题,应该是非0
       if length(I)>=2*K
           break;
       else
           %正则化标准意思是选择各列向量与残差内积绝对值的最大值不能比最小值大两倍以上
           %(comparable coordinates)且能量最大的一组(with the maximal energy)
           [row,col]=size(S);%记录S的行数
           uu=[0;u];
           index1=S.*J';%存储所有子集在u中对应的索引
           temp2=uu(index1+1);%存储所有子集中的元素值,(2^K)个子集,一行看成一个向量
           abstemp=abs(temp2);
           max1=max(abstemp');
           abstemp1=abstemp';
           %abstemp1(find(abstemp1==0))=maxnum;
           min1=min(abstemp1(find(abstemp1~=0)));
           min1(1)=0;
           diff=abs(max1)-2*abs(min1);%最大值大于最小值两倍以上,则>0%2倍关系可以换成其他的条件吗?
           index2=find(diff>0);%寻找满足最大值小于最小值两倍以上的索引
           if isempty(find(diff>0,1))
              disp('未发现满足正则条件的集合');
              break;
           end
           temp3=temp2(index2,:);%保留满足条件的子集
           temp4=temp3.*temp3;%中间变量,存储元素平方
           temp5=sqrt(sum(temp4,2));%按行求和,再开方,即求二范数
           [maxnorm,index3]=max(temp5);%会不会存在多个最大值?
           temp6=temp3(index3(1),:);%存储第一个最大范数的坐标向量
           temp7=temp2-temp6;%若temp6是temp2中某一行,则那一行全为0
           temp8=temp7;
           temp8(1,:)=[];
           J0=find(all(temp8== 0,2));%all(temp7== 0,2)寻找全是0的行 
           J0=index1(J0+1,:);
           if isempty(intersect(I,J0))~=0
               T=intersect(I,J0);%查找重复元素,并删除
               for jj=1:length(T)
                   T1=find(J0==T(jj));
                   J0(T1)=[];
               end
           end
           I=[I J0];%更新索引集合
           A1=A(:,I);
           y2=inv(A1'*A1)*A1'*y;%最小二乘法求y2
           r=y-A1*y2;
       end
    end
    y3=zeros(N,1);
    y3(I)=y2;%储存恢复结果
    plot(t,x,'r',t,y3,'g  *');
    legend('red 代表原始信号','green 代表恢复的信号');
    length(find(y3~=0));%在能正确恢复信号时,y3中非0系数是2K
    % A1(:,1)'*A1(:,10)
    % A(:,1)'*A(:,200)
    

    Needell D , Vershynin R . Uniform Uncertainty Principle and signal recovery via Regularized Orthogonal Matching Pursuit[J]. 2007.

    展开全文
  • 有很详细的注释,关于信号的稀疏分解和重建的,不错的资料...
  • 压缩感知正交匹配追踪算法,在MATLAB中的程序。本人已经调试 完毕
  • 正交匹配追踪算法OMP(Orthogonal Matching Pursuit),压缩感知的一种典型贪婪算法,该资源给出了OMP的MATLAB的仿真代码
  • 正交匹配追踪算法OMP

    千次阅读 2019-09-01 22:34:08
    浅谈压缩感知(九):正交匹配追踪算法OMP </h1> <div class="clear"></div> <div class="postBody"> 主要内容: OMP算法介绍 OMP的MATLAB实现...

    浅谈压缩感知(九):正交匹配追踪算法OMP

            </h1>
            <div class="clear"></div>
            <div class="postBody">
    

    主要内容:

    1. OMP算法介绍
    2. OMP的MATLAB实现
    3. OMP中的数学知识

    一、OMP算法介绍

    来源:http://blog.csdn.net/scucj/article/details/7467955

    1、信号的稀疏表示(sparse representation of signals)

    给定一个过完备字典矩阵,其中它的每列表示一种原型信号的原子。给定一个信号y,它可以被表示成这些原子的稀疏线性组合。信号 y 可以被表达为 y = Dx ,或者字典矩阵中所谓过完备性,指的是原子的个数远远大于信号y的长度(其长度很显然是n),即n<<k

    2、MP算法(匹配追踪算法)

    2.1 算法描述

    作为对信号进行稀疏分解的方法之一,将信号在完备字典库上进行分解。

    假定被表示的信号为y,其长度为n。假定H表示Hilbert空间,在这个空间H里,由一组向量构成字典矩阵D,其中每个向量可以称为原子(atom),其长度与被表示信号 y 的长度n相同,而且这些向量已作为归一化处理,即|,也就是单位向量长度为1MP算法的基本思想:从字典矩阵D(也称为过完备原子库中),选择一个与信号 y 最匹配的原子(也就是某列),构建一个稀疏逼近,并求出信号残差,然后继续选择与信号残差最匹配的原子,反复迭代,信号y可以由这些原子来线性和,再加上最后的残差值来表示。很显然,如果残差值在可以忽略的范围内,则信号y就是这些原子的线性组合。如果选择与信号y最匹配的原子?如何构建稀疏逼近并求残差?如何进行迭代?我们来详细介绍使用MP进行信号分解的步骤:[1] 计算信号 y 与字典矩阵中每列(原子)的内积,选择绝对值最大的一个原子,它就是与信号 y 在本次迭代运算中最匹配的。用专业术语来描述:令信号,从字典矩阵中选择一个最为匹配的原子,满足r0 表示一个字典矩阵的列索引。这样,信号 y 就被分解为在最匹配原子的垂直投影分量和残值两部分,即:[2]对残值R1f进行步骤[1]同样的分解,那么第K步可以得到:

     其中 满足。可见,经过K步分解后,信号 y 被分解为:,其中

    2.2 继续讨论

    (1)为什么要假定在Hilbert空间中?Hilbert空间就是定义了完备的内积空。很显然,MP中的计算使用向量的内积运算,所以在在Hilbert空间中进行信号分解理所当然了。什么是完备的内积空间?篇幅有限就请自己搜索一下吧。

    (2)为什么原子要事先被归一化处理了,即上面的描述。内积常用于计算一个矢量在一个方向上的投影长度,这时方向的矢量必须是单位矢量。MP中选择最匹配的原子是,是选择内积最大的一个,也就是信号(或是残值)在原子(单位的)垂直投影长度最长的一个,比如第一次分解过程中,投影长度就是,三个向量,构成一个三角形,且正交(不能说垂直,但是可以想象二维空间这两个矢量是垂直的)。

    (3)MP算法是收敛的,因为正交,由这两个可以得出,得出每一个残值比上一次的小,故而收敛。

    2.3 MP算法的缺点

    如上所述,如果信号(残值)在已选择的原子进行垂直投影是非正交性的,这会使得每次迭代的结果并不少最优的而是次最优的,收敛需要很多次迭代。举个例子说明一下:在二维空间上,有一个信号 y D=[x1, x2]来表达,MP算法迭代会发现总是在x1x2上反复迭代,即,这个就是信号(残值)在已选择的原子进行垂直投影的非正交性导致的。再用严谨的方式描述[1]可能容易理解:Hilbert空间H中,,定义,就是它是这些向量的张成中的一个,MP构造一种表达形式:;这里的Pvf表示 fV上的一个正交投影操作,那么MP算法的第 k 次迭代的结果可以表示如下(前面描述时信号为y,这里变成f了,请注意)

    如果  是最优的k项近似值,当且仅当。由于MP仅能保证,所以一般情况下是次优的。这是什么意思呢?k个项的线性表示,这个组合的值作为近似值,只有在第k个残差和正交,才是最优的。如果第k个残值与正交,意味这个残值与fk的任意一项都线性无关,那么第k个残值在后面的分解过程中,不可能出现fk中已经出现的项,这才是最优的。而一般情况下,不能满足这个条件,MP一般只能满足第k个残差和xk正交,这也就是前面为什么提到"信号(残值)在已选择的原子进行垂直投影是非正交性的"的原因。如果第k个残差和fk不正交,那么后面的迭代还会出现fk中已经出现的项,很显然fk就不是最优的,这也就是为什么说MP收敛就需要更多次迭代的原因。不是说MP一定得到不到最优解,而且其前面描述的特性导致一般得到不到最优解而是次优解。那么,有没有办法让第k个残差与正交,方法是有的,这就是下面要谈到的OMP算法。

    3、OMP算法

    3.1 算法描述

    OMP算法的改进之处在于:在分解的每一步对所选择的全部原子进行正交化处理,这使得在精度要求相同的情况下,OMP算法的收敛速度更快。

    那么在每一步中如何对所选择的全部原子进行正交化处理呢?在正式描述OMP算法前,先看一点基础思想。

    先看一个 k  阶模型,表示信号 f 经过 k 步分解后的情况,似乎很眼熟,但要注意它与MP算法不同之处,它的残值与前面每个分量正交,这就是为什么这个算法多了一个正交的原因,MP中仅与最近选出的的那一项正交。

    (1)

     k + 1 阶模型如下:

    (2)

    应用 k + 1阶模型减去k 阶模型,得到如下:

    (3)

    我们知道,字典矩阵D的原子是非正交的,引入一个辅助模型,它是表示对前k个项的依赖,描述如下:

    (4)

    和前面描述类似,span(x1, ...xk)之一上的正交投影操作,后面的项是残值。这个关系用数学符号描述:

    请注意,这里的 a b 的上标表示第 k 步时的取值。

    (4)带入(3)中,有:

    5

    如果一下两个式子成立,(5)必然成立。

    (6)

    (7)

    ,有

    其中

    ak的值是由求法很简单,通过对(7)左右两边添加作内积消减得到:

    后边的第二项因为它们正交,所以为0,所以可以得出ak的第一部分。对于,在(4)左右两边中与作内积,可以得到ak的第二部分。

    对于(4),可以求出,求的步骤请参见参考文件的计算细节部分。为什么这里不提,因为后面会介绍更简单的方法来计算。

    3.2 收敛性证明

    通过(7),由于正交,将两个残值移到右边后求二范的平方,并将ak的值代入可以得到:

    可见每一次残差比上一次残差小,可见是收敛的。

    3.3 算法步骤

    整个OMP算法的步骤如下:

    由于有了上面的来龙去脉,这个算法就相当好理解了。

    到这里还不算完,后来OMP的迭代运算用另外一种方法可以计算得知,有位同学的论文[2]描述就非常好,我就直接引用进来:

    对比中英文描述,本质都是一样,只是有细微的差别。这里顺便贴出网一哥们写的OMP算法的代码,源出处不得而知,共享给大家。

    二、OMP的MATLAB实现

     1、一维信号重建

    代码:

    复制代码
    %  1-D信号压缩传感的实现(正交匹配追踪法Orthogonal Matching Pursuit)
    %  测量数M>=K*log(N/K),K是稀疏度,N信号长度,可以近乎完全重构
    %  编程人--香港大学电子工程系 沙威  Email: wsha@eee.hku.hk
    %  编程时间:2008年11月18日
    %  文档下载: http://www.eee.hku.hk/~wsha/Freecode/freecode.htm 
    %  参考文献:Joel A. Tropp and Anna C. Gilbert 
    %  Signal Recovery From Random Measurements Via Orthogonal Matching
    %  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
    %  DECEMBER 2007.
    

    clc;clear

    %% 1. 时域测试信号生成
    K
    =7; % 稀疏度(做FFT可以看出来)
    N
    =256; % 信号长度
    M
    =64; % 测量数(M>=Klog(N/K),至少40,但有出错的概率)
    f1
    =50; % 信号频率1
    f2
    =100; % 信号频率2
    f3
    =200; % 信号频率3
    f4
    =400; % 信号频率4
    fs
    =800; % 采样频率
    ts
    =1/fs; % 采样间隔
    Ts
    =1:N; % 采样序列
    x
    =0.3cos(2pif1Ts
    ts)+0.6cos(2pif2Tsts)+0.1cos(2pif3Tsts)+0.9cos(2pif4Ts*ts); % 完整信号,由4个信号叠加而来

    %% 2. 时域信号压缩传感
    Phi
    =randn(M,N); % 测量矩阵(高斯分布白噪声)64256的扁矩阵,Phi也就是文中说的D矩阵
    s
    =Phi
    x.; % 获得线性测量 ,s相当于文中的y矩阵

    %% 3. 正交匹配追踪法重构信号(本质上是L_1范数最优化问题)
    %匹配追踪:找到一个其标记看上去与收集到的数据相关的小波;在数据中去除这个标记的所有印迹;不断重复直到我们能用小波标记“解释”收集到的所有数据。

    m=2K; % 算法迭代次数(m>=K),设x是K-sparse的
    Psi
    =fft(eye(N,N))/sqrt(N); % 傅里叶正变换矩阵
    T
    =Phi
    Psi; % 恢复矩阵(测量矩阵*正交反变换矩阵)

    hat_y
    =zeros(1,N); % 待重构的谱域(变换域)向量
    Aug_t
    =[]; % 增量矩阵(初始值为空矩阵)
    r_n
    =s; % 残差值

    for times=1:m; % 迭代次数(有噪声的情况下,该迭代次数为K)
    for col=1:N; % 恢复矩阵的所有列向量
    product(col)
    =abs(T(:,col)*r_n); % 恢复矩阵的列向量和残差的投影系数(内积值)
    end
    [val,pos]
    =max(product); % 最大投影系数对应的位置,即找到一个其标记看上去与收集到的数据相关的小波
    Aug_t
    =[Aug_t,T(:,pos)]; % 矩阵扩充

    T(:,pos)</span>=zeros(M,<span style="color: #800080;">1</span>);                          %<span style="color: #000000;">  选中的列置零(实质上应该去掉,为了简单我把它置零),在数据中去除这个标记的所有印迹
    aug_y</span>=(Aug_t<span style="color: #800000;">'</span><span style="color: #800000;">*Aug_t)^(-1)*Aug_t</span><span style="color: #800000;">'</span>*s;           %<span style="color: #000000;">  最小二乘,使残差最小
    r_n</span>=s-Aug_t*aug_y;                            %<span style="color: #000000;">  残差
    pos_array(times)</span>=pos;                         %<span style="color: #000000;">  纪录最大投影系数的位置
    

    end
    hat_y(pos_array)=aug_y; % 重构的谱域向量
    hat_x
    =real(Psi*hat_y.); % 做逆傅里叶变换重构得到时域信号

    %% 4. 恢复信号和原始信号对比
    figure(
    1);
    hold on;
    plot(hat_x,
    k.-) % 重建信号
    plot(x,
    r) % 原始信号
    legend(
    Recovery,Original)
    norm(hat_x.
    -x)/norm(x) % 重构误差

    复制代码

    结果:

    2、二维信号重建

    代码:

    复制代码
    %  本程序实现图像LENA的压缩传感
    %  程序作者:沙威,香港大学电气电子工程学系,wsha@eee.hku.hk
    %  算法采用正交匹配法,参考文献 Joel A. Tropp and Anna C. Gilbert 
    %  Signal Recovery From Random Measurements Via Orthogonal Matching
    %  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
    %  DECEMBER 2007.
    %  该程序没有经过任何优化
    

    function Wavelet_OMP

    clc;clear

    % 读文件
    X
    =imread(lena256.bmp);
    X
    =double(X);
    [a,b]
    =size(X);

    % 小波变换矩阵生成
    ww
    =DWT(a);

    % 小波变换让图像稀疏化(注意该步骤会耗费时间,但是会增大稀疏度)
    X1
    =ww*sparse(X)*ww;
    % X1=X;
    X1
    =full(X1);

    % 随机矩阵生成
    M
    =190;
    R
    =randn(M,a);
    % R=mapminmax(R,0,255);
    % R
    =round®;

    % 测量值
    Y
    =R*X1;

    % OMP算法
    % 恢复矩阵
    X2
    =zeros(a,b);
    % 按列循环
    for i=1:b
    % 通过OMP,返回每一列信号对应的恢复值(小波域)
    rec
    =omp(Y(:,i),R,a);
    % 恢复值矩阵,用于反变换
    X2(:,i)
    =rec;
    end

    % 原始图像
    figure(
    1);
    imshow(uint8(X));
    title(
    原始图像);

    % 变换图像
    figure(
    2);
    imshow(uint8(X1));
    title(
    小波变换后的图像);

    % 压缩传感恢复的图像
    figure(
    3);
    % 小波反变换
    X3
    =ww*sparse(X2)*ww;
    % X3=X2;
    X3
    =full(X3);
    imshow(uint8(X3));
    title(
    恢复的图像);

    % 误差(PSNR)
    % MSE误差
    errorx
    =sum(sum(abs(X3-X).^2));
    % PSNR
    psnr
    =10log10(255255/(errorx/a/b))

    % OMP的函数
    % s-测量;T-观测矩阵;N-向量大小
    function hat_y
    =omp(s,T,N)
    Size
    =size(T); % 观测矩阵大小
    M
    =Size(1); % 测量
    hat_y
    =zeros(1,N); % 待重构的谱域(变换域)向量
    Aug_t
    =[]; % 增量矩阵(初始值为空矩阵)
    r_n
    =s; % 残差值

    for times=1:M; % 迭代次数(稀疏度是测量的1/4)
    for col=1:N; % 恢复矩阵的所有列向量
    product(col)
    =abs(T(:,col)*r_n); % 恢复矩阵的列向量和残差的投影系数(内积值)
    end
    [val,pos]
    =max(product); % 最大投影系数对应的位置
    Aug_t
    =[Aug_t,T(:,pos)]; % 矩阵扩充
    T(:,pos)
    =zeros(M,1); % 选中的列置零(实质上应该去掉,为了简单我把它置零)
    aug_y
    =(Aug_t*Aug_t)^(-1)*Aug_ts; % 最小二乘,使残差最小
    r_n
    =s-Aug_t
    aug_y; % 残差
    pos_array(times)
    =pos; % 纪录最大投影系数的位置

    </span><span style="color: #0000ff;">if</span> (norm(r_n)&lt;<span style="color: #800080;">0.9</span>)                            %<span style="color: #000000;">  残差足够小
        </span><span style="color: #0000ff;">break</span><span style="color: #000000;">;
    end
    

    end
    hat_y(pos_array)=aug_y; % 重构的向量

    复制代码

    结果:

    三、OMP算法中的数学知识

    算法会有一些线性代数中的概念和知识,主要是关于子空间、最小二乘、投影矩阵等,详情请参考

    最小二乘的几何意义及投影矩阵

    <div id="blog_post_info">
    
    4
    0
    <div class="clear"></div>
    <div id="post_next_prev">
    
    <a href="https://www.cnblogs.com/AndyJee/p/5045668.html" class="p_n_p_prefix">« </a> 上一篇:    <a href="https://www.cnblogs.com/AndyJee/p/5045668.html" title="发布于 2015-12-14 16:29">浅谈压缩感知(八):两篇科普文章</a>
    <br>
    <a href="https://www.cnblogs.com/AndyJee/p/5048235.html" class="p_n_p_prefix">» </a> 下一篇:    <a href="https://www.cnblogs.com/AndyJee/p/5048235.html" title="发布于 2015-12-15 14:30">浅谈压缩感知(十):范数与稀疏性</a>
    
    posted @ 2015-12-15 09:29 AndyJee 阅读(28779) 评论(1) 编辑 收藏
    </div><!--end: topics 文章、评论容器-->
    

      
    			</div>
    

    #1楼

        <span id="comment-maxId" style="display:none">3718412</span>
        <span id="comment-maxDate" style="display:none">2017/6/20 下午12:59:28</span>
    

    2017-06-20 12:59

            <a id="a_comment_author_3718412" href="https://home.cnblogs.com/u/1178954/" target="_blank">数理小先生</a>
    
    		</div>
    		<div class="feedbackCon">
    
    楼主在吗?急需要求助
    		</div>
    	</div>
    
        <div id="google_ads_iframe_/1090369/C2_0__container__" style="border: 0pt none;"><iframe id="google_ads_iframe_/1090369/C2_0" title="3rd party ad content" name="google_ads_iframe_/1090369/C2_0" width="468" height="60" scrolling="no" marginwidth="0" marginheight="0" frameborder="0" srcdoc="" style="border: 0px; vertical-align: bottom;" data-google-container-id="2" data-load-complete="true"></iframe></div></div>
    </div>
    <div id="under_post_kb">
    
    展开全文
  • 压缩感知中迂回式匹配追踪算法 2014年9期 计算机研究与发展 配套 论文提出算法和对比算法代码 (matlab
  • 针对实际中未知稀疏度信号的重建问题,提出了一种自适应的弱选择压缩采样匹配追踪算法。该算法将自适应思想、弱选择思想与CoSaMP算法相结合,在预选阶段后利用限制性弱选择策略对候选集进行二次筛选,通过双迭代阈值...
  • 打包了匹配追踪算法OMP和自适应匹配追踪SAMP算法。基于MATLAB实现,亲测可以直接运行,结果正确,适用于科研、学习、工程应用。
  • omp正交匹配追踪算法

    2013-11-19 16:54:22
    该资源详细描述了OMP代码的matlab程序和c语言程序(矩阵的求逆采用LU分解法),并且对两者结果进行了比较,恢复的信号可以精确到小数点5位,误差非常小,测量矩阵采用随机高斯矩阵,程序里面还有matlab和c语言版对...
  • matlab代码,数据源是VCChen的MIG25飞机的ISAR仿真数据,数据通过两维FFT只可以直接成像,本算法基于压缩感知,采用OMP,实现ISAR成像算法,在欠采样的条件下,也可以很好的成像
  • 浅谈压缩感知(九):正交匹配追踪算法OMP 主要内容: OMP算法介绍 OMP的MATLAB实现 OMP中的数学知识 一、OMP算法介绍 来源:http://blog.csdn.net/scucj/article/details/7467955 1、信号的稀疏表示(sparse...

    浅谈压缩感知(九):正交匹配追踪算法OMP

    主要内容:

    1. OMP算法介绍
    2. OMP的MATLAB实现
    3. OMP中的数学知识

    一、OMP算法介绍

    来源:http://blog.csdn.net/scucj/article/details/7467955

    1、信号的稀疏表示(sparse representation of signals)

    给定一个过完备字典矩阵,其中它的每列表示一种原型信号的原子。给定一个信号y,它可以被表示成这些原子的稀疏线性组合。信号 y 可以被表达为 y = Dx ,或者字典矩阵中所谓过完备性,指的是原子的个数远远大于信号y的长度(其长度很显然是n),即n<<k

    2、MP算法(匹配追踪算法)

    2.1 算法描述

    作为对信号进行稀疏分解的方法之一,将信号在完备字典库上进行分解。

    假定被表示的信号为y,其长度为n。假定H表示Hilbert空间,在这个空间H里,由一组向量构成字典矩阵D,其中每个向量可以称为原子(atom),其长度与被表示信号 y 的长度n相同,而且这些向量已作为归一化处理,即|,也就是单位向量长度为1MP算法的基本思想:从字典矩阵D(也称为过完备原子库中),选择一个与信号 y 最匹配的原子(也就是某列),构建一个稀疏逼近,并求出信号残差,然后继续选择与信号残差最匹配的原子,反复迭代,信号y可以由这些原子来线性和,再加上最后的残差值来表示。很显然,如果残差值在可以忽略的范围内,则信号y就是这些原子的线性组合。如果选择与信号y最匹配的原子?如何构建稀疏逼近并求残差?如何进行迭代?我们来详细介绍使用MP进行信号分解的步骤:[1] 计算信号 y 与字典矩阵中每列(原子)的内积,选择绝对值最大的一个原子,它就是与信号 y 在本次迭代运算中最匹配的。用专业术语来描述:令信号,从字典矩阵中选择一个最为匹配的原子,满足r0 表示一个字典矩阵的列索引。这样,信号 y 就被分解为在最匹配原子的垂直投影分量和残值两部分,即:[2]对残值R1f进行步骤[1]同样的分解,那么第K步可以得到:

     其中 满足。可见,经过K步分解后,信号 y 被分解为:,其中

    2.2 继续讨论

    (1)为什么要假定在Hilbert空间中?Hilbert空间就是定义了完备的内积空。很显然,MP中的计算使用向量的内积运算,所以在在Hilbert空间中进行信号分解理所当然了。什么是完备的内积空间?篇幅有限就请自己搜索一下吧。

    (2)为什么原子要事先被归一化处理了,即上面的描述。内积常用于计算一个矢量在一个方向上的投影长度,这时方向的矢量必须是单位矢量。MP中选择最匹配的原子是,是选择内积最大的一个,也就是信号(或是残值)在原子(单位的)垂直投影长度最长的一个,比如第一次分解过程中,投影长度就是,三个向量,构成一个三角形,且正交(不能说垂直,但是可以想象二维空间这两个矢量是垂直的)。

    (3)MP算法是收敛的,因为正交,由这两个可以得出,得出每一个残值比上一次的小,故而收敛。

    2.3 MP算法的缺点

    如上所述,如果信号(残值)在已选择的原子进行垂直投影是非正交性的,这会使得每次迭代的结果并不少最优的而是次最优的,收敛需要很多次迭代。举个例子说明一下:在二维空间上,有一个信号 y  D=[x1, x2]来表达,MP算法迭代会发现总是在x1x2上反复迭代,即,这个就是信号(残值)在已选择的原子进行垂直投影的非正交性导致的。再用严谨的方式描述[1]可能容易理解:Hilbert空间H中,,定义,就是它是这些向量的张成中的一个,MP构造一种表达形式:;这里的Pvf表示 fV上的一个正交投影操作,那么MP算法的第 k 次迭代的结果可以表示如下(前面描述时信号为y,这里变成f了,请注意)

    如果  是最优的k项近似值,当且仅当。由于MP仅能保证,所以一般情况下是次优的。这是什么意思呢?k个项的线性表示,这个组合的值作为近似值,只有在第k个残差和正交,才是最优的。如果第k个残值与正交,意味这个残值与fk的任意一项都线性无关,那么第k个残值在后面的分解过程中,不可能出现fk中已经出现的项,这才是最优的。而一般情况下,不能满足这个条件,MP一般只能满足第k个残差和xk正交,这也就是前面为什么提到"信号(残值)在已选择的原子进行垂直投影是非正交性的"的原因。如果第k个残差和fk不正交,那么后面的迭代还会出现fk中已经出现的项,很显然fk就不是最优的,这也就是为什么说MP收敛就需要更多次迭代的原因。不是说MP一定得到不到最优解,而且其前面描述的特性导致一般得到不到最优解而是次优解。那么,有没有办法让第k个残差与正交,方法是有的,这就是下面要谈到的OMP算法。

    3、OMP算法

    3.1 算法描述

    OMP算法的改进之处在于:在分解的每一步对所选择的全部原子进行正交化处理,这使得在精度要求相同的情况下,OMP算法的收敛速度更快。

    那么在每一步中如何对所选择的全部原子进行正交化处理呢?在正式描述OMP算法前,先看一点基础思想。

    先看一个 k  阶模型,表示信号 f 经过 k 步分解后的情况,似乎很眼熟,但要注意它与MP算法不同之处,它的残值与前面每个分量正交,这就是为什么这个算法多了一个正交的原因,MP中仅与最近选出的的那一项正交。

    (1)

     k + 1 阶模型如下:

    (2)

    应用 k + 1阶模型减去阶模型,得到如下:

    (3)

    我们知道,字典矩阵D的原子是非正交的,引入一个辅助模型,它是表示对前k个项的依赖,描述如下:

    (4)

    和前面描述类似,span(x1, ...xk)之一上的正交投影操作,后面的项是残值。这个关系用数学符号描述:

    请注意,这里的 a  b 的上标表示第 k 步时的取值。

    (4)带入(3)中,有:

    5

    如果一下两个式子成立,(5)必然成立。

    (6)

    (7)

    ,有

    其中

    ak的值是由求法很简单,通过对(7)左右两边添加作内积消减得到:

    后边的第二项因为它们正交,所以为0,所以可以得出ak的第一部分。对于,在(4)左右两边中与作内积,可以得到ak的第二部分。

    对于(4),可以求出,求的步骤请参见参考文件的计算细节部分。为什么这里不提,因为后面会介绍更简单的方法来计算。

    3.2 收敛性证明

    通过(7),由于正交,将两个残值移到右边后求二范的平方,并将ak的值代入可以得到:

    可见每一次残差比上一次残差小,可见是收敛的。

    3.3 算法步骤

    整个OMP算法的步骤如下:

    由于有了上面的来龙去脉,这个算法就相当好理解了。

    到这里还不算完,后来OMP的迭代运算用另外一种方法可以计算得知,有位同学的论文[2]描述就非常好,我就直接引用进来:

    对比中英文描述,本质都是一样,只是有细微的差别。这里顺便贴出网一哥们写的OMP算法的代码,源出处不得而知,共享给大家。

    二、OMP的MATLAB实现

     1、一维信号重建

    代码:

    复制代码
    %  1-D信号压缩传感的实现(正交匹配追踪法Orthogonal Matching Pursuit)
    %  测量数M>=K*log(N/K),K是稀疏度,N信号长度,可以近乎完全重构
    %  编程人--香港大学电子工程系 沙威  Email: wsha@eee.hku.hk
    %  编程时间:2008年11月18日
    %  文档下载: http://www.eee.hku.hk/~wsha/Freecode/freecode.htm 
    %  参考文献:Joel A. Tropp and Anna C. Gilbert 
    %  Signal Recovery From Random Measurements Via Orthogonal Matching
    %  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
    %  DECEMBER 2007.
    
    clc;clear
    
    %%  1. 时域测试信号生成
    K=7;      %  稀疏度(做FFT可以看出来)
    N=256;    %  信号长度
    M=64;     %  测量数(M>=K*log(N/K),至少40,但有出错的概率)
    f1=50;    %  信号频率1
    f2=100;   %  信号频率2
    f3=200;   %  信号频率3
    f4=400;   %  信号频率4
    fs=800;   %  采样频率
    ts=1/fs;  %  采样间隔
    Ts=1:N;   %  采样序列
    x=0.3*cos(2*pi*f1*Ts*ts)+0.6*cos(2*pi*f2*Ts*ts)+0.1*cos(2*pi*f3*Ts*ts)+0.9*cos(2*pi*f4*Ts*ts);  %  完整信号,由4个信号叠加而来
    
    %%  2.  时域信号压缩传感
    Phi=randn(M,N);                                   %  测量矩阵(高斯分布白噪声)64*256的扁矩阵,Phi也就是文中说的D矩阵
    s=Phi*x.';                                        %  获得线性测量 ,s相当于文中的y矩阵
    
    %%  3.  正交匹配追踪法重构信号(本质上是L_1范数最优化问题)
    %匹配追踪:找到一个其标记看上去与收集到的数据相关的小波;在数据中去除这个标记的所有印迹;不断重复直到我们能用小波标记“解释”收集到的所有数据。
    
    m=2*K;                                            %  算法迭代次数(m>=K),设x是K-sparse的
    Psi=fft(eye(N,N))/sqrt(N);                        %  傅里叶正变换矩阵
    T=Phi*Psi';                                       %  恢复矩阵(测量矩阵*正交反变换矩阵)
    
    hat_y=zeros(1,N);                                 %  待重构的谱域(变换域)向量                     
    Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
    r_n=s;                                            %  残差值
    
    for times=1:m;                                    %  迭代次数(有噪声的情况下,该迭代次数为K)
        for col=1:N;                                  %  恢复矩阵的所有列向量
            product(col)=abs(T(:,col)'*r_n);          %  恢复矩阵的列向量和残差的投影系数(内积值) 
        end
        [val,pos]=max(product);                       %  最大投影系数对应的位置,即找到一个其标记看上去与收集到的数据相关的小波
        Aug_t=[Aug_t,T(:,pos)];                       %  矩阵扩充    
        
        T(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零),在数据中去除这个标记的所有印迹
        aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使残差最小
        r_n=s-Aug_t*aug_y;                            %  残差
        pos_array(times)=pos;                         %  纪录最大投影系数的位置
    end
    hat_y(pos_array)=aug_y;                           %  重构的谱域向量
    hat_x=real(Psi'*hat_y.');                         %  做逆傅里叶变换重构得到时域信号
    
    %%  4.  恢复信号和原始信号对比
    figure(1);
    hold on;
    plot(hat_x,'k.-')                                 %  重建信号
    plot(x,'r')                                       %  原始信号
    legend('Recovery','Original')
    norm(hat_x.'-x)/norm(x)                           %  重构误差
    复制代码

    结果:

    2、二维信号重建

    代码:

    复制代码
    %  本程序实现图像LENA的压缩传感
    %  程序作者:沙威,香港大学电气电子工程学系,wsha@eee.hku.hk
    %  算法采用正交匹配法,参考文献 Joel A. Tropp and Anna C. Gilbert 
    %  Signal Recovery From Random Measurements Via Orthogonal Matching
    %  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
    %  DECEMBER 2007.
    %  该程序没有经过任何优化
    
    function Wavelet_OMP
    
    clc;clear
    
    %  读文件
    X=imread('lena256.bmp');
    X=double(X);
    [a,b]=size(X);
    
    %  小波变换矩阵生成
    ww=DWT(a);
    
    %  小波变换让图像稀疏化(注意该步骤会耗费时间,但是会增大稀疏度)
    X1=ww*sparse(X)*ww';
    % X1=X;
    X1=full(X1);
    
    %  随机矩阵生成
    M=190;
    R=randn(M,a);
    % R=mapminmax(R,0,255);
    % R=round(R);
    
    %  测量值
    Y=R*X1;
    
    %  OMP算法
    %  恢复矩阵
    X2=zeros(a,b); 
    %  按列循环
    for i=1:b 
        %  通过OMP,返回每一列信号对应的恢复值(小波域)
        rec=omp(Y(:,i),R,a);
        %  恢复值矩阵,用于反变换
        X2(:,i)=rec;
    end
    
    
    %  原始图像
    figure(1);
    imshow(uint8(X));
    title('原始图像');
    
    %  变换图像
    figure(2);
    imshow(uint8(X1));
    title('小波变换后的图像');
    
    %  压缩传感恢复的图像
    figure(3);
    %  小波反变换
    X3=ww'*sparse(X2)*ww; 
    % X3=X2;
    X3=full(X3);
    imshow(uint8(X3));
    title('恢复的图像');
    
    %  误差(PSNR)
    %  MSE误差
    errorx=sum(sum(abs(X3-X).^2));        
    %  PSNR
    psnr=10*log10(255*255/(errorx/a/b))   
    
    
    %  OMP的函数
    %  s-测量;T-观测矩阵;N-向量大小
    function hat_y=omp(s,T,N)
    Size=size(T);                                     %  观测矩阵大小
    M=Size(1);                                        %  测量
    hat_y=zeros(1,N);                                 %  待重构的谱域(变换域)向量                     
    Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
    r_n=s;                                            %  残差值
    
    for times=1:M;                                  %  迭代次数(稀疏度是测量的1/4)
        for col=1:N;                                  %  恢复矩阵的所有列向量
            product(col)=abs(T(:,col)'*r_n);          %  恢复矩阵的列向量和残差的投影系数(内积值) 
        end
        [val,pos]=max(product);                       %  最大投影系数对应的位置
        Aug_t=[Aug_t,T(:,pos)];                       %  矩阵扩充
        T(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零)
        aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使残差最小
        r_n=s-Aug_t*aug_y;                            %  残差
        pos_array(times)=pos;                         %  纪录最大投影系数的位置
        
        if (norm(r_n)<0.9)                            %  残差足够小
            break;
        end
    end
    hat_y(pos_array)=aug_y;                           %  重构的向量
    复制代码

    结果:

    三、OMP算法中的数学知识

    算法会有一些线性代数中的概念和知识,主要是关于子空间、最小二乘、投影矩阵等,详情请参考

    (数学)最小二乘的几何意义及投影矩阵

    主要内容:

    1. 什么是最小二乘
    2. 最小二乘的几何意义
    3. 正交投影矩阵

     

    什么是最小二乘?

    假设我们手上有n组成对的数据,{(xi,yi):i=1…n},为了探究y变量与x变量的关系,我们希望用一个多项式来匹配它,可是多项式中的系数怎么确定呢?拿来拼凑肯定是不行的,最小二乘法告诉我们,这个多项式的系数应该让每个点的误差的平方之和最小。

    (百度百科)最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

    最小二乘的几何意义

    最小二乘的几何意义:最小二乘法中的几何意义是高维空间中的一个向量在低维子空间的投影。

    从上面的定义中,我们很难想象到最小二乘的几何意义,那么我们通过一个简单的例子来推导一下:

    我们根据定义中的误差平方之和最小化来拟合直线:

    每个点的误差表示:

    最小误差的平方和:

    要求解上面的最小化问题,我们可以通过求导的方式得到,最好是转化为矩阵表达形式:AX=b (这里x表示上述的系数a)

    求得结果为:

    如果通过超定方程的解法,很容易就可以得到上面结果。

    先来说说向量表达形式:

    小括号中表示:它是两个向量 [1, ... , 1]T 和 [x1, ... , xn]T 的线性组合,换句话说,它是这两个向量构成的二维子空间(想成一个平面就可以)的任意一点。

    那么上面式子的几何含义:表示向量 [y1, ... , yn]T(表示空间中的一点) 到这个二维子空间任意一点的距离;(向量的长度)

    最小化上面式子的平方(向量长度的最小化)的几何含义:寻找在 [1, ... , 1]T 和 [x1, ... , xn]T 构成的二维子空间上的一个点,使得向量 [y1, ... , yn]T 到这个点的距离最小。怎么找这个点呢?只要做一个几何投影就好了。(如下图)

    如上图所示,在三维空间中给定一个向量 u,以及由向量 v1,v2 构成的一个二维平面,向量 p 为 u 到这个平面的投影,它是 v1,v2 的线性组合:

    利用投影的垂直性质,我们可以得到关于系数C的两个方程:

    令 V = [v1, v2], p = c1v1 + c2v2,将上述式子合并并转化为矩阵形式(更容易扩展到高维空间),得到:

     

    因此系数c的表达式为:

    有没有发现很熟悉?和式子 一模一样有木有!!!

    好了,我们回到原来的例子,看看几何关系中的投影点和被投影的空间分别代表什么。

    把图中的 u 替换成 [y1, ... , yn]T ,把 v1,v2 分别替换成 [1, ... , 1]T 和 [x1, ... , xn]T, 系数 c1 和 c2 也就是我们要求的 a0,a1。

    所以,最小二乘法的几何意义是高维空间的一个向量(由y数据决定)在低维子空间(由x数据以及多项式的次数决定)的投影。

    正交投影矩阵

    上面提到了最小二乘的几何意义就是空间中的投影,其实投影在线性代数中也是存在其数学公式的,可以联系以下数学知识来理解最小二乘的几何意义。

    张成子空间:

    张成子空间的投影矩阵:

    最小二乘的投影解释:


    展开全文
  •  正交匹配追踪算法的本质思想是,以贪婪迭代的方式选择测量矩阵Φ的列,使得在每次迭代中所选择的列与当前的冗余向量最大程度地相关,从测量向量中减去相关部分并反复迭代,直到迭代次数达到稀疏度K,强制迭代停止...

    压缩感知之OMP恢复算法

    1、基本思想

      y=Φx   x=Ψθ
      正交匹配追踪算法的本质思想是,以贪婪迭代的方式选择测量矩阵Φ的列,使得在每次迭代中所选择的列与当前的冗余向量最大程度地相关,从测量向量中减去相关部分并反复迭代,直到迭代次数达到稀疏度K,强制迭代停止。

    2、算法步骤

      输入:(1)M*N的感知矩阵A,其中M远远小于N,A=Φ*Ψ。
         (2)长度为M的数据向量b,即测量值y。
      输出:长度为N的重建向量xˆ,满足y=Ax。
      初始化:残差r0=y,重建信号x0=0,索引集Λ0=Φ,迭代次数n=2*K,计数器k=0。
      步骤1:计算残差和感知矩阵A的每一列的投影系数(内积值)ck=ATrk1
      步骤2:找出ck中元素最大的元素ck=max{ck}以及对应的位置pos;
      步骤3:更新索引集 Λk=Λk1{pos},以及原子集合AΛK=AΛk1{A(:pos)}
      步骤4:利用最小二乘求得近似解 xk=(ATΛkAΛk1)1ATΛky
      步骤5:更新余量rk=yAxk
      步骤6:判断迭代是否满足停止条件,满足则停止xˆ=xk,r=rk, 输出xˆ,r ,否则转步骤1。

    3、仿真验证

      3.1 一维时间稀疏信号
      首先进行一维时间稀疏信号的恢复,信号长度为512,稀疏度选取10、20、30、40、50,matlab代码如下:

    clc;clear;close all
    %%  1. 时域测试信号生成
    CNT = 100;                                       %对于每组(K,M,N),重复迭代次数  
    N=512;                                           %信号长度
    K_set= [10,20,30,40,50];                             %信号x的稀疏度集合 
    Percentage = zeros(length(K_set),N);             %存储恢复成功概率
    for kk=1:length(K_set)
        K=K_set(kk);                                      %本次稀疏度
        M_set=1:5:N;                                          %测量数每隔五个取一次
        PercentageK = zeros(1,length(M_set));                 %存储此稀疏度K下不同M的恢复成功概率  
        for mm=1:length(M_set)
            M=M_set(mm);                                      %本次观测次数
            P=0;
            for cnt=1:CNT                             %每个观测值个数均运行CNT次 
                Index_K=randperm(N);                              %将1-N随机打乱 行向量
                x=zeros(N,1);
                x(Index_K(1:K))=5*randn(K,1);                     %x为K稀疏的,且位置是随机的
    
    %%  2.  时域信号压缩传感
                Phi=randn(M,N);                                   %  测量矩阵(高斯分布白噪声)
                Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(M),1]); %正则化
                y=Phi*x;                                        %  获得线性测量 
    
    %%  3.  正交匹配追踪法重构信号(本质上是L_1范数最优化问题)
    n=2*K;                                            %  算法迭代次数(m>=K)
    Psi=eye(N);                                       %  单位矩阵为正变换矩阵
    A=Phi*Psi;                                        %  恢复矩阵(测量矩阵*正交反变换矩阵)
    
    hat_y=zeros(1,N);                                 %  待重构的变换域里的向量                     
    Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
    r0=y;                                             %  残差值
    
    for times=1:n;                                    %  迭代次数(有噪声的情况下,该迭代次数为K)
        for col=1:N;                                  %  恢复矩阵的所有列向量    步骤1
            product(col)=abs(A(:,col)'*r0);           %  恢复矩阵的列向量和残差的投影系数(内积值) 
        end
        [val,pos]=max(product);                       %  最大投影系数对应的位置  步骤2
        Aug_t=[Aug_t,A(:,pos)];                       %  矩阵扩充               步骤3 更新原子集合
        A(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零)
        aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*y;           %  最小二乘,使残差最小     步骤4 求近似解
        r0=y-Aug_t*aug_y;                             %  残差                   步骤5 更新余量
        pos_array(times)=pos;                         %  纪录最大投影系数的位置  步骤3 更新索引集
    end
    hat_y(pos_array)=aug_y;                           %  重构的变换域里的向量
    hat_x=real(Psi'*hat_y.');                         %  做逆变换重构得到时域信号
                if norm(hat_x-x)<1e-6                   %如果残差小于1e-6则认为恢复成功  
                    P = P + 1;  
                end
            end
            PercentageK(mm) = P/CNT;                  %计算恢复概率
        end
            Percentage(kk,1:length(M_set)) = PercentageK;
    end
    save MtoPercentage100 %运行一次不容易,把变量全部存储下来
    %}
    %% 绘制概率图
    S = ['-ks';'-ko';'-kd';'-kv';'-k*'];  
    figure;  
    for kk = 1:length(K_set)  
        K = K_set(kk);  
        M_set = 1:5:N;  
        L_Mset = length(M_set);  
        plot(M_set,Percentage(kk,1:L_Mset),S(kk,:));%绘出x的恢复信号  
        hold on;  
    end  
    hold off;  
    xlim([0 256]);  
    legend('K=10','K=20','K=30','K=40','K=50');  
    xlabel('Number of measurements(M)');  
    ylabel('Percentage recovered');  
    title('Percentage of input signals recovered correctly(N=256)(Gaussian)');  
    
    %%  4.  恢复信号和原始信号对比
    % figure(1);
    % hold on;
    % plot(hat_x,'k.-')                                 %  重建信号
    % plot(x,'r')                                       %  原始信号
    % legend('Recovery','Original')
    % norm(hat_x.'-x)/norm(x)                           %  重构误差

    程序代码改变测量数得到恢复概率,如下图所示


    《压缩感知及应用》书上的p68页图3-7如下
    这里写图片描述
    结果要比《压缩感知及应用》书上p68页结果感觉略好,原因不详。
    3.2 二维lena图像恢复
    仿照《压缩感知及应用》这本书,p120页,进行lena图像恢复,压缩比使用0.3、0.4、0.5。所程序如下:

    img=imread('lena256.bmp');      %读文件
    img=double(img);             
    [n,b]=size(img);             %文件为n行b列
    figure(1)
    subplot(2,2,1)
    imshow(img,[])
    weizhi=1;
    Pm=zeros(3,1);           %放功率信噪比
    Tm=zeros(3,1);           %放时间
    
    for CR=0.3:0.1:0.5
        disp(CR);
        %  测量矩阵
        m=floor(n*CR);                              
        Phi=randn(m,n);                     %  测量矩阵生成
        for i=1:n                           %  测量矩阵归一化
            Phi(:,i) = Phi(:,i) / norm(Phi(:,i));%正则化测量矩阵φ
        end 
    
        %  对图像进行欠采样
        y=Phi*img;       
    
        %% CS重建( 已知测量值y,测量矩阵Phi,稀疏基Psi'(小波变换矩阵) )
        Psi=DWT(n);                         %  小波变换矩阵生成(要求大小是2的整数幂次)%傅里叶正变换矩阵dftmtx(n)
        A = Phi*Psi';                       %  y = Phi*X0 = Phi*Psi'*s(s是稀疏系数)此处X0=Psi's
        y = y*Psi';                         %  此处不明白,为什么还要乘psi' 如果不乘,恢复效果很差
    
        %OMP算法
        tic
        fprintf('OMP...\r\n');
        ReS3 = zeros(n,b);
        for i = 1:b                        
            rec = omp(y(:,i),A,n);           %对y的的每一列进行重构,恢复变换域矩阵
            ReS3(:,i) = rec;
        end
        T3 = toc;
        ReS3 = Psi'*sparse(ReS3)*Psi;       %%%%%%%%%%%%%%%         
        ReImg3 = full(ReS3);
    
        weizhi=weizhi+1;
        subplot(2,2,weizhi)
        imshow(ReImg3,[]);
        %% 计算峰值噪声(PSNR)、用时
    
        %   OMP
        errorx=sum(sum(abs(ReImg3-img).^2));      %  MSE误差
        psnr=10*log10(255*255/(errorx/n/b));%  
        Pm(weizhi-1,1)=psnr;
        Tm(weizhi-1,1)=T3;
    
    end
    %% 显示结果
    
    CR=0.3:0.1:0.5;
    figure;
    plot(CR,Pm(:,1),'ro');
    hold on
    plot(CR,Pm(:,1),'r');
    xlabel('压缩比');
    ylabel('峰值信噪比/dB');
    axis([0.2,0.6,5,30]);
    
    figure;
    plot(CR,Tm(:,1),'g*')
    hold on
    plot(CR,Tm(:,1),'g')
    
    xlabel('压缩比');
    ylabel('运行时间/s');
    axis([0.2,0.6,0,1+max(Tm(3,:))]);
    
    
    function hat_y=omp(s,T,N)
    %  OMP的函数
    %  s-测量;T-观测矩阵;N-向量大小
    
    Size=size(T);                                     %  观测矩阵大小
    M=Size(1);                                        %  测量
    hat_y=zeros(1,N);                                 %  待重构的谱域(变换域)向量                     
    Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
    r_n=s;                                            %  残差值
    
    for times=1:M;                                    %  迭代次数(不会超过测量值)
    
        for col=1:N;                                  %  恢复矩阵的所有列向量
            product(col)=abs(T(:,col)'*r_n);          %  恢复矩阵的列向量和残差的投影系数(内积值) 
        end
        [val,pos]=max(product);                       %  最大投影系数对应的位置
        Aug_t=[Aug_t,T(:,pos)];                       %  矩阵扩充
        T(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零)
        aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使残差最小
        r_n=s-Aug_t*aug_y;                            %  残差
        pos_array(times)=pos;                         %  纪录最大投影系数的位置
    
        if (abs(aug_y(end))^2/norm(aug_y)<0.0005)       %  自适应截断误差(***需要调整经验值)
            break;
        end
    end
    
    hat_y(pos_array)=aug_y;                           %  重构的向量
    
    %  程序作者:沙威,香港大学电气电子工程学系,wsha@eee.hku.hk
    %  参考文献:小波分析理论与MATLAB R2007实现,葛哲学,沙威,第20章  小波变换在矩阵方程求解中的应用(沙威、陈明生编写).
    
    %  构造正交小波变换矩阵,图像大小N*N,N=2^P,P是整数。
    
    function ww=DWT(N)
    
    [h,g]= wfilters('sym8','d');       %  获得symlets8小波的低通分解滤波器和高通分解滤波器的系数;
                                       %  采用Symlets8的小波分解可得到较好的压缩比率及较高的信号恢复质量
    
    % N=256;                           %  矩阵维数(大小为2的整数幂次)
    L=length(h);                       %  滤波器长度
    rank_max=log2(N);                  %  最大层数 8 ?
    rank_min=double(int8(log2(L)))+1;  %  最小层数 5 ?
    ww=1;   %  预处理矩阵
    
    %  矩阵构造
    for jj=rank_min:rank_max
    
        nn=2^jj;
    
        %  构造向量
        p1_0=sparse([h,zeros(1,nn-L)]);
        p2_0=sparse([g,zeros(1,nn-L)]);
    
        %  向量圆周移位
        for ii=1:nn/2
            p1(ii,:)=circshift(p1_0',2*(ii-1))'; %对1*m矩阵来说,相当于右移
            p2(ii,:)=circshift(p2_0',2*(ii-1))';
        end
    
        %  构造正交矩阵
        w1=[p1;p2];
        mm=2^rank_max-length(w1);
        w=sparse([w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)]);
        ww=ww*w;                     % ?
    
        clear p1;clear p2;
    end

    结果如下,第一幅为原图,然后依次为压缩比0.3、0.4、0.5,
    这里写图片描述
    psnr和所用时间分别如下
    这里写图片描述
    这里写图片描述
    书上的图像如下,依次为原图,压缩比为0.5、0.4、0.3
    这里写图片描述
    所得实际结果感觉较之稍差。

    4、结果分析

      由于使用的测量矩阵为高斯随机矩阵, 所以每次结果会有偏差,这是部分原因,其余应该还有其它原因。

    5、参考文献

      1、《压缩感知及应用》,闫敬文 刘蕾 屈小波
      2、沙威老师压缩感知代码
      3、彬彬有礼的博客,地址如下
      http://blog.csdn.net/jbb0523/article/details/45130793

    展开全文
  • 在分析分段正交匹配追踪StOMP算法迭代过程中多原子匹配方法的基础上,为了进一步减少算法迭代次数,提高重构精度,提出了基于互相关向量的自适应极限因子选取和迭代结束条件重设的优化方法。通过MATLAB编程实验,在...
  • 这是本人编写的基础类的压缩感知重构算法,主要是匹配追踪算法
  • 多路径匹配追踪深度优先(MMP-DF),包含深度优先程序代码,方便理解该算法的原子搜索方式
  • 压缩感知的稀疏重构中广泛应用的正交匹配追踪(OMP)算法matlab程序,该算法由香港大学电子工程系 沙威老师开发,代码注释详细,便于读者理解。已测试,可以正常运行。读者通过代码可以加深对该算法以及压缩感知、...

空空如也

空空如也

1 2 3 4 5 ... 8
收藏数 158
精华内容 63
关键字:

matlab匹配追踪算法

matlab 订阅