omp可运行程序（m格式）

omp可运行程序(m格式)

function [y] = BCS（ xK ）
%BCS 贝叶斯压缩感知法（Bayesian Compressive Sensing）
%   x：原始信号
%   K：稀疏度，从FFT可以看出
%   y：重构信号
%  测量数M>=K*log（N/K）K是稀疏度N信号长度可以近乎完全重构
NN = length（x）;
% M = K*log（N/k）;
M = 64;

PHI = randn（MNN）;

fx = fft（x）;        %  求信号的频谱
fx_r = real（fx）;    %  频谱的实部，幅频信号
fx_i = imag（fx）;    %  频谱的虚部，相频信号    
t = PHI*fx_r‘;
sigma2 = std（t）^2/1e2;%标准差
eta = 1e-8;

% 寻找初始值
[NM] = size（PHI）;
PHIt = PHI‘*t;
PHI2 = sum（PHI.^2）‘;
ratio = （PHIt.^2）./PHI2;
[maxrindex] = max（ratio）;
alpha = PHI2（index）/（maxr-sigma2）;
% 计算初始均值、方差
phi = PHI（:index）;
Hessian = alpha + phi‘*phi/sigma2;
Sig = 1/Hessian;
mu = Sig*PHIt（index）/sigma2;
left = PHI‘*phi/sigma2;
S = PHI2/sigma2-Sig*left.^2;
Q = PHIt/sigma2-Sig*PHIt（index）/sigma2*left;
%
for count = 1:10000

    s = S; q = Q;
    s（index） = alpha.*S（index）./（alpha-S（index））;
    q（index） = alpha.*Q（index）./（alpha-S（index））;
    theta = q.^2-s;

    % 根据最大似然估计选择下一次迭代的alpha
    ml = -inf*ones（1M）;
    ig0 = find（theta>0）;
    % 重新估计的下标集
    [irefoowhich] = intersect（ig0index）;
    if ~isempty（ire）
        Alpha = s（ire）.^2./theta（ire）;
        delta = （alpha（which）-Alpha）./（Alpha.*alpha（which））;
        ml（ire） = Q（ire）.^2.*delta./（S（ire）.*delta+1）-log（1+S（ire）.*delta）;
    end
    % 保留的下标集
    iad = setdiff（ig0ire）;
    if ~isempty（iad）
        ml（iad） = （Q（iad）.^2-S（iad））./S（iad）+log（S（iad）./（Q（iad）.^2））;
    end
    is0 = setdiff（[1:M]ig0）;
    % 去除的下标集
    [idefoowhich] = intersect（is0index）;
    if ~isempty（ide）
        ml（ide） = Q（ide）.^2./（S（ide）-alpha（which））-log（1-S（ide）./alpha（which））;
    end

    [ML（count）idx] = max（ml）;
    % 检查是否满足停止迭代条件
    if count > 2 & abs（ML（count）-ML（count-1）） < abs（ML（count）-ML（1））*eta
        break;
    end

    % 更新 alphas
    which = find（index==idx）;
    if theta（idx） > 0
        if ~isempty（which） % 重新估计
            Alpha = s（idx）^2/theta（idx）;
            Sigii = Sig（whichwhich）; mui = mu（which）; Sigi = Sig（:which）;
            delta = Alpha-alpha（which）;
            ki = delta/（1+Sigii*delta）;
            mu = mu-ki*mui*Sigi;
            Sig = Sig-ki*Sigi*Sigi‘;
            comm = PHI‘*（phi*Sigi）/sigma2;
            S = S + ki*comm.^2;
            Q = Q + ki*mui*comm;
            %
            alpha（which） = Alpha;
        else % 保留
            Alpha = s（idx）^2/theta（idx）;
            phii = PHI（:idx）; Sigii = 1/（Alpha+S（idx））; mui = Sigii*Q（idx）;
            comm1 = Sig*（phi‘*phii）/sigma2;
            ei = phii-phi*comm1;
            off = -Sigii*comm1;
            Sig = [Sig+Sigii*comm1*comm1‘ off; off‘ Sigii];
            mu = [mu-mui*comm1; mui];
            comm2 = PHI‘*ei/sigma2;
            S = S - Sigii*comm2.^2;
            Q = Q - mui*comm2;
            index = [index;idx];
            alpha = [alpha;Alpha];
            phi = [phiphii];
        end
    else
        if ~isempty（which） % 去除
            Sigii = Sig（whichwhich）; mui = mu（which

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2017-06-17 21:27  可运行程序（m格式）\
     文件       19433  2017-06-17 21:32  可运行程序（m格式）\1.fig
     文件      186051  2017-06-17 21:32  可运行程序（m格式）\1.tif
     文件       20464  2017-06-17 21:32  可运行程序（m格式）\2.fig
     文件      178469  2017-06-17 21:32  可运行程序（m格式）\2.tif
     文件       28695  2017-06-17 21:32  可运行程序（m格式）\3.fig
     文件      186031  2017-06-17 21:32  可运行程序（m格式）\3.tif
     文件       20627  2017-06-17 21:32  可运行程序（m格式）\4.fig
     文件      189167  2017-06-17 21:32  可运行程序（m格式）\4.tif
     文件       28840  2017-06-17 21:32  可运行程序（m格式）\5.fig
     文件      186289  2017-06-17 21:32  可运行程序（m格式）\5.tif
     文件       16919  2017-06-17 21:32  可运行程序（m格式）\6.fig
     文件      108603  2017-06-17 21:32  可运行程序（m格式）\6.tif
     文件        3692  2017-01-15 11:03  可运行程序（m格式）\BCS.m
     文件        3527  2017-01-15 11:05  可运行程序（m格式）\BCS2.m
     文件        2154  2017-01-14 16:52  可运行程序（m格式）\OMP.m
     文件        2184  2017-01-14 21:43  可运行程序（m格式）\OMP2.m
     文件         281  2017-01-10 19:28  可运行程序（m格式）\lfm_sign.m
     文件         348  2017-01-11 00:04  可运行程序（m格式）\n_fre_sign.m
     文件        3149  2017-01-15 10:59  可运行程序（m格式）\rebuild.m
     文件        4784  2017-01-15 10:59  可运行程序（m格式）\rebuild_test.m
     文件        4773  2017-01-15 10:59  可运行程序（m格式）\rebuild_test2.m