/
Experiment_Effect_p_Gaussian.m
150 lines (114 loc) · 5.04 KB
/
Experiment_Effect_p_Gaussian.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
% Effect of the sparsity rate given sparse signal is Gaussian distributed
%
% This code demonstrates the effect of sparsity rate "p" on the
% performance of the algorithm described in "Structure-Based Bayesian
% Sparse Reconstruction, Ahmed A. Quadeer & Tareq Y. Al-Naffouri." and
% compares it with other sparse reconstruction algorithms for the case when
% the sparse signal is Gaussian distributed.
%
% Coded by: Ahmed Abdul Quadeer
% E-mail: ahmedaq@gmail.com
% Last change: Dec. 12, 2012
% Copyright (c) Ahmed Abdul Quadeer, Tareq Y. Al-Naffouri, 2012
clear all;close all;clc
N = 1000; % Length of sparse signal
M = N/4; % Number of Measurements
s = 1:2:20; % Sparsity level
F = 1/sqrt(N)*exp(-sqrt(-1)*(2*pi)*(0:1:N-1)'*(0:1:N-1)/N); %DFT Matrix
ITER = 5000; % Number of iterations
p = s./N; % Sparsity rate
N0 = 1; % Noice variance
%Initialization
nmse_l1_ref = zeros(1,length(s));
nmse_fbmp = zeros(1,length(s));
nmse_ocmg = zeros(1,length(s));
nmse_omp_ref = zeros(1,length(s));
t_l1_ref = zeros(1,length(s));
t_fbmp = zeros(1,length(s));
t_ocmg = zeros(1,length(s));
t_omp_ref = zeros(1,length(s));
for ss = 1:1:length(s)
%Initialization
error_l1_ref = zeros(1,ITER);
error_fbmp = zeros(1,ITER);
error_ocmg = zeros(1,ITER);
error_omp_ref = zeros(1,ITER);
time_l1_ref = zeros(1,ITER);
time_fbmp = zeros(1,ITER);
time_ocmg = zeros(1,ITER);
time_omp_ref = zeros(1,ITER);
for iter = 1:ITER
%% Constructing the sparse signal
x = zeros(N,1);
impulse_place = randi(N,1,s(ss)); % Takes s random values out of N
sigma_imp = 1000; % Sparse signal variance
x(impulse_place) = sqrt(sigma_imp/2)*(randn(size(impulse_place)) + 1i*randn(size(impulse_place)));
%% Measured Signal
AA = eye(N);
index= 1:M;
A = AA(index,:);
Psi = A*F; % Measurement matrix
n = sqrt(N0/2)*randn(M,1)+1i*randn(M,1); % Noise vecror
y = Psi*x + n; % Measurement vector
%% L1 norm minimization using CVX (Given b and A, find x_hat such that A*x_hat = b)
tic
cvx_begin
cvx_quiet(true) % For supressing results displayed in command window
variable x_hat(N) complex; % Defining the estimation variable with its length in brackets
epsilon = sqrt(N0*(M+2*sqrt(2*M)));
minimize(norm(x_hat,1)); % Problem Definition
subject to
norm(y-Psi*x_hat,2) <= epsilon;
cvx_end
J = BLC_Gaussian(x_hat,y,p(ss),Psi',N0,sigma_imp); % Support Refinement (x|Supp. Gaussian)
x_hat_ref = Refinement_MMSE(J,y,Psi',N0,sigma_imp); % Amplitude Refinement (MMSE)
time_l1_ref(iter) = toc;
error_l1_ref(iter) = norm(x-x_hat_ref)^2/norm(x)^2;
%% fbmp
sig2s = [0; sigma_imp]; % Sparse coefficient variances [off;on]
mus = [0; 0]; % sparse coefficient means [off;on;...;on]
D = 10; % Greedy branches
stop = 0; % Stop criteria
tic
x_fbmp = fbmpc_fxn_reduced(y, Psi, p(ss), N0, sig2s, mus, D, stop);
time_fbmp(iter) = toc;
error_fbmp(iter) = norm(x-x_fbmp)^2/norm(x)^2;
%% OC
tic
x_ocmg = OC_Gaussian_seq(y,Psi,sigma_imp,N0,p(ss),index);
time_ocmg(iter) = toc;
error_ocmg(iter) = norm(x-x_ocmg)^2/norm(x)^2;
%% OMP
tic
sup_omp = greed_omp(y,Psi,N);
J_omp = BLC_Gaussian(sup_omp,y,p(ss),Psi',N0,sigma_imp); % Support Refinement
x_omp_ref = Refinement_MMSE(J_omp,y,Psi',N0,sigma_imp); % Amplitude Refinement
time_omp_ref(iter) = toc;
error_omp_ref(iter) = norm(x-x_omp_ref)^2/norm(x)^2;
end
%Update
nmse_l1_ref(ss) = 10*log10(sum(error_l1_ref)/ITER)
nmse_fbmp(ss) = 10*log10(sum(error_fbmp)/ITER)
nmse_ocmg(ss) = 10*log10(sum(error_ocmg)/ITER)
nmse_omp_ref(ss) = 10*log10(sum(error_omp_ref)/ITER)
t_l1_ref(ss) = sum(time_l1_ref)/ITER
t_fbmp(ss) = sum(time_fbmp)/ITER
t_ocmg(ss) = sum(time_ocmg)/ITER
t_omp_ref(ss) = sum(time_omp_ref)/ITER
end
%% Plots
figure
plot(p,nmse_l1_ref,'g:d',p,nmse_fbmp,'r-.o',p,nmse_ocmg,'b-s',...
p,nmse_omp_ref,'k--p','LineWidth',2)
legend('CR','FBMP','OC','OMP')
axis([0 0.02 -25 0])
xlabel('\bf \it p')
ylabel('\bf NMSE (dB)')
grid on
figure
semilogy(p,t_l1_ref,'g:d',p,t_fbmp,'r-.o',p,t_ocmg,'b-s',...
p,t_omp_ref,'k--p','LineWidth',2)
legend('CR','FBMP','OC','OMP')
xlabel('\bf \it p')
ylabel('\bf Mean run-time')
grid on