压缩感知算法预测

❤️学妹发来求助消息，说压缩感知算法程序调不出来。一开始以为是数据没有归一化，后面发现这并不是导致算法效果不好的原因。
通过一波疯狂调试后发现，出问题的主要原因在于自己输入的仿真样本与示例中给的样本特征不一致，需要调整仿真样本特征，并均衡样本（这点很重要）。看到预想结果后学妹开心及了！

如果觉得博主分享的不错，希望能留下您的一键❤️三连❤️（点赞+评论+收藏） ，您的支持就是我前进的动力❤️，您的三连对我特别重要。

一、背景前言：

基于相关向量机的叶贝斯压缩感知主要步骤如下：
1、首先，待重构信号 x 为可压缩信号，即在稀疏基表示下为 N 维稀疏列向量 s，此稀疏向量 s 中为数极少的非零元素可表征此信号的主要特征。由于噪声使得稀疏信号 s 中掺杂部分微小幅值。故：
     y=Φs=Φ(s_s + s_e)=Φs_s + n_e      (1)
其中，n_e中的元素近似服从均值为零的高斯分布。
2、在压缩感知中的观测噪声n_m满足均值为零的高斯分布，故式（1）可转换为
    y=Φs_s + n_e + n_m = Φs_s + n      (2)
其中，n近似为一个均值为零，方差为σ²的高斯噪声。
3、y的高斯似然模型为：
    p(y|s_s,σ²)=(2πσ²)^-k/2exp(-$\frac{1}{2{\sigma }^2}$||y-Φs||²)     (3)
因此，压缩感知中的稀疏权重 S_s转换问题就变成了带有S_s稀疏约束的线性回归问题。
假设Φ已知，那么基于y的压缩感知测量中需要估计的量是稀疏权重s_s以及噪声方差σ²。
4、相关向量机采用层次先验稀疏模型来促进稀疏性并满足共轭先验，对x中的每一个元素均假设为均值为0的高斯先验
    p(x|α)= ${\prod }_i^n(N(x_i|0,{\alpha }_i^{-1})$     (4)
    p(α|a,b)= ${\prod }_{i=1}^n(\Gamma ({\alpha }_i|a,b)$     (5)
5、对超参数α进行边缘分布计算，得到x先验
    p(x|a,b)= ${\prod }_{i=1}^N{\int }_0^{\infty }(N(x_i|0,{{\alpha }_i}^{-1})\ast \Gamma ({\alpha }_i|a,b)d{\alpha }_i$      (6)
6、类似的，为σ_n²选择 $\Gamma ({\alpha }_0|c,d)$先验
7、在给定观测信号y和测量矩阵Φ，根据贝叶斯公式， x的后验是多变量的高斯分布，其均值和方差为：
    $\mu ={\alpha }_0\sum {\Phi }^{T}y$     (7)
    $\sum =({\alpha }_0{\Phi }^T\Phi +A{)}^{-1}$     (8)
其中$A=diag({\alpha }_0,{\alpha }_1...{\alpha }_N)$
8、通过type-II ML逼近过程计算参数${\alpha }_i和{\alpha }_0$

二、 Matlab代码

1、Matlab数据生成函数（data.m）

% function [input, output ] = data( N )
%     input=(rand(N,1)-0.5)*4;     %生成【-2,2】的均匀分布随机数
%     output=sin(input*5);%+randn(N,1)/4;     %randn表示生成正态分布的随机数
% end

%===============================
function [input, output ] = data( N )
    f1=50; 
    f2=100; 
    f3=200; 
    fs = 3200;
    ts=1/fs;%采样频率为3200，1s内可以采集3200个点

    temp = randperm(fs);%产生1~3200的随机不重复整数（3200个）
    Tn = 200;%200s内
    input=(temp(1:N)'*ts*Tn);

    output=0.3*cos(2*pi*f1*input*ts)+0.6*cos(2*pi*f2*input*ts)+0.1*cos(2*pi*f3*input*ts);

2、Matlab高斯内核函数（GussianKernal.m）

function GK = GussianKernal( r1, r2 )

    ler1=size(r1,1);
    ler2=size(r2,1);
    GK=ones(ler1,ler2);

    for j=2:ler2
        for i = 1:ler1 
            GK(i,j)=exp(-(r1(i)-r2(j))^2);
        end
    end
 
end

3、主函数（SBLsinc.m）

clc;
clear;
%训练集样本数和测试集样本数的数据个数需要在一个采样周期内，即：Ntrain和Ntest ≤ fs=3200
Ntrain = 800;     %训练集样本数
Ntest = 2000;    %测试集样本数
[x,t] = data(Ntrain);  %并生成训练数据集
%t = mapminmax(t); %数据归一化
f =  GussianKernal(x,x);  %生成高斯核函数?
%  f =  GussianKernal(x',x');
m=size(f,2);
alpha=rand(1,m);
beta=rand();    %初始化学习参数  表示alpha0

for ep = 1:3000   %循环次数
    A = diag(alpha);
    F = pinv(beta*f'*f+A);   
    u = beta*F*f'*t;  %原程序
    % u = beta*F*f'*t';
    gamma = 1-alpha.*diag(F)';   %更新参数
    alpha_old=alpha;
    beta_old=beta;

    index=abs(alpha)<1e3;

    % a=gamma(index)   %测试学习
    alpha(index)=gamma(index)./(u(index)'.^2);
    beta=(Ntrain-sum(gamma))/((t-f*u)'*(t-f*u));%将数值小于1000的列标去除，只更新此部分
    % a=(Ntrain-sum(gamma));
    % b=((t-f*u)'*(t-f*u));
    % beta=(Ntrain-sum(gamma))/((t-f*u)*(t-f*u)');
    % beta(beta>10^10) = 10^10;
    
    err=max(abs(alpha(index)-alpha_old(index))./abs(alpha(index)))+abs(beta-beta_old)/abs(beta);
    fprintf('=======>第%d次迭代，误差为%d\n',ep,err);
        if err<0.000001
            %gamma = 1-alpha.*diag(F)';   %更新参数
            break;
        end %判断如果满足收敛条件，停止迭代
end

tpre=f(:,index)*u(index);  %计算训练数据集
%  tpre=u;

figure(1);
plot(x,t,'r+');
hold on;
plot(x,tpre,'b.');
hold on;
plot(x(index(2:end)),t(index(2:end)),'ko'); %分别绘制原始数据，预测数据和相关向量
% plot(x,t,'ko');

title('Testing Data1');
xlabel('x');
ylabel('t');
legend('训练数据','预测数据','相关向量');%标注
wc=abs(norm(tpre)-norm(t))/norm(t);

[xte,tte]=data(Ntest);%生成500个测试数据集
ftest =  GussianKernal(xte,x);%生成关联矩阵
ttepre=ftest(:,index)*u(index);%预测
WC=abs(norm(ttepre)-norm(tte))/norm(tte);
% ttepre=u;
figure(2)
plot(xte,tte,'r+');
%plot(xte,tte,'r+');
hold on;
plot(xte,ttepre,'b.');
hold on;
plot(x(index(2:end)),t(index(2:end)),'ko');%分别绘制原始数据，预测数据和相关向量（由200个数据绘出）
%  plot(x,t,'ko');

title('Testing Data2');
xlabel('x');
ylabel('t');
legend('训练数据','预测数据','相关向量');%标注

三、测试效果

1、示例结果

（1）迭代次数：

（2）初始数据打印

（3）预测结果打印

2、个人数据结果

（1）迭代次数：

（2）初始数据打印

（3）预测结果打印