精华内容
参与话题
问答
  • 压缩感知重构算法之正交匹配追踪(OMP)

    万次阅读 多人点赞 2015-04-19 17:26:12
    题目:压缩感知重构算法之正交匹配追踪(OMP)  前面经过几篇的基础铺垫,本篇给出正交匹配追踪(OMP)算法的MATLAB函数代码,并且给出单次测试例程代码、测量数M与重构成功概率关系曲线绘制例程代码、信号稀疏度K与...

    题目:压缩感知重构算法之正交匹配追踪(OMP)

            前面经过几篇的基础铺垫,本篇给出正交匹配追踪(OMP)算法的MATLAB函数代码,并且给出单次测试例程代码、测量数M与重构成功概率关系曲线绘制例程代码、信号稀疏度K与重构成功概率关系曲线绘制例程代码。

    0、符号说明如下:

            压缩观测y=Φx,其中y为观测所得向量M×1,x为原信号N×1(M<<N)。x一般不是稀疏的,但在某个变换域Ψ是稀疏的,即x=Ψθ,其中θ为K稀疏的,即θ只有K个非零项。此时y=ΦΨθ,令A=ΦΨ,则y=

            (1) y为观测所得向量,大小为M×1

            (2)x为原信号,大小为N×1

            (3)θ为K稀疏的,是信号在x在某变换域的稀疏表示

            (4)Φ称为观测矩阵、测量矩阵、测量基,大小为M×N

            (5)Ψ称为变换矩阵、变换基、稀疏矩阵、稀疏基、正交基字典矩阵,大小为N×N

            (6)A称为测度矩阵、传感矩阵、CS信息算子,大小为M×N

    上式中,一般有K<<M<<N,后面三个矩阵各个文献的叫法不一,以后我将Φ称为测量矩阵、将Ψ称为稀疏矩阵、将A称为传感矩阵

    1、OMP重构算法流程:



    2、正交匹配追踪(OMP)MATLAB代码(CS_OMP.m)

    function [ theta ] = CS_OMP( y,A,t )
    %CS_OMP Summary of this function goes here
    %Version: 1.0 written by jbb0523 @2015-04-18
    %   Detailed explanation goes here
    %   y = Phi * x
    %   x = Psi * theta
    %	y = Phi*Psi * theta
    %   令 A = Phi*Psi, 则y=A*theta
    %   现在已知y和A,求theta
        [y_rows,y_columns] = size(y);
        if y_rows<y_columns
            y = y';%y should be a column vector
        end
        [M,N] = size(A);%传感矩阵A为M*N矩阵
        theta = zeros(N,1);%用来存储恢复的theta(列向量)
        At = zeros(M,t);%用来迭代过程中存储A被选择的列
        Pos_theta = zeros(1,t);%用来迭代过程中存储A被选择的列序号
        r_n = y;%初始化残差(residual)为y
        for ii=1:t%迭代t次,t为输入参数
            product = A'*r_n;%传感矩阵A各列与残差的内积
            [val,pos] = max(abs(product));%找到最大内积绝对值,即与残差最相关的列
            At(:,ii) = A(:,pos);%存储这一列
            Pos_theta(ii) = pos;%存储这一列的序号
            A(:,pos) = zeros(M,1);%清零A的这一列,其实此行可以不要,因为它与残差正交
            %y=At(:,1:ii)*theta,以下求theta的最小二乘解(Least Square)
            theta_ls = (At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;%最小二乘解
            %At(:,1:ii)*theta_ls是y在At(:,1:ii)列空间上的正交投影
            r_n = y - At(:,1:ii)*theta_ls;%更新残差        
        end
        theta(Pos_theta)=theta_ls;%恢复出的theta
    end


    3、OMP单次重构测试代码(CS_Reconstuction_Test.m)

            代码中,直接构造一个K稀疏的信号,所以稀疏矩阵为单位阵。

    %压缩感知重构算法测试
    clear all;close all;clc;
    M = 64;%观测值个数
    N = 256;%信号x的长度
    K = 10;%信号x的稀疏度
    Index_K = randperm(N);
    x = zeros(N,1);
    x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的
    Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta
    Phi = randn(M,N);%测量矩阵为高斯矩阵
    A = Phi * Psi;%传感矩阵
    y = Phi * x;%得到观测向量y
    %% 恢复重构信号x
    tic
    theta = CS_OMP(y,A,K);
    x_r = Psi * theta;% x=Psi * theta
    toc
    %% 绘图
    figure;
    plot(x_r,'k.-');%绘出x的恢复信号
    hold on;
    plot(x,'r');%绘出原信号x
    hold off;
    legend('Recovery','Original')
    fprintf('\n恢复残差:');
    norm(x_r-x)%恢复残差

    运行结果如下:(信号为随机生成,所以每次结果均不一样)

    1)图:


    2)Command Windows

    Elapsed time is 0.849710 seconds.
    恢复残差:
    ans =
      5.5020e-015

    4、测量数M与重构成功概率关系曲线绘制例程代码

    %压缩感知重构算法测试CS_Reconstuction_MtoPercentage.m
    %   绘制参考文献中的Fig.1
    %   参考文献: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.
    %   Elapsed time is 1171.606254 seconds.(@20150418night)
    clear all;close all;clc;
    %% 参数配置初始化
    CNT = 1000;%对于每组(K,M,N),重复迭代次数
    N = 256;%信号x的长度
    Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta
    K_set = [4,12,20,28,36];%信号x的稀疏度集合
    Percentage = zeros(length(K_set),N);%存储恢复成功概率
    %% 主循环,遍历每组(K,M,N)
    tic
    for kk = 1:length(K_set)
        K = K_set(kk);%本次稀疏度
        M_set = K:5:N;%M没必要全部遍历,每隔5测试一个就可以了
        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);
                x = zeros(N,1);
                x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的                
                Phi = randn(M,N);%测量矩阵为高斯矩阵
                A = Phi * Psi;%传感矩阵
                y = Phi * x;%得到观测向量y
                theta = CS_OMP(y,A,K);%恢复重构信号theta
                x_r = Psi * theta;% x=Psi * theta
                if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功
                    P = P + 1;
                end
           end
           PercentageK(mm) = P/CNT*100;%计算恢复概率
        end
        Percentage(kk,1:length(M_set)) = PercentageK;
    end
    toc
    save MtoPercentage1000 %运行一次不容易,把变量全部存储下来
    %% 绘图
    S = ['-ks';'-ko';'-kd';'-kv';'-k*'];
    figure;
    for kk = 1:length(K_set)
        K = K_set(kk);
        M_set = K: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=4','K=12','K=20','K=28','K=36');
    xlabel('Number of measurements(M)');
    ylabel('Percentage recovered');
    title('Percentage of input signals recovered correctly(N=256)(Gaussian)');

            本程序在联想ThinkPadE430C笔记本(4GB DDR3内存,i5-3210)上运行共耗时1171.606254秒,程序中将所有数据均通过“save MtoPercentage1000”存储了下来,以后可以再对数据进行分析,只需“load MtoPercentage1000”即可。

            程序运行结果比文献[1]的Fig.1要好,原因不详。

    本程序运行结果:


    文献[1]中的Fig.1:


    5、信号稀疏度K与重构成功概率关系曲线绘制例程代码

    %压缩感知重构算法测试CS_Reconstuction_KtoPercentage.m
    %   绘制参考文献中的Fig.2
    %   参考文献: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.
    %   Elapsed time is 1448.966882 seconds.(@20150418night)
    clear all;close all;clc;
    %% 参数配置初始化
    CNT = 1000;%对于每组(K,M,N),重复迭代次数
    N = 256;%信号x的长度
    Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta
    M_set = [52,100,148,196,244];%测量值集合
    Percentage = zeros(length(M_set),N);%存储恢复成功概率
    %% 主循环,遍历每组(K,M,N)
    tic
    for mm = 1:length(M_set)
        M = M_set(mm);%本次测量值个数
        K_set = 1:5:ceil(M/2);%信号x的稀疏度K没必要全部遍历,每隔5测试一个就可以了
        PercentageM = zeros(1,length(K_set));%存储此测量值M下不同K的恢复成功概率
        for kk = 1:length(K_set)
           K = K_set(kk);%本次信号x的稀疏度K
           P = 0;
           for cnt = 1:CNT %每个观测值个数均运行CNT次
                Index_K = randperm(N);
                x = zeros(N,1);
                x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的                
                Phi = randn(M,N);%测量矩阵为高斯矩阵
                A = Phi * Psi;%传感矩阵
                y = Phi * x;%得到观测向量y
                theta = CS_OMP(y,A,K);%恢复重构信号theta
                x_r = Psi * theta;% x=Psi * theta
                if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功
                    P = P + 1;
                end
           end
           PercentageM(kk) = P/CNT*100;%计算恢复概率
        end
        Percentage(mm,1:length(K_set)) = PercentageM;
    end
    toc
    save KtoPercentage1000test %运行一次不容易,把变量全部存储下来
    %% 绘图
    S = ['-ks';'-ko';'-kd';'-kv';'-k*'];
    figure;
    for mm = 1:length(M_set)
        M = M_set(mm);
        K_set = 1:5:ceil(M/2);
        L_Kset = length(K_set);
        plot(K_set,Percentage(mm,1:L_Kset),S(mm,:));%绘出x的恢复信号
        hold on;
    end
    hold off;
    xlim([0 125]);
    legend('M=52','M=100','M=148','M=196','M=244');
    xlabel('Sparsity level(K)');
    ylabel('Percentage recovered');
    title('Percentage of input signals recovered correctly(N=256)(Gaussian)');

            本程序在联想ThinkPadE430C笔记本(4GB DDR3内存,i5-3210)上运行共耗时1448.966882秒,程序中将所有数据均通过“save KtoPercentage1000”存储了下来,以后可以再对数据进行分析,只需“load KtoPercentage1000”即可。

            程序运行结果比文献[1]的Fig.2要好,原因不详。

    本程序运行结果:


    文献[1]中的Fig.2:


    6、参考文献

    【1】Joel A. Tropp and Anna C. Gilbert. Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit[J]. IEEETransactions on Information Theory, VOL. 53, NO. 12, DECEMBER 2007.

    【2】Y.C.Pati, R.Rezaiifar,and P.S.Krishnaprasad. Orthogonal Matching Pursuit-Recursive FunctionApproximation with Applications to wavelet decomposition, Proc. 27thAnnu. Asilomar Conf. Signals, Systems, and Computers, Pacific Grove, CA, Nov.1993,vol.1,pp40-44.

    【3】沙威.CS_OMP,http://www.eee.hku.hk/~wsha/Freecode/Files/CS_OMP.zip

    展开全文
  • MP算法和OMP算法及其思想

    万次阅读 多人点赞 2012-04-17 03:09:18
    主要介绍MP(Matching Pursuits)算法OMP(Orthogonal Matching Pursuit)算法[1],这两个算法虽然在90年代初就提出来了,但作为经典的算法,国内文献(可能有我没有搜索到)都仅描述了算法步骤和简单的应用,并未对其...

    主要介绍MP(Matching Pursuits)算法和OMP(Orthogonal Matching Pursuit)算法[1],这两个算法虽然在90年代初就提出来了,但作为经典的算法,国内文献(可能有我没有搜索到)都仅描述了算法步骤和简单的应用,并未对其进行详尽的分析,国外的文献还是分析的很透彻,所以我结合自己的理解,来分析一下写到博客里,算作笔记。

    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相同,而且这些向量已作为归一化处理,即|,也就是单位向量长度为1。MP算法的基本思想:从字典矩阵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算法迭代会发现总是在x1和x2上反复迭代,即,这个就是信号(残值)在已选择的原子进行垂直投影的非正交性导致的。再用严谨的方式描述[1]可能容易理解:在Hilbert空间H中,,定义,就是它是这些向量的张成中的一个,MP构造一种表达形式:;这里的Pvf表示 f在V上的一个正交投影操作,那么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算法的代码,源出处不得而知,共享给大家。

    再贴另外一个洋牛paper[3]中关于OMP的描述,之所以引入,是因为它描述的非常严谨,但是也有点苦涩难懂,不过有了上面的基础,就容易多了。


    它的描述中的Sweep步骤就是寻找与当前残差最大的内积时列在字典矩阵D中的索引,它的这个步骤描述说明为什么要选择内积最大的以及如何选择。见下图,说的非常清晰。


    它的算法步骤Update Provisional Solution中求很简单,就是在 b = Ax 已知 A和b求x, 在x的最小二范就是A的伪逆与b相乘,即:


    看上去头疼,其实用matlab非常简单,看看上面的matlab的代码就明白了。

    我们可以看得出来,算法流程清晰明了,还是很好理解的。这正是OMP算法的魅力,作为工具使用简单,背后却隐藏着很有趣的思想。

    写这篇博客的目的,是因为搜索了一下,MP和OMP没有人很详细的介绍。文献[1]讲的很清楚的,大家有兴趣可以找来看看。不要被老板发现我居然在搜中文文献还写中文博客。


    参考文献:

    [1] Orthogonal Matching Pursuit:Recursive Function Approximat ion with Applications to Wavelet Decomposition
    [2]http://wenku.baidu.com/view/22f3171614791711cc7917e4.html

    [3] From Sparse Solutions of Systems of Equations to Sparse Modeling of Signals and Images

    展开全文
  • OMP算法

    千次阅读 2016-12-07 09:49:54
    OMP算法思想假设信号s∈Rd\mathbf{s}\in\mathbb{R}^d的稀疏度为mm, {x1,...,xN}\{\mathbf{x}_1,...,\mathbf{x}_N\}是NN个基向量. Φ∈RN×d\mathbf{\Phi}\in\mathbb{R}^{N\times d}是由这些基向量构成的测量矩阵。...

    OMP算法思想

    假设信号sRd的稀疏度为m, {x1,...,xN}N个基向量. ΦRN×d是由这些基向量构成的测量矩阵。通过测量矩阵,我们获得一个N维向量的测量值

    v=Φs

    OMP算法用于从测量值v中恢复信号s. 算法的思想是以贪婪的方式选取Φ中与v的剩余项最相关的列,减去这些列对v的贡献后,对剩余项迭代,我们希望在m次迭代后,算法可以选出正确的列来

    OMP算法表述

    OMP

    简单的几个例子

    例子一 稀疏信号的恢复

    -代码片

    %压缩感知重构算法测试  
    clear all;close all;clc;  
    M = 64;%观测值个数  
    N = 256;%信号x的长度  
    K = 10;%信号x的稀疏度  
    Index_K = randperm(N);  
    x = zeros(N,1);  
    x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的  
    Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta  
    Phi = randn(M,N);%测量矩阵为高斯矩阵  
    A = Phi * Psi;%传感矩阵  
    y = Phi * x;%得到观测向量y  
    %% 恢复重构信号x  
    tic  
    theta = CS_OMP(y,A,K);  
    x_r = Psi * theta;% x=Psi * theta  
    toc  
    %% 绘图  
    figure;  
    plot(x_r,'k.-');%绘出x的恢复信号  
    hold on;  
    plot(x,'r');%绘出原信号x  
    hold off;  
    legend('Recovery','Original')  
    fprintf('\n恢复残差:');  
    norm(x_r-x)%恢复残差  

    function [ theta ] = CS_OMP( y,A,t )  
    %CS_OMP Summary of this function goes here  
    %Version: 1.0 written by jbb0523 @2015-04-18  
    %   Detailed explanation goes here  
    %   y = Phi * x  
    %   x = Psi * theta  
    %   y = Phi*Psi * theta  
    %   令 A = Phi*Psi, 则y=A*theta  
    %   现在已知y和A,求theta  
        [y_rows,y_columns] = size(y);  
        if y_rows<y_columns  
            y = y';%y should be a column vector  
        end  
        [M,N] = size(A);%传感矩阵A为M*N矩阵  
        theta = zeros(N,1);%用来存储恢复的theta(列向量)  
        At = zeros(M,t);%用来迭代过程中存储A被选择的列  
        Pos_theta = zeros(1,t);%用来迭代过程中存储A被选择的列序号  
        r_n = y;%初始化残差(residual)为y  
        for ii=1:t%迭代t次,t为输入参数  
            product = A'*r_n;%传感矩阵A各列与残差的内积  
            [val,pos] = max(abs(product));%找到最大内积绝对值,即与残差最相关的列  
            At(:,ii) = A(:,pos);%存储这一列  
            Pos_theta(ii) = pos;%存储这一列的序号  
            A(:,pos) = zeros(M,1);%清零A的这一列,其实此行可以不要,因为它与残差正交  
            %y=At(:,1:ii)*theta,以下求theta的最小二乘解(Least Square)  
            theta_ls = (At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;%最小二乘解  
            %At(:,1:ii)*theta_ls是y在At(:,1:ii)列空间上的正交投影  
            r_n = y - At(:,1:ii)*theta_ls;%更新残差          
        end  
        theta(Pos_theta)=theta_ls;%恢复出的theta  
    end 

    -运行结果
    1) 图
    这里写图片描述
    2) 结果
    Elapsed time is 0.849710 seconds.
    恢复残差:
    ans =
    5.5020e-015

    ###例子二 时域信号压缩传感
    -代码片

    %  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)                           %  重构误差

    -运行结果
    这里写图片描述

    ##参考文献
    [1] Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition,” in Proceedings of 27th Asilomar Conference on Signals, Systems and Computers, 1993, pp. 40–44 vol.1.
    [2] J. A. Tropp and A. C. Gilbert, “Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit,” IEEE Transactions on Information Theory, vol. 53, no. 12, pp. 4655–4666, Dec. 2007.
    [3]OMP算法
    [4] OMP matlab
    [5] OMP free source code

    一点感想

    1. 如果你有时间、有条件、有氛围去学习和研究,请抓住机会,别等有一天你想搞点研究的时候发现你根本没有时间、没有条件、没有氛围……
    2. 这位博主的博客写得很不错
      彬彬有礼的博客
    展开全文
  • MP算法与OMP算法

    千次阅读 2014-01-17 13:40:01
    MP算法是一种贪心算法(greedy),每次迭代选取与当前样本残差最接近的原子,直至残差满足一定条件。 首先解决两个问题,怎么定义“最接近原子”,怎么计算残差? 选择最接近残差的原子:MP里定义用向量内积原子与...

    稀疏编码的一般最优化公式为:


    其中的零范数为非凸优化。那么如何解这么一个非凸优化问题呢?其中一个常用的解法就是MP算法。


    MP算法

    MP算法是一种贪心算法(greedy),每次迭代选取与当前样本残差最接近的原子,直至残差满足一定条件。


    求解方法


    首先解决两个问题,怎么定义“最接近原子”,怎么计算残差?

    选择最接近残差的原子:MP里定义用向量内积原子与残差的距离,我们用R表示残差,di表示原子,则:

    Max[Dist(R,di)]=max[<R,di>];

    残差更新:R=R-<R,di>I;继续选择下一个,直至收敛;


    需要注意的是,MP算法中要求字典原子||di||=1,上面的公式才成立。

    我们用二维空间上的向量来表示,用如下的图来表述上面的过程:


    上图中d1,d2,d3表示归一化的原子,红色向量r表示当前残差;

    进过内积计算,<r,d3>最大,于是r分解为d3方向以及垂直于d3方向的两个向量(<r,d3>d3及r-<r,d3>d3),把d3方向的分量(<r,d3>d3)加入到已经求得的重构项中,那么绿色向量(r-<r,d3>d3)变为新的残差。

    再一轮迭代得到如下:


    R往d1方向投影分解,绿色向量成为新的残差。

     

    具体算法:

     


    收敛性

    从上面的向量图我们可以清楚地看出,k+1的残差Rk+1是k步残差Rk的分量。根据直角三角形斜边大于直角边,|Rk+1|<=|Rk|,则算法收敛。


    注意事项:

    1.上面也讲过,字典的原子是归一化的,也就是||di||=1,因为我们选取max<R,di>时,如果di长度不统一,不能得出最好的投影。

    2.如果我们的字典只有两个向量d1,d2,那么MP算法会在这两个向量间交叉迭代投影,也就是f=a1d1+a2d2+a3d1+a4d2+…..;也就是之前投影过的原子方向,之后还有可能投影。换句话说,MP的方向选择不是最优的,是次优的。

    如下图:


    这也是其改进版本OMP要改进的地方。

     

    OMP算法

    也就是正交的MP算法。

    MP算法的次最优性来源其残差只与当前投影方向垂直,这样在接下来的投影中,很有可能会再次投影到原来的方向。

    于是,在投影时,如果我们使得残差Rk+1与x1-xk+1的所有向量垂直,则可以克服这个问题,如下:

     


    求解方法

    假设我们已经得到了第k步的最优解:


    我们要继续更新到第k+1步,目标是得到:


    需要注意的是,我们下一步更新时,之前原子的系数 也要更新,否则不能满足约束。

    于是我们需要求得如何更新之前原子系数 ,以及如何求得下一个投影方向 。

     





    收敛性:

    同样根据勾股定理,得到如下:

    于是算法收敛。

     

    具体步骤:


     

    最后,贴一个sparse求解的工具包,里面包含了MP,OMP算法的代码:

    http://spams-devel.gforge.inria.fr/

     

    参考文献:

    http://lear.inrialpes.fr/people/mairal/tutorial_iccv09/

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

    OrthogonalMatching Pursuit:Recursive Function Approximation with Applications to WaveletDecomposition

    展开全文
  • OMP算法MATLAB

    2017-12-03 20:27:03
    结合http://blog.csdn.net/jbb0523/article/details/45130793,由此转载
  • OMP算法笔记

    千次阅读 2019-06-20 14:56:47
    OMP算法笔记欢迎使用Markdown编辑器新的改变功能快捷键合理的创建标题,有助于目录的生成如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容居中、居左、居右...
  • 压缩感知的稀疏重构中广泛应用的正交匹配追踪(OMP算法matlab程序,该算法由香港大学电子工程系 沙威老师开发,代码注释详细,便于读者理解。已测试,可以正常运行。读者通过代码可以加深对该算法以及压缩感知、...
  • OMP算法MATLAB程序

    热门讨论 2012-02-20 18:49:21
    很好用的OMP算法程序,经过本人多次使用 ,很实用
  • OMP算法原理

    2012-03-05 10:37:33
    OMP算法原理,基于稀疏二进制序列的低密度奇偶校验码
  • OMP算法: % OMP的函数 % s-测量;T-观测矩阵;N-向量大小 function hat_y=omp_fun(s,T,K) N = size(T,2); Size=size(T); % 观测矩阵大小 M=Size(1); % 测量 ...
  • MP算法和OMP算法及其思想与实现

    万次阅读 2017-01-09 16:29:44
    主要介绍MP(Matching Pursuits)算法OMP(Orthogonal Matching Pursuit)算法[1],这两个算法虽然在90年代初就提出来了,但作为经典的算法,国内文献(可能有我没有搜索到)都仅描述了算法步骤和简单的应用,并未对其...
  • 压缩感知重构算法之OMP算法python实现

    万次阅读 热门讨论 2016-03-18 15:08:55
    压缩感知重构算法之OMP算法python实现 压缩感知重构算法之CoSaMP算法python实现 压缩感知重构算法之SP算法python实现 压缩感知重构算法之IHT算法python实现 压缩感知重构算法之OLS算法python实现 压缩感知重构...
  • 稀疏分解中的MP与OMP算法

    万次阅读 多人点赞 2017-04-12 15:11:57
    主要介绍MP与OMP算法的思想与流程,解释为什么需要引入正交? !!今天发现一个重大问题,是在读了博主的正交匹配追踪(OMP)在稀疏分解与压缩感知重构中的异同,...
  • C语言实现OMP算法

    2017-12-28 17:03:34
    C语言实现的OMP(正交匹配追踪算法),1024长度数据恢复在1秒以内,重构效果很好
  • 二维图像OMP算法代码

    2015-10-24 09:38:35
    使用matlab代码实现的二维图像压缩感知OMP算法实现过程
  • OMP算法的C语言代码

    2013-11-04 20:46:43
    该资源详细的写了OMP的实现过程,程序采用C语言写成,有助于理解OMP的整个过程,并且有详细的代码注释!
  • OMP算法学习笔记

    千次阅读 2016-12-06 16:01:39
    OMP算法
  • OMP算法代码学习

    2017-10-23 14:03:00
    正交匹配追踪(OMP)算法的MATLAB函数代码并给出单次测试例程代码 测量数M与重构成功概率关系曲线绘制例程代码 信号稀疏度K与重构成功概率关系曲线绘制例程代码 参考来源:...
  • 压缩感知 (Compressed Sensing) 是一种利用信号普遍存在低维结构的先验知识,只需极少的采样点,就能以极...本实验主要涉及两部分代码,一部分为压缩感知的信号采样与重建(见test.m),另一部分为OMP算法(见OMP.m)。
  • 压缩感知的稀疏重构中广泛应用的正交匹配追踪(OMP算法matlab程序,该算法由香港大学电子工程系 沙威老师开发,代码注释详细,便于读者理解。已测试,可以正常运行。读者通过代码可以加深对该算法以及压缩感知、...
  • 压缩感知之OMP算法及DFT

    千次阅读 2015-05-09 23:20:50
    (转载)题目:Rachel_Zhang的“压缩感知”之 =================引言==================== 这段代码是香港大学沙威的代码,可以去他的空间看看:http://www.eee.hku.hk/~wsha/Freecode/freecode.htm,里面还有...
  • 压缩感知中OMP算法的C/C++实现

    千次阅读 2018-03-30 22:17:14
    压缩感知中OMP算法的C/C++实现 背景介绍 算法实现部分 总结 阅读之前注意: 本文阅读建议用时:30min 本文阅读结构如下表: 项目 下属项目 测试用例数量 背景介绍 ...
  • OMP_FACE——人脸识别系统 该软件包实现了基于稀疏表示的面部识别方法 程序相对便捷且易上手 主脚本中包含具体的一个例子 通常,通常遵循以下使用顺序即可实现人脸识别功能: 选择训练数据的数据库途径。 选择测试...
  • OMP算法的matlab实现

    千次阅读 2017-02-12 23:23:20
    function [x,res] = omp(b,A,k) % compute k-sparse approximation to b with matrix A using Matching pursuit [m,N] = size(A); x = zeros(N,1); res = b; support = []; A_get=zeros(m,k); for i=1:kcorr = (A')...
  • omp算法(matlab)稀疏表示中用来求最优解。比较好,有demo
  • 正交匹配追踪 OMP 算法原理分析

    万次阅读 2013-11-07 14:00:07
    % 1-D信号压缩传感的实现(正交匹配追踪法Orthogonal Matching Pursuit) % 测量数M>=K*log(N/K),K是稀疏度,N信号长度,可以近乎完全重构 % 编程人--香港大学电子工程系 沙威 Email: wsha@eee.hku.hk ...

空空如也

1 2 3 4 5 ... 20
收藏数 7,815
精华内容 3,126
关键字:

omp算法