Add "Macroscopic Analysis of VAMP in a MMS"
This commit is contained in:
Ksyer
@@ -1,77 +1,77 @@
|
|||||||
% Input:y,A,lambda,tau,Kit
|
% Input:y,A,lambda,tau,Kit
|
||||||
% Output:x_hat_wl,x_hat_d
|
% Output:x_hat_wl,x_hat_d
|
||||||
|
|
||||||
function [x_hat_wl, x_hat_d] = cVAMP(y, A, lambda, tau, Kit)
|
function [x_hat_wl, x_hat_d] = cVAMP(y, A, lambda, tau, Kit)
|
||||||
% Initialization
|
% Initialization
|
||||||
gamma = 768 ./ 1024;
|
gamma = 768 ./ 1024;
|
||||||
k = 0;
|
k = 0;
|
||||||
p = ctranspose(A) * y;
|
p = ctranspose(A) * y;
|
||||||
h_1 = p;
|
h_1 = p;
|
||||||
Q_1 = gamma;
|
Q_1 = gamma;
|
||||||
tau_d = 1;
|
tau_d = 1;
|
||||||
|
|
||||||
% while
|
% while
|
||||||
while (k < Kit) && (tau_d > tau)
|
while (k < Kit) && (tau_d > tau)
|
||||||
% Factorized Part
|
% Factorized Part
|
||||||
x_1 = ST(h_1, lambda, Q_1); % \hat{x}_1^{(k)}
|
x_1 = ST(h_1, lambda, Q_1); % \hat{x}_1^{(k)}
|
||||||
chi_1 = F1(x_1, lambda, Q_1); % \chi_1^{(k)}
|
chi_1 = F1(x_1, lambda, Q_1); % \chi_1^{(k)}
|
||||||
% Message Passing
|
% Message Passing
|
||||||
h_2 = x_1 ./ chi_1 - h_1; % h_2^{(k)}
|
h_2 = x_1 ./ chi_1 - h_1; % h_2^{(k)}
|
||||||
Q_2 = 1 ./ chi_1 - Q_1; % \hat{Q}_2^{(k)}
|
Q_2 = 1 ./ chi_1 - Q_1; % \hat{Q}_2^{(k)}
|
||||||
% Gaussian Part
|
% Gaussian Part
|
||||||
t1 = (p + h_2) ./ Q_2;
|
t1 = (p + h_2) ./ Q_2;
|
||||||
t2 = ctranspose(A) * (A * (p + h_2)) / ((Q_2 + 1) * Q_2);
|
t2 = ctranspose(A) * (A * (p + h_2)) / ((Q_2 + 1) * Q_2);
|
||||||
x_2 = t1 + t2; % \hat{x}_2^{(k)}
|
x_2 = t1 + t2; % \hat{x}_2^{(k)}
|
||||||
chi_2 = gamma ./ (Q_2 + 1) + (1 - gamma) ./ Q_2;
|
chi_2 = gamma ./ (Q_2 + 1) + (1 - gamma) ./ Q_2;
|
||||||
% Message Passing
|
% Message Passing
|
||||||
h_1_next = x_2 ./ chi_2 - h_2;
|
h_1_next = x_2 ./ chi_2 - h_2;
|
||||||
Q_1_next = 1 ./ chi_2 - Q_2;
|
Q_1_next = 1 ./ chi_2 - Q_2;
|
||||||
tau_d = norm(h_1_next - h_1) ./ norm(h_1_next);
|
tau_d = norm(h_1_next - h_1) ./ norm(h_1_next);
|
||||||
k = k + 1;
|
k = k + 1;
|
||||||
% output
|
% output
|
||||||
x_hat_wl = x_1;
|
x_hat_wl = x_1;
|
||||||
x_hat_d = h_1_next ./ Q_1_next;
|
x_hat_d = h_1_next ./ Q_1_next;
|
||||||
% next
|
% next
|
||||||
h_1 = h_1_next;
|
h_1 = h_1_next;
|
||||||
Q_1 = Q_1_next;
|
Q_1 = Q_1_next;
|
||||||
end
|
end
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
% SoftThreshold function
|
% SoftThreshold function
|
||||||
function x = ST(h_1, lambda, Q_1)
|
function x = ST(h_1, lambda, Q_1)
|
||||||
[N, M] = size(h_1);
|
[N, M] = size(h_1);
|
||||||
x = zeros(N, M);
|
x = zeros(N, M);
|
||||||
|
|
||||||
for i = 1:N
|
for i = 1:N
|
||||||
% sign = h_1(i) ./ abs(h_1(i));
|
% sign = h_1(i) ./ abs(h_1(i));
|
||||||
diff = abs(h_1(i)) - lambda(i);
|
diff = abs(h_1(i)) - lambda(i);
|
||||||
x(i) = sign(h_1(i)) .* (diff ./ Q_1) .* SF(diff);
|
x(i) = sign(h_1(i)) .* (diff ./ Q_1) .* SF(diff);
|
||||||
end
|
end
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
% Heaviside's step function
|
% Heaviside's step function
|
||||||
function v = SF(a)
|
function v = SF(a)
|
||||||
% if a > 0
|
% if a > 0
|
||||||
% v = 1;
|
% v = 1;
|
||||||
% elseif a == 0
|
% elseif a == 0
|
||||||
% v = 0; % at zero points
|
% v = 0; % at zero points
|
||||||
% else
|
% else
|
||||||
% v = 0;
|
% v = 0;
|
||||||
% end
|
% end
|
||||||
v = heaviside(a);
|
v = heaviside(a);
|
||||||
end
|
end
|
||||||
|
|
||||||
% SoftThreshold function
|
% SoftThreshold function
|
||||||
function v = F1(x_1, lambda, Q_1)
|
function v = F1(x_1, lambda, Q_1)
|
||||||
[N, M] = size(x_1);
|
[N, M] = size(x_1);
|
||||||
count = 0;
|
count = 0;
|
||||||
|
|
||||||
for i = 1:N
|
for i = 1:N
|
||||||
temp = Q_1 .* abs(x_1(i)) + lambda(i);
|
temp = Q_1 .* abs(x_1(i)) + lambda(i);
|
||||||
count = count + (2 - lambda(i) ./ temp) .* SF(abs(x_1(i)));
|
count = count + (2 - lambda(i) ./ temp) .* SF(abs(x_1(i)));
|
||||||
end
|
end
|
||||||
|
|
||||||
v = count ./ (2 .* N .* Q_1);
|
v = count ./ (2 .* N .* Q_1);
|
||||||
end
|
end
|
||||||
|
|||||||
@@ -1,46 +1,46 @@
|
|||||||
% Input:y,A,lambda,tau,Kit
|
% Input:y,A,lambda,tau,Kit
|
||||||
% Output:x_hat_wl,x_hat_d
|
% Output:x_hat_wl,x_hat_d
|
||||||
clear;
|
clear;
|
||||||
clc;
|
clc;
|
||||||
% rng(1); % 随机种子
|
% rng(1); % 随机种子
|
||||||
|
|
||||||
%% test_稀疏向量
|
%% test_稀疏向量
|
||||||
% 设定稀疏度
|
% 设定稀疏度
|
||||||
k = 100; % 设定稀疏度
|
k = 100; % 设定稀疏度
|
||||||
|
|
||||||
% 构造感知矩阵D
|
% 构造感知矩阵D
|
||||||
m = 768; % 感知矩阵行数
|
m = 768; % 感知矩阵行数
|
||||||
n = 1024; % 感知矩阵列数 （n>>m）
|
n = 1024; % 感知矩阵列数 （n>>m）
|
||||||
|
|
||||||
% D = randn(m,n); % 生成满足高斯分布的感知矩阵 64*256
|
% D = randn(m,n); % 生成满足高斯分布的感知矩阵 64*256
|
||||||
|
|
||||||
F = dftmtx(n);
|
F = dftmtx(n);
|
||||||
row_indices = randperm(n, m);
|
row_indices = randperm(n, m);
|
||||||
D = F(row_indices, :);
|
D = F(row_indices, :);
|
||||||
|
|
||||||
% 构造稀疏信号X——共n个元素，其中k个元素不为0
|
% 构造稀疏信号X——共n个元素，其中k个元素不为0
|
||||||
X = zeros(n, 1);
|
X = zeros(n, 1);
|
||||||
index = randperm(n, k);
|
index = randperm(n, k);
|
||||||
val = randn(1, k);
|
val = randn(1, k);
|
||||||
X(index) = val;
|
X(index) = val;
|
||||||
|
|
||||||
% 得到观测矩阵（压缩后）
|
% 得到观测矩阵（压缩后）
|
||||||
A = D * X;
|
A = D * X;
|
||||||
|
|
||||||
%%
|
%%
|
||||||
% % 通过cVMAP算法完成恢复X，得到恢复后信号
|
% % 通过cVMAP算法完成恢复X，得到恢复后信号
|
||||||
lambda = ones(n, 1) ./ 10;
|
lambda = ones(n, 1) ./ 10;
|
||||||
[x_hat_wl, x_hat_d] = cVAMP(A, D, lambda, 1e-4, 200);
|
[x_hat_wl, x_hat_d] = cVAMP(A, D, lambda, 1e-4, 200);
|
||||||
|
|
||||||
%%
|
%%
|
||||||
% 显示结果
|
% 显示结果
|
||||||
figure;
|
figure;
|
||||||
subplot(3, 1, 1)
|
subplot(3, 1, 1)
|
||||||
stem(X);
|
stem(X);
|
||||||
title('origin signal')
|
title('origin signal')
|
||||||
subplot(3, 1, 2)
|
subplot(3, 1, 2)
|
||||||
stem(x_hat_wl);
|
stem(x_hat_wl);
|
||||||
title('restored signal')
|
title('restored signal')
|
||||||
subplot(3, 1, 3)
|
subplot(3, 1, 3)
|
||||||
stem(X - x_hat_wl);
|
stem(X - x_hat_wl);
|
||||||
title('differ')
|
title('differ')
|
||||||
|
|||||||
Reference in New Issue
Block a user