###### ao~

2020-09-17 13:54 阅读 61

# 求大神帮忙看一下哪里出了问题，怎么运行不了？

%% 参数配置初始化
CNT = 10;%对于每组(K,M,N)，重复迭代次数
N = 256;%信号x的长度
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);%本次测量值个数
K1 = 1:5:ceil(M/2);%信号x的稀疏度K没必要全部遍历，每隔5测试一个就可以了
K2 = 1:5:ceil(M/2);
K_set = [K1,K2];
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次
%x1= gene_x(N/2,K1);
%x2= 100 * gene_x(N/2,K2);
x1= gene_x(N/2,K1(cnt));
x2= 100 * gene_x(N/2,K2(cnt));
x = [x1;x2];
Phi = randn(M,N/2);%测量矩阵为高斯矩阵
Phi = orth(Phi')'; %列归一化矩阵
A1 = Phi;
A2 = Phi;
A = [A1,A2]; %A分为两片
y = A * x;

``````             x_r= zeros(N/2,1);%存储分片恢复信号
x_r= {x_r;x_r};
``````

%parfor i=1:2 用法不对，应该对原来的函数用parfor
for i=1:2
x_r{i}=nnomp(x(1+(i-1)*N/2:i*N/2),K(i),M,A(:,1+(i-1)*N/2:i*N/2));
end
toc
x_r=[x_r{1};x_r{2}];
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)');

function x = gene_x(N,K)
% 稀疏度K，x长度N
%生成随机稀疏x

%N=256;
%k=10
x=zeros(N,1);
Index_K = randperm(N);
% x(Index_K(1:K)) = 5 * abs(randn(K,1)); %x为K稀疏的,且位置是随机的,均值为0,方差为5的高斯噪声
x(Index_K(1:K)) = rand(K,1);
% x(Index_K(K/2+1:K)) =100 * rand(K/2,1);
end

function theta=nnomp(x,K,M,Phi,epsilon,loopmax)
%theta恢复，pos_array_位置，x信号，K稀疏度，M测量数
if nargin < 6
epsilon = 1e-5;
end
if nargin < 5
loopmax = 10000;
end
loop = 0;
N=length(x);
Phi =randn(M,N); %观测矩阵为均值为0，方差为1的高斯矩阵
Phi = orth(Phi')'; %列归一化矩阵
% e = 1e-3;

y= Phi*x; %获得线性测量

%% nnomp重构信号

``````[y_rows,y_columns] = size(y);
if y_rows<y_columns
y = y';%y should be a column vector
end
theta = zeros(N,1);%用来存储恢复的theta(列向量)
Pos_theta = [];%用来迭代过程中存储A被选择的列序号
r_n = y;%初始化残差为y
times = 1;
``````

while 1
loop = loop + 1;
% for kk=1:K%最多迭代K次
%(1) Identification
for col = 1:N
product(col) = Phi(:,col)'*r_n;
end
% product = Phi'*r_n;%传感矩阵A各列与残差的内积
a = max(product,0);
[val,pos]=sort(a,'descend');
Js = pos(1:K);%选出内积值最大的K列
%(2) Support Merger
Is = union(Pos_theta,Js);%Pos_theta与Js并集
%(3) Estimation
%At的行数要大于列数，此为最小二乘的基础(列线性无关)
if length(Is)<=M
At = Phi(:,Is);%将A的这几列组成矩阵At
else%At的列数大于行数，列必为线性相关的,At'*At将不可逆
if kk == 1
theta_ls = 0;
end
break;%跳出for循环
end
%y=At*theta，以下求theta的非负最小二乘解(Least Square)
theta_ls = lsqnonneg(At,y);
%(4) Pruning
[val,pos] = sort(theta_ls,'descend');
%(5) Sample Update
Pos_theta = Is(pos(1:K));
theta_ls = theta_ls(pos(1:K));
%At(:,pos(1:K))*theta_ls是y在At(:,pos(1:K))列空间上的正交投影
r_n = y - At(:,pos(1:K))*theta_ls;%更新残差

``````    times = times +1;
if norm(r_n) < epsilon
break
end
if loop >= loopmax
fprintf('loop > %d\n',loopmax);
break
end
``````

% if norm(r_n)<1e-5%Repeat the steps until r=0
% break;%跳出for循环
% end
end
theta(Pos_theta)=theta_ls;%恢复出的theta
end

% xr_out = theta;

• 点赞
• 写回答
• 关注问题
• 收藏
• 复制链接分享