Compare commits
2
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Ksyer
|
17386d80c1 | ||
|
Ksyer
|
b91a609ca5 |
@@ -33,7 +33,6 @@ Manifest.toml
|
||||
*.asv
|
||||
*.pdf
|
||||
|
||||
0 - example
|
||||
*.mat
|
||||
*.fig
|
||||
*.tif
|
||||
|
||||
Submodule 0 - example/CS-Recovery-Algorithms deleted from 24ed57443d
Submodule 0 - example/CompressedSensing.jl deleted from addf3f147a
@@ -0,0 +1,82 @@
|
||||
% Input: y,A,lambda,tau,Kit
|
||||
% Output: x_hat_wl,x_hat_d
|
||||
|
||||
% Main structure of cVAMP
|
||||
function [x_hat_wl, x_hat_d] = cVAMPro(y, A, lambda, tau, Kit)
|
||||
|
||||
% Initialization
|
||||
[M, N] = size(A);
|
||||
gamma = M / N;
|
||||
k = 0;
|
||||
p = ctranspose(A) * y;
|
||||
h_1 = p;
|
||||
Q_1 = gamma;
|
||||
tau_d = 1;
|
||||
|
||||
% Iteration
|
||||
while ((k < Kit) && (tau_d > tau))
|
||||
% Factorized Part
|
||||
x_1 = ST(h_1, lambda, Q_1);
|
||||
chi_1 = F1(x_1, lambda, Q_1);
|
||||
% Message Passing
|
||||
h_2 = x_1 / chi_1 - h_1;
|
||||
Q_2 = 1 / chi_1 - Q_1;
|
||||
% Gaussian Part
|
||||
t1 = (p + h_2) / Q_2;
|
||||
t2 = ctranspose(A) * (A * (p + h_2)) / ((Q_2 + 1) * Q_2);
|
||||
x_2 = t1 - t2;
|
||||
chi_2 = gamma / (Q_2 + 1) + (1 - gamma) / Q_2;
|
||||
% Message Passing
|
||||
h_1_next = x_2 ./ chi_2 - h_2;
|
||||
Q_1_next = 1 / chi_2 - Q_2;
|
||||
tau_d = norm(h_1_next - h_1, Inf) / norm(h_1_next, Inf);
|
||||
k = k + 1;
|
||||
% output
|
||||
x_hat_wl = x_1;
|
||||
x_hat_d = h_1_next / Q_1_next;
|
||||
% next
|
||||
h_1 = h_1_next;
|
||||
Q_1 = Q_1_next;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
% SoftThreshold function
|
||||
function x = ST(h_1, lambda, Q_1)
|
||||
[N, M] = size(h_1);
|
||||
x = zeros(N, M);
|
||||
|
||||
for i = 1:N
|
||||
sign = h_1(i) ./ abs(h_1(i));
|
||||
diff = abs(h_1(i)) - lambda(i);
|
||||
x(i) = sign .* (diff ./ Q_1) .* SF(diff);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
% Heaviside's step function
|
||||
function v = SF(a)
|
||||
|
||||
if a > 0
|
||||
v = 1;
|
||||
elseif a == 0
|
||||
v = 0; % at zero points
|
||||
else
|
||||
v = 0;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
% Calculation of chi_1
|
||||
function chi_1 = F1(x_1, lambda, Q_1)
|
||||
[N, M] = size(x_1);
|
||||
count = 0;
|
||||
|
||||
for i = 1:N
|
||||
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) * (abs(x_1(i)) > 1e-4);
|
||||
end
|
||||
|
||||
chi_1 = count / (2 * N * Q_1);
|
||||
end
|
||||
@@ -0,0 +1,122 @@
|
||||
clear;
|
||||
close all;
|
||||
clc;
|
||||
|
||||
load test_Pd_Pfa_vamp_cal_stat.mat;
|
||||
|
||||
Fontsize = 18;
|
||||
plot_width = 800;
|
||||
plot_height = 600;
|
||||
Linewidth = 2;
|
||||
Markersize = 8;
|
||||
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
loglog(P_fa, P_fa_CROD(1,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
loglog(P_fa, P_fa_CROD(2,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_CROD(3,:), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_CROD(4,:), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_CROD(5,:), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('SNR = 0', 'SNR = 2', 'SNR = 4', 'SNR = 6', 'SNR = 8');
|
||||
xlabel('P_{fa} Set');
|
||||
ylabel('Actual P_{fa}');
|
||||
ylim([8e-5,1]);
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('cVAMPro P_{fa}')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(3);
|
||||
semilogx(P_fa_CROD(1,:), P_d_CROD(1,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
semilogx(P_fa_CROD(2,:), P_d_CROD(2,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_CROD(3,:), P_d_CROD(3,:), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_(3,:), P_d_CROD(4,:), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_CROD(3,:), P_d_CROD(5,:), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('SNR = 0', 'SNR = 2', 'SNR = 4', 'SNR = 6', 'SNR = 8');
|
||||
xlabel('P_{fa}');
|
||||
ylabel('P_{d}');
|
||||
xlim([8e-5,1]);
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('cVAMPro ROC')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(5);
|
||||
semilogy(SNR, P_fa_CROD(:,1), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
semilogy(SNR, P_fa_CROD(:,3), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_CROD(:,5), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_CROD(:,7), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_CROD(:,9), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('Pfa = 1e-4', 'Pfa = 1e-3', 'Pfa = 1e-2', 'Pfa = 1e-1', 'Pfa = 1');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{fa}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('cVAMPro')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(6);
|
||||
plot(SNR, P_d_CROD(:,1), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(SNR, P_d_CROD(:,3), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_CROD(:,5), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_CROD(:,7), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_CROD(:,9), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('Pfa = 1e-4', 'Pfa = 1e-3', 'Pfa = 1e-2', 'Pfa = 1e-1', 'Pfa = 1');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{d}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('cVAMPro')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
@@ -0,0 +1,11 @@
|
||||
function [ stat ] = stat_window(signal, win_size)
|
||||
n = length(signal);
|
||||
stat = zeros(n, 1);
|
||||
for i = 1: n
|
||||
for j = 1: win_size
|
||||
if i + j - 1 <= n
|
||||
stat(i, 1) = stat(i, 1) + signal(i + j - 1, 1);
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,280 @@
|
||||
clc;
|
||||
clear;
|
||||
close all;
|
||||
|
||||
|
||||
%% parameter setting
|
||||
|
||||
% SNR = 5: 1: 6;
|
||||
% SNR = 0: 2: 10;
|
||||
% SNR = 0: 5: 25;
|
||||
SNR = [0,5,8,10,12,15,20,25];
|
||||
len_SNR = length(SNR);
|
||||
SNR_t2 = 25;
|
||||
|
||||
rep_time = 2000;
|
||||
|
||||
% P_fa = 1e-1;
|
||||
P_fa = [1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1];
|
||||
len_P_fa = length(P_fa);
|
||||
|
||||
target_scattering = [0.8,1,0.9];
|
||||
win_size = length(target_scattering);
|
||||
|
||||
alpha_prop = 0.8;
|
||||
% alpha_prop = [0.6, 0.4, 0.2];
|
||||
|
||||
|
||||
lambda = 0.1;
|
||||
% Sparsity
|
||||
sigma_0 = 0.1;
|
||||
% noise sigma
|
||||
sigma_n = sigma_0;
|
||||
%% vamp set
|
||||
delta_VAMP = 1e-6;
|
||||
iter_max = 500;
|
||||
lambda_val = 0.1;
|
||||
|
||||
gamma = 0.5;
|
||||
|
||||
%% 参数设置
|
||||
B=5e5; %信号带宽10MHz
|
||||
Tp=100e-6; %脉宽100us
|
||||
fs=2*B; %采样频率
|
||||
Ts= 1 / fs; %采样周期
|
||||
K = B / Tp; %线性调频率
|
||||
fc = 1e8; %载波频率
|
||||
|
||||
Tr = 1e-3;
|
||||
|
||||
t = 0:1/fs:Tr-1/fs;
|
||||
t2 = 0:1/fs/2:Tr-1/fs/2;
|
||||
|
||||
c = 3e8; % 光速
|
||||
|
||||
distance_max = (Tr-Tp) * c / 2;
|
||||
|
||||
n = Tr * fs;
|
||||
m = round(n * gamma);
|
||||
|
||||
|
||||
%% experiment
|
||||
P_fa_CROD_cnt = zeros(len_SNR, len_P_fa, rep_time);
|
||||
|
||||
P_d_CROD_cnt = zeros(len_SNR, len_P_fa, rep_time);
|
||||
|
||||
|
||||
% h_thd = -log(P_fa);
|
||||
h_thd = chi2inv(1 - P_fa, 6) / 2;
|
||||
|
||||
%%
|
||||
parfor rep = 1: rep_time
|
||||
% for rep = 1: rep_time
|
||||
|
||||
for cnt_SNR = 1: len_SNR
|
||||
|
||||
% 设置目标位置
|
||||
% target_index = randi([2,n - length(target_scattering)]);
|
||||
target_index = 300;
|
||||
target2_index = round(n * 0.4);
|
||||
% 根据 SNR 设置散射点强度
|
||||
alpha = alpha_prop * sqrt(10^(SNR(cnt_SNR)/10) * sigma_n^2);
|
||||
alpha2 = alpha_prop * sqrt(10^(SNR_t2/10) * sigma_n^2);
|
||||
% 设置 x
|
||||
x = zeros(n ,1);
|
||||
x(target_index:target_index+2,1) = target_scattering * alpha;
|
||||
x(target2_index:target2_index+2,1) = target_scattering * alpha2;
|
||||
% plot(x)
|
||||
|
||||
%% 生成矩阵 A
|
||||
A_idx = randperm(n);
|
||||
A_idx = A_idx(1: m);
|
||||
A_idx = sort(A_idx);
|
||||
A = dftmtx(n);
|
||||
A = A(A_idx, :);
|
||||
A = A / sqrt(n);
|
||||
% % 匹配滤波放大倍数
|
||||
% multiple = A(:,1)' * A(:,1);
|
||||
|
||||
|
||||
%% 生成 y
|
||||
noise = random('Normal', 0, sigma_n/sqrt(2), m, 1) + 1j * random('Normal', 0, sigma_n/sqrt(2), m, 1);
|
||||
y = A * x + noise;
|
||||
|
||||
|
||||
%% cVAMPro 求解
|
||||
lambda = zeros(n,1) + lambda_val;
|
||||
[x_LASSO, x_hat_d_ro] = cVAMPro(y, A, lambda, delta_VAMP, iter_max);
|
||||
|
||||
|
||||
%%
|
||||
% CROD
|
||||
rho_active = sum(abs(x_LASSO) > 1e-3)/n;
|
||||
Q_hat = (gamma - rho_active)/(1 - rho_active);
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat = (gamma-Rho)/(1-Rho);
|
||||
x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat;
|
||||
% plot(abs(x_d_CROD))
|
||||
|
||||
|
||||
%%
|
||||
RSS = sum(abs(y - A * x_LASSO).^2)/m;
|
||||
chi = Rho*(1 - Rho)/(gamma - Rho);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_CROD = sqrt(2*chi_hat) / Q_hat;
|
||||
stat_CROD = abs(x_d_CROD / sigma_CROD).^2;
|
||||
|
||||
% window
|
||||
stat_extend = stat_window(stat_CROD, win_size);
|
||||
|
||||
figure(100)
|
||||
plot(stat_CROD);hold on;
|
||||
plot(stat_extend);hold off;
|
||||
|
||||
n_index = ones(size(stat_extend));
|
||||
t_index = zeros(size(stat_extend));
|
||||
|
||||
% distance_node = [target_index, target_index + 1, target_index + 2];
|
||||
distance_node = target_index;
|
||||
distance_area = target_index-2: 1: target_index+2;
|
||||
distance_area2 = target2_index-2: 1: target2_index+2;
|
||||
|
||||
n_index(distance_area) = 0;
|
||||
n_index(distance_area2) = 0;
|
||||
t_index(distance_node) = 1;
|
||||
|
||||
for cnt_h_th = 1: len_P_fa
|
||||
|
||||
figure(12);
|
||||
plot(stat_extend)
|
||||
hold on;
|
||||
kdline = zeros(size(stat_extend)) + h_thd(cnt_h_th);
|
||||
plot(kdline)
|
||||
hold off;
|
||||
|
||||
P_fa_CROD_cnt(cnt_SNR, cnt_h_th, rep) = sum(stat_extend(n_index>0) > h_thd(cnt_h_th)) / sum(n_index);
|
||||
P_d_CROD_cnt(cnt_SNR, cnt_h_th, rep) = sum(stat_extend(t_index>0) > h_thd(cnt_h_th)) / sum(t_index);
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
fprintf('%d\n', rep);
|
||||
|
||||
end
|
||||
|
||||
P_fa_CROD = mean(P_fa_CROD_cnt, 3);
|
||||
|
||||
P_d_CROD = mean(P_d_CROD_cnt, 3);
|
||||
|
||||
display(['设定虚警率',num2str(P_fa(1))])
|
||||
display(['实验虚警率',num2str(P_fa_CROD(1,1))])
|
||||
display(['实验检测率',num2str(P_d_CROD(1,1))])
|
||||
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
loglog(P_fa,P_fa_CROD(1,:), 'linewidth', 2);
|
||||
hold on;
|
||||
grid on;
|
||||
loglog(P_fa,P_fa_CROD(2,:), 'linewidth', 2);
|
||||
loglog(P_fa,P_fa_CROD(3,:), 'linewidth', 2);
|
||||
loglog(P_fa,P_fa_CROD(4,:), 'linewidth', 2);
|
||||
loglog(P_fa,P_fa_CROD(5,:), 'linewidth', 2);
|
||||
loglog(P_fa,P_fa_CROD(6,:), 'linewidth', 2);
|
||||
legend('SNR = 0','SNR = 2','SNR = 4','SNR = 6','SNR = 8','SNR = 10');
|
||||
xlabel('P_fa');
|
||||
ylabel('P_fa_CROD');
|
||||
|
||||
% figure(2);
|
||||
% loglog(P_fa_CROD_thry(1,:), P_d_CROD_thry(1,:), 'linewidth', 2);
|
||||
% hold on;
|
||||
% grid on;
|
||||
% loglog(P_fa_CROD_thry(2,:), P_d_CROD_thry(2,:), 'linewidth', 2);
|
||||
% loglog(P_fa_CROD_thry(3,:), P_d_CROD_thry(3,:), 'linewidth', 2);
|
||||
% legend('SNR = 0', 'SNR = 10', 'SNR = 20');
|
||||
% xlabel('P_{fa}');
|
||||
% ylabel('P_{d}');
|
||||
%
|
||||
% % figure(4);
|
||||
% % loglog(P_fa_CROD_thry(1,:), P_d_CROD_thry2(1,:), 'linewidth', 2);
|
||||
% % hold on;
|
||||
% % grid on;
|
||||
% % loglog(P_fa_CROD_thry(2,:), P_d_CROD_thry2(2,:), 'linewidth', 2);
|
||||
% % loglog(P_fa_CROD_thry(3,:), P_d_CROD_thry2(3,:), 'linewidth', 2);
|
||||
% % legend('SNR = 0', 'SNR = 10', 'SNR = 20');
|
||||
% % xlabel('P_{fa}');
|
||||
% % ylabel('P_{d}');
|
||||
%
|
||||
%
|
||||
% figure(3);
|
||||
% loglog(P_fa_CROD(1,:), P_d_CROD(1,:), 'linewidth', 2);
|
||||
% hold on;
|
||||
% grid on;
|
||||
% loglog(P_fa_CROD(2,:), P_d_CROD(2,:), 'linewidth', 2);
|
||||
% loglog(P_fa_CROD(3,:), P_d_CROD(3,:), 'linewidth', 2);
|
||||
% loglog(P_fa_CROD_thry(1,:), P_d_CROD_thry(1,:), 'linewidth', 2);
|
||||
% loglog(P_fa_CROD_thry(2,:), P_d_CROD_thry(2,:), 'linewidth', 2);
|
||||
% loglog(P_fa_CROD_thry(3,:), P_d_CROD_thry(3,:), 'linewidth', 2);
|
||||
% legend('rSNR = 0', 'rSNR = 10', 'rSNR = 20', ...
|
||||
% 'tSNR = 0', 'tSNR = 10', 'tSNR = 20');
|
||||
% xlabel('P_{fa}');
|
||||
% ylabel('P_{d}');
|
||||
|
||||
|
||||
save test_Pd_Pfa_vamp_cal_stat.mat ...
|
||||
SNR...
|
||||
P_fa...
|
||||
P_fa_CROD...
|
||||
P_d_CROD...
|
||||
lambda...
|
||||
m...
|
||||
n...
|
||||
h_thd...
|
||||
sigma_n...
|
||||
target_scattering;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,27 @@
|
||||
function [z] = FISTA(y, A, lambda, delta)
|
||||
|
||||
x_pre = A'*y;
|
||||
t = 1;
|
||||
z = x_pre;
|
||||
z_pre = z;
|
||||
t_pre = t;
|
||||
N = size(A, 2);
|
||||
diff = 1;
|
||||
E = eig(A'*A);
|
||||
L = E(end);
|
||||
temp1 = A'*y/L;
|
||||
temp2 = eye(N) - A'*A/L;
|
||||
k = 0;
|
||||
while((diff > delta) && (k < 1000))
|
||||
temp = temp1 + temp2 * z_pre;
|
||||
x = sft_thd(temp, lambda/L);
|
||||
t = 0.5*(1 + sqrt(1+4*t_pre*t_pre));
|
||||
z = x + (x - x_pre) * (t_pre-1) / t;
|
||||
diff = mean(abs(z_pre - z));
|
||||
x_pre = x;
|
||||
z_pre = z;
|
||||
t_pre = t;
|
||||
k = k + 1;
|
||||
end
|
||||
|
||||
end
|
||||
+66
@@ -0,0 +1,66 @@
|
||||
clear;
|
||||
close all;
|
||||
clc;
|
||||
|
||||
load test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.mat;
|
||||
|
||||
Fontsize = 18;
|
||||
plot_width = 800;
|
||||
plot_height = 600;
|
||||
Linewidth = 2;
|
||||
Markersize = 8;
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
plot(SNR, P_fa_CROD, '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(SNR, P_fa_CAMP, '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_fa_SDL, '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_fa_ROD, '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{fa}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman');
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
|
||||
figure(2);
|
||||
plot(SNR, P_d_CROD, '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(SNR, P_d_CAMP, '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_SDL, '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_ROD, '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{d}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman');
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
+247
@@ -0,0 +1,247 @@
|
||||
clc;
|
||||
clear;
|
||||
close all;
|
||||
|
||||
%% parameter setting
|
||||
m = 128;
|
||||
n = 256;
|
||||
|
||||
SNR = 0: 1: 15;
|
||||
len_SNR = length(SNR);
|
||||
|
||||
rep_time = 1e4;
|
||||
|
||||
P_fa = 1e-2;
|
||||
|
||||
p0 = 0.1;
|
||||
|
||||
lambda = 0.1;
|
||||
|
||||
sigma_0 = 0.1;
|
||||
|
||||
sigma_n = sigma_0;
|
||||
|
||||
%% experiment
|
||||
gamma = m/n;
|
||||
|
||||
P_fa_CROD_cnt = zeros(len_SNR, rep_time);
|
||||
P_fa_CAMP_cnt = zeros(len_SNR, rep_time);
|
||||
P_fa_SDL_cnt = zeros(len_SNR, rep_time);
|
||||
P_fa_ROD_cnt = zeros(len_SNR, rep_time);
|
||||
P_fa_LASSO_cnt = zeros(len_SNR, rep_time);
|
||||
|
||||
P_d_CROD_cnt = zeros(len_SNR, rep_time);
|
||||
P_d_CAMP_cnt = zeros(len_SNR, rep_time);
|
||||
P_d_SDL_cnt = zeros(len_SNR, rep_time);
|
||||
P_d_ROD_cnt = zeros(len_SNR, rep_time);
|
||||
P_d_LASSO_cnt = zeros(len_SNR, rep_time);
|
||||
|
||||
x_idx = rand(n, 1);
|
||||
if p0 == 0
|
||||
thd = -1;
|
||||
x_l0 = sum(x_idx > thd);
|
||||
else
|
||||
thd = sort(x_idx);
|
||||
thd = thd(round(n*p0));
|
||||
x_l1 = sum(x_idx <= thd);
|
||||
x_l0 = sum(x_idx > thd);
|
||||
end
|
||||
|
||||
x = zeros(n, 1);
|
||||
x_temp = sigma_0 * exp(1j*random('Uniform', 0, 2*pi, n, 1));
|
||||
x(x_idx <= thd) = x_temp(x_idx <= thd);
|
||||
x = x * sqrt(n/m);
|
||||
|
||||
h_thd = -log(P_fa);
|
||||
|
||||
parfor rep = 1: rep_time
|
||||
|
||||
A_idx = randperm(n);
|
||||
A_idx = A_idx(1: m);
|
||||
A_idx = sort(A_idx);
|
||||
A = dftmtx(n);
|
||||
A = A(A_idx, :);
|
||||
A = A / sqrt(n);
|
||||
|
||||
w = random('Normal', 0, sigma_0/sqrt(2), m, 1) + 1j * random('Normal', 0, sigma_0/sqrt(2), m, 1);
|
||||
% w = w / sqrt(10^(SNR(cnt_SNR)/10));
|
||||
|
||||
for cnt_SNR = 1: len_SNR
|
||||
|
||||
x1 = x * sqrt(10^(SNR(cnt_SNR)/10));
|
||||
|
||||
y = A * x1 + w;
|
||||
|
||||
% x_LASSO = LASSO_cvx(y, A, lambda);
|
||||
x_LASSO = FISTA(y, A, lambda, 1e-5);
|
||||
|
||||
% CROD
|
||||
rho_active = sum(abs(x_LASSO) > 1e-3)/n;
|
||||
Q_hat = (gamma - rho_active)/(1 - rho_active);
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat = (gamma-Rho)/(1-Rho);
|
||||
x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat;
|
||||
RSS = sum(abs(y - A * x_LASSO).^2)/m;
|
||||
chi = Rho*(1 - Rho)/(gamma - Rho);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_CROD = sqrt(2*chi_hat) / Q_hat;
|
||||
stat_CROD = abs(x_d_CROD / sigma_CROD).^2;
|
||||
|
||||
P_fa_CROD_cnt(cnt_SNR, rep) = sum(stat_CROD(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_CROD_cnt(cnt_SNR, rep) = sum(stat_CROD(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% CAMP
|
||||
Q_hat1 = gamma - rho_active;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat1*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat1 = (gamma-Rho);
|
||||
x_d_CAMP = x_LASSO + A'*(y - A*x_LASSO)/Q_hat1;
|
||||
sigma_CAMP = 1/sqrt(log(2))*median(abs(x_d_CAMP));
|
||||
stat_CAMP = abs(x_d_CAMP / sigma_CAMP).^2;
|
||||
|
||||
P_fa_CAMP_cnt(cnt_SNR, rep) = sum(stat_CAMP(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_CAMP_cnt(cnt_SNR, rep) = sum(stat_CAMP(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% SDL
|
||||
Q_hat2 = (gamma - rho_active);
|
||||
x_d_SDL = x_LASSO + A'*(y - A*x_LASSO)/Q_hat2;
|
||||
sigma_SDL = sqrt(gamma)/sqrt(log(2))/(gamma - rho_active)*median(abs(y - A*x_LASSO));
|
||||
stat_SDL = abs(x_d_SDL / sigma_SDL).^2;
|
||||
|
||||
P_fa_SDL_cnt(cnt_SNR, rep) = sum(stat_SDL(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_SDL_cnt(cnt_SNR, rep) = sum(stat_SDL(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% ROD
|
||||
Q_hat3 = (gamma - rho_active)/(1 - rho_active);
|
||||
x_d_ROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat3;
|
||||
chi = rho_active*(1 - rho_active)/(gamma - rho_active);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_ROD = sqrt(2*chi_hat2) / Q_hat3;
|
||||
stat_ROD = abs(x_d_ROD / sigma_ROD).^2;
|
||||
|
||||
P_fa_ROD_cnt(cnt_SNR, rep) = sum(stat_ROD(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_ROD_cnt(cnt_SNR, rep) = sum(stat_ROD(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% LASSO
|
||||
% stat_LASSO = abs(x_LASSO).^2;
|
||||
% if p0 ~= 0
|
||||
% stat_H1_LASSO_cnt(rep, :) = stat_LASSO(x_idx <= thd);
|
||||
% end
|
||||
% stat_H0_LASSO_cnt(rep, :) = stat_LASSO(x_idx > thd);
|
||||
|
||||
|
||||
end
|
||||
|
||||
fprintf('%d\n', rep);
|
||||
|
||||
end
|
||||
|
||||
P_fa_CROD = mean(P_fa_CROD_cnt, 2);
|
||||
P_fa_CAMP = mean(P_fa_CAMP_cnt, 2);
|
||||
P_fa_SDL = mean(P_fa_SDL_cnt, 2);
|
||||
P_fa_ROD = mean(P_fa_ROD_cnt, 2);
|
||||
|
||||
P_d_CROD = mean(P_d_CROD_cnt, 2);
|
||||
P_d_CAMP = mean(P_d_CAMP_cnt, 2);
|
||||
P_d_SDL = mean(P_d_SDL_cnt, 2);
|
||||
P_d_ROD = mean(P_d_ROD_cnt, 2);
|
||||
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
plot(SNR, P_fa_CROD, 'linewidth', 2);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(SNR, P_fa_CAMP, 'linewidth', 2);
|
||||
plot(SNR, P_fa_SDL, 'linewidth', 2);
|
||||
plot(SNR, P_fa_ROD, 'linewidth', 2);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{fa}');
|
||||
|
||||
|
||||
figure(2);
|
||||
plot(SNR, P_d_CROD, 'linewidth', 2);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(SNR, P_d_CAMP, 'linewidth', 2);
|
||||
plot(SNR, P_d_SDL, 'linewidth', 2);
|
||||
plot(SNR, P_d_ROD, 'linewidth', 2);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{d}');
|
||||
|
||||
|
||||
save test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.mat ...
|
||||
SNR...
|
||||
P_fa_CROD...
|
||||
P_fa_CAMP...
|
||||
P_fa_SDL...
|
||||
P_fa_ROD...
|
||||
P_d_CROD...
|
||||
P_d_CAMP...
|
||||
P_d_SDL...
|
||||
P_d_ROD;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,27 @@
|
||||
function [z] = FISTA(y, A, lambda, delta)
|
||||
|
||||
x_pre = A'*y;
|
||||
t = 1;
|
||||
z = x_pre;
|
||||
z_pre = z;
|
||||
t_pre = t;
|
||||
N = size(A, 2);
|
||||
diff = 1;
|
||||
E = eig(A'*A);
|
||||
L = E(end);
|
||||
temp1 = A'*y/L;
|
||||
temp2 = eye(N) - A'*A/L;
|
||||
k = 0;
|
||||
while((diff > delta) && (k < 1000))
|
||||
temp = temp1 + temp2 * z_pre;
|
||||
x = sft_thd(temp, lambda/L);
|
||||
t = 0.5*(1 + sqrt(1+4*t_pre*t_pre));
|
||||
z = x + (x - x_pre) * (t_pre-1) / t;
|
||||
diff = mean(abs(z_pre - z));
|
||||
x_pre = x;
|
||||
z_pre = z;
|
||||
t_pre = t;
|
||||
k = k + 1;
|
||||
end
|
||||
|
||||
end
|
||||
+66
@@ -0,0 +1,66 @@
|
||||
clear;
|
||||
close all;
|
||||
clc;
|
||||
|
||||
load test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.mat;
|
||||
|
||||
Fontsize = 18;
|
||||
plot_width = 800;
|
||||
plot_height = 600;
|
||||
Linewidth = 2;
|
||||
Markersize = 8;
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
plot(p0_total, P_fa_CROD, '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(p0_total, P_fa_CAMP, '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(p0_total, P_fa_SDL, '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(p0_total, P_fa_ROD, '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('signal density');
|
||||
ylabel('P_{fa}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman');
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
|
||||
figure(2);
|
||||
plot(p0_total, P_d_CROD, '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(p0_total, P_d_CAMP, '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(p0_total, P_d_SDL, '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(p0_total, P_d_ROD, '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('signal density');
|
||||
ylabel('P_{d}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman');
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
+249
@@ -0,0 +1,249 @@
|
||||
clc;
|
||||
clear;
|
||||
close all;
|
||||
|
||||
%% parameter setting
|
||||
m = 128;
|
||||
n = 256;
|
||||
|
||||
SNR = 13;
|
||||
|
||||
rep_time = 1e4;
|
||||
|
||||
P_fa = 1e-2;
|
||||
|
||||
p0_total = 0.02: 0.02: 0.2;
|
||||
len_p0 = length(p0_total);
|
||||
|
||||
lambda = 0.1;
|
||||
|
||||
sigma_0 = 0.1;
|
||||
|
||||
sigma_n = sigma_0;
|
||||
|
||||
%% experiment
|
||||
gamma = m/n;
|
||||
|
||||
P_fa_CROD_cnt = zeros(len_p0, rep_time);
|
||||
P_fa_CAMP_cnt = zeros(len_p0, rep_time);
|
||||
P_fa_SDL_cnt = zeros(len_p0, rep_time);
|
||||
P_fa_ROD_cnt = zeros(len_p0, rep_time);
|
||||
P_fa_LASSO_cnt = zeros(len_p0, rep_time);
|
||||
|
||||
P_d_CROD_cnt = zeros(len_p0, rep_time);
|
||||
P_d_CAMP_cnt = zeros(len_p0, rep_time);
|
||||
P_d_SDL_cnt = zeros(len_p0, rep_time);
|
||||
P_d_ROD_cnt = zeros(len_p0, rep_time);
|
||||
P_d_LASSO_cnt = zeros(len_p0, rep_time);
|
||||
|
||||
h_thd = -log(P_fa);
|
||||
|
||||
parfor rep = 1: rep_time
|
||||
|
||||
A_idx = randperm(n);
|
||||
A_idx = A_idx(1: m);
|
||||
A_idx = sort(A_idx);
|
||||
A = dftmtx(n);
|
||||
A = A(A_idx, :);
|
||||
A = A / sqrt(n);
|
||||
|
||||
w = random('Normal', 0, sigma_0/sqrt(2), m, 1) + 1j * random('Normal', 0, sigma_0/sqrt(2), m, 1);
|
||||
% w = w / sqrt(10^(SNR(cnt_SNR)/10));
|
||||
|
||||
for cnt_p0 = 1: len_p0
|
||||
|
||||
p0 = p0_total(cnt_p0);
|
||||
|
||||
x_idx = rand(n, 1);
|
||||
if p0 == 0
|
||||
thd = -1;
|
||||
x_l0 = sum(x_idx > thd);
|
||||
else
|
||||
thd = sort(x_idx);
|
||||
thd = thd(round(n*p0));
|
||||
x_l1 = sum(x_idx <= thd);
|
||||
x_l0 = sum(x_idx > thd);
|
||||
end
|
||||
|
||||
x = zeros(n, 1);
|
||||
x_temp = sigma_0 * exp(1j*random('Uniform', 0, 2*pi, n, 1));
|
||||
x(x_idx <= thd) = x_temp(x_idx <= thd);
|
||||
x = x * sqrt(n/m);
|
||||
|
||||
x1 = x * sqrt(10^(SNR/10));
|
||||
|
||||
y = A * x1 + w;
|
||||
|
||||
% x_LASSO = LASSO_cvx(y, A, lambda);
|
||||
x_LASSO = FISTA(y, A, lambda, 1e-5);
|
||||
|
||||
% CROD
|
||||
rho_active = sum(abs(x_LASSO) > 1e-3)/n;
|
||||
Q_hat = (gamma - rho_active)/(1 - rho_active);
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat = (gamma-Rho)/(1-Rho);
|
||||
x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat;
|
||||
RSS = sum(abs(y - A * x_LASSO).^2)/m;
|
||||
chi = Rho*(1 - Rho)/(gamma - Rho);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_CROD = sqrt(2*chi_hat) / Q_hat;
|
||||
stat_CROD = abs(x_d_CROD / sigma_CROD).^2;
|
||||
|
||||
P_fa_CROD_cnt(cnt_p0, rep) = sum(stat_CROD(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_CROD_cnt(cnt_p0, rep) = sum(stat_CROD(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% CAMP
|
||||
Q_hat1 = gamma - rho_active;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat1*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat1 = (gamma-Rho);
|
||||
x_d_CAMP = x_LASSO + A'*(y - A*x_LASSO)/Q_hat1;
|
||||
sigma_CAMP = 1/sqrt(log(2))*median(abs(x_d_CAMP));
|
||||
stat_CAMP = abs(x_d_CAMP / sigma_CAMP).^2;
|
||||
|
||||
P_fa_CAMP_cnt(cnt_p0, rep) = sum(stat_CAMP(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_CAMP_cnt(cnt_p0, rep) = sum(stat_CAMP(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% SDL
|
||||
Q_hat2 = (gamma - rho_active);
|
||||
x_d_SDL = x_LASSO + A'*(y - A*x_LASSO)/Q_hat2;
|
||||
sigma_SDL = sqrt(gamma)/sqrt(log(2))/(gamma - rho_active)*median(abs(y - A*x_LASSO));
|
||||
stat_SDL = abs(x_d_SDL / sigma_SDL).^2;
|
||||
|
||||
P_fa_SDL_cnt(cnt_p0, rep) = sum(stat_SDL(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_SDL_cnt(cnt_p0, rep) = sum(stat_SDL(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% ROD
|
||||
Q_hat3 = (gamma - rho_active)/(1 - rho_active);
|
||||
x_d_ROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat3;
|
||||
chi = rho_active*(1 - rho_active)/(gamma - rho_active);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_ROD = sqrt(2*chi_hat2) / Q_hat3;
|
||||
stat_ROD = abs(x_d_ROD / sigma_ROD).^2;
|
||||
|
||||
P_fa_ROD_cnt(cnt_p0, rep) = sum(stat_ROD(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_ROD_cnt(cnt_p0, rep) = sum(stat_ROD(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% LASSO
|
||||
% stat_LASSO = abs(x_LASSO).^2;
|
||||
% if p0 ~= 0
|
||||
% stat_H1_LASSO_cnt(rep, :) = stat_LASSO(x_idx <= thd);
|
||||
% end
|
||||
% stat_H0_LASSO_cnt(rep, :) = stat_LASSO(x_idx > thd);
|
||||
|
||||
|
||||
end
|
||||
|
||||
fprintf('%d\n', rep);
|
||||
|
||||
end
|
||||
|
||||
P_fa_CROD = mean(P_fa_CROD_cnt, 2);
|
||||
P_fa_CAMP = mean(P_fa_CAMP_cnt, 2);
|
||||
P_fa_SDL = mean(P_fa_SDL_cnt, 2);
|
||||
P_fa_ROD = mean(P_fa_ROD_cnt, 2);
|
||||
|
||||
P_d_CROD = mean(P_d_CROD_cnt, 2);
|
||||
P_d_CAMP = mean(P_d_CAMP_cnt, 2);
|
||||
P_d_SDL = mean(P_d_SDL_cnt, 2);
|
||||
P_d_ROD = mean(P_d_ROD_cnt, 2);
|
||||
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
plot(p0_total, P_fa_CROD, 'linewidth', 2);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(p0_total, P_fa_CAMP, 'linewidth', 2);
|
||||
plot(p0_total, P_fa_SDL, 'linewidth', 2);
|
||||
plot(p0_total, P_fa_ROD, 'linewidth', 2);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('signal density');
|
||||
ylabel('P_{fa}');
|
||||
|
||||
|
||||
figure(2);
|
||||
plot(p0_total, P_d_CROD, 'linewidth', 2);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(p0_total, P_d_CAMP, 'linewidth', 2);
|
||||
plot(p0_total, P_d_SDL, 'linewidth', 2);
|
||||
plot(p0_total, P_d_ROD, 'linewidth', 2);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('signal density');
|
||||
ylabel('P_{d}');
|
||||
|
||||
|
||||
save test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.mat ...
|
||||
p0_total...
|
||||
P_fa_CROD...
|
||||
P_fa_CAMP...
|
||||
P_fa_SDL...
|
||||
P_fa_ROD...
|
||||
P_d_CROD...
|
||||
P_d_CAMP...
|
||||
P_d_SDL...
|
||||
P_d_ROD;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,27 @@
|
||||
function [z] = FISTA(y, A, lambda, delta)
|
||||
|
||||
x_pre = A'*y;
|
||||
t = 1;
|
||||
z = x_pre;
|
||||
z_pre = z;
|
||||
t_pre = t;
|
||||
N = size(A, 2);
|
||||
diff = 1;
|
||||
E = eig(A'*A);
|
||||
L = E(end);
|
||||
temp1 = A'*y/L;
|
||||
temp2 = eye(N) - A'*A/L;
|
||||
k = 0;
|
||||
while((diff > delta) && (k < 1000))
|
||||
temp = temp1 + temp2 * z_pre;
|
||||
x = sft_thd(temp, lambda/L);
|
||||
t = 0.5*(1 + sqrt(1+4*t_pre*t_pre));
|
||||
z = x + (x - x_pre) * (t_pre-1) / t;
|
||||
diff = mean(abs(z_pre - z));
|
||||
x_pre = x;
|
||||
z_pre = z;
|
||||
t_pre = t;
|
||||
k = k + 1;
|
||||
end
|
||||
|
||||
end
|
||||
+66
@@ -0,0 +1,66 @@
|
||||
clear;
|
||||
close all;
|
||||
clc;
|
||||
|
||||
load test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.mat;
|
||||
|
||||
Fontsize = 18;
|
||||
plot_width = 800;
|
||||
plot_height = 600;
|
||||
Linewidth = 2;
|
||||
Markersize = 8;
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
plot(gamma_total, P_fa_CROD, '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(gamma_total, P_fa_CAMP, '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(gamma_total, P_fa_SDL, '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(gamma_total, P_fa_ROD, '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('compression rate');
|
||||
ylabel('P_{fa}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman');
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
|
||||
figure(2);
|
||||
plot(gamma_total, P_d_CROD, '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(gamma_total, P_d_CAMP, '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(gamma_total, P_d_SDL, '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(gamma_total, P_d_ROD, '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('compression rate');
|
||||
ylabel('P_{d}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
%set(gca, 'FontSize', Fontsize, 'fontname', 'Times New Roman');
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
+251
@@ -0,0 +1,251 @@
|
||||
clc;
|
||||
clear;
|
||||
close all;
|
||||
|
||||
%% parameter setting
|
||||
n = 256;
|
||||
|
||||
SNR = 13;
|
||||
|
||||
rep_time = 1e4;
|
||||
|
||||
P_fa = 1e-2;
|
||||
|
||||
p0 = 0.1;
|
||||
|
||||
gamma_total = (4: 12)/16;
|
||||
len_gamma = length(gamma_total);
|
||||
|
||||
lambda = 0.1;
|
||||
|
||||
sigma_0 = 0.1;
|
||||
|
||||
sigma_n = sigma_0;
|
||||
|
||||
%% experiment
|
||||
|
||||
P_fa_CROD_cnt = zeros(len_gamma, rep_time);
|
||||
P_fa_CAMP_cnt = zeros(len_gamma, rep_time);
|
||||
P_fa_SDL_cnt = zeros(len_gamma, rep_time);
|
||||
P_fa_ROD_cnt = zeros(len_gamma, rep_time);
|
||||
P_fa_LASSO_cnt = zeros(len_gamma, rep_time);
|
||||
|
||||
P_d_CROD_cnt = zeros(len_gamma, rep_time);
|
||||
P_d_CAMP_cnt = zeros(len_gamma, rep_time);
|
||||
P_d_SDL_cnt = zeros(len_gamma, rep_time);
|
||||
P_d_ROD_cnt = zeros(len_gamma, rep_time);
|
||||
P_d_LASSO_cnt = zeros(len_gamma, rep_time);
|
||||
|
||||
h_thd = -log(P_fa);
|
||||
|
||||
parfor rep = 1: rep_time
|
||||
|
||||
|
||||
|
||||
for cnt_gamma = 1: len_gamma
|
||||
|
||||
gamma = gamma_total(cnt_gamma);
|
||||
m = round(gamma*n);
|
||||
|
||||
A_idx = randperm(n);
|
||||
A_idx = A_idx(1: m);
|
||||
A_idx = sort(A_idx);
|
||||
A = dftmtx(n);
|
||||
A = A(A_idx, :);
|
||||
A = A / sqrt(n);
|
||||
|
||||
w = random('Normal', 0, sigma_0/sqrt(2), m, 1) + 1j * random('Normal', 0, sigma_0/sqrt(2), m, 1);
|
||||
|
||||
x_idx = rand(n, 1);
|
||||
if p0 == 0
|
||||
thd = -1;
|
||||
x_l0 = sum(x_idx > thd);
|
||||
else
|
||||
thd = sort(x_idx);
|
||||
thd = thd(round(n*p0));
|
||||
x_l1 = sum(x_idx <= thd);
|
||||
x_l0 = sum(x_idx > thd);
|
||||
end
|
||||
|
||||
x = zeros(n, 1);
|
||||
x_temp = sigma_0 * exp(1j*random('Uniform', 0, 2*pi, n, 1));
|
||||
x(x_idx <= thd) = x_temp(x_idx <= thd);
|
||||
x = x * sqrt(n/m);
|
||||
|
||||
x1 = x * sqrt(10^(SNR/10));
|
||||
|
||||
y = A * x1 + w;
|
||||
|
||||
% x_LASSO = LASSO_cvx(y, A, lambda);
|
||||
x_LASSO = FISTA(y, A, lambda, 1e-5);
|
||||
|
||||
% CROD
|
||||
rho_active = sum(abs(x_LASSO) > 1e-3)/n;
|
||||
Q_hat = (gamma - rho_active)/(1 - rho_active);
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat = (gamma-Rho)/(1-Rho);
|
||||
x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat;
|
||||
RSS = sum(abs(y - A * x_LASSO).^2)/m;
|
||||
chi = Rho*(1 - Rho)/(gamma - Rho);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_CROD = sqrt(2*chi_hat) / Q_hat;
|
||||
stat_CROD = abs(x_d_CROD / sigma_CROD).^2;
|
||||
|
||||
P_fa_CROD_cnt(cnt_gamma, rep) = sum(stat_CROD(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_CROD_cnt(cnt_gamma, rep) = sum(stat_CROD(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% CAMP
|
||||
Q_hat1 = gamma - rho_active;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat1*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat1 = (gamma-Rho);
|
||||
x_d_CAMP = x_LASSO + A'*(y - A*x_LASSO)/Q_hat1;
|
||||
sigma_CAMP = 1/sqrt(log(2))*median(abs(x_d_CAMP));
|
||||
stat_CAMP = abs(x_d_CAMP / sigma_CAMP).^2;
|
||||
|
||||
P_fa_CAMP_cnt(cnt_gamma, rep) = sum(stat_CAMP(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_CAMP_cnt(cnt_gamma, rep) = sum(stat_CAMP(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% SDL
|
||||
Q_hat2 = (gamma - rho_active);
|
||||
x_d_SDL = x_LASSO + A'*(y - A*x_LASSO)/Q_hat2;
|
||||
sigma_SDL = sqrt(gamma)/sqrt(log(2))/(gamma - rho_active)*median(abs(y - A*x_LASSO));
|
||||
stat_SDL = abs(x_d_SDL / sigma_SDL).^2;
|
||||
|
||||
P_fa_SDL_cnt(cnt_gamma, rep) = sum(stat_SDL(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_SDL_cnt(cnt_gamma, rep) = sum(stat_SDL(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% ROD
|
||||
Q_hat3 = (gamma - rho_active)/(1 - rho_active);
|
||||
x_d_ROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat3;
|
||||
chi = rho_active*(1 - rho_active)/(gamma - rho_active);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat2 = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_ROD = sqrt(2*chi_hat2) / Q_hat3;
|
||||
stat_ROD = abs(x_d_ROD / sigma_ROD).^2;
|
||||
|
||||
P_fa_ROD_cnt(cnt_gamma, rep) = sum(stat_ROD(x_idx > thd) > h_thd) / x_l0;
|
||||
P_d_ROD_cnt(cnt_gamma, rep) = sum(stat_ROD(x_idx <= thd) > h_thd) / x_l1;
|
||||
|
||||
% LASSO
|
||||
% stat_LASSO = abs(x_LASSO).^2;
|
||||
% if p0 ~= 0
|
||||
% stat_H1_LASSO_cnt(rep, :) = stat_LASSO(x_idx <= thd);
|
||||
% end
|
||||
% stat_H0_LASSO_cnt(rep, :) = stat_LASSO(x_idx > thd);
|
||||
|
||||
|
||||
end
|
||||
|
||||
fprintf('%d\n', rep);
|
||||
|
||||
end
|
||||
|
||||
P_fa_CROD = mean(P_fa_CROD_cnt, 2);
|
||||
P_fa_CAMP = mean(P_fa_CAMP_cnt, 2);
|
||||
P_fa_SDL = mean(P_fa_SDL_cnt, 2);
|
||||
P_fa_ROD = mean(P_fa_ROD_cnt, 2);
|
||||
|
||||
P_d_CROD = mean(P_d_CROD_cnt, 2);
|
||||
P_d_CAMP = mean(P_d_CAMP_cnt, 2);
|
||||
P_d_SDL = mean(P_d_SDL_cnt, 2);
|
||||
P_d_ROD = mean(P_d_ROD_cnt, 2);
|
||||
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
plot(gamma_total, P_fa_CROD, 'linewidth', 2);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(gamma_total, P_fa_CAMP, 'linewidth', 2);
|
||||
plot(gamma_total, P_fa_SDL, 'linewidth', 2);
|
||||
plot(gamma_total, P_fa_ROD, 'linewidth', 2);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('compression rate');
|
||||
ylabel('P_{fa}');
|
||||
|
||||
|
||||
figure(2);
|
||||
plot(gamma_total, P_d_CROD, 'linewidth', 2);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(gamma_total, P_d_CAMP, 'linewidth', 2);
|
||||
plot(gamma_total, P_d_SDL, 'linewidth', 2);
|
||||
plot(gamma_total, P_d_ROD, 'linewidth', 2);
|
||||
legend('CROD', 'CAMP', 'SDL-test', 'ROD');
|
||||
xlabel('compression rate');
|
||||
ylabel('P_{d}');
|
||||
|
||||
|
||||
save test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.mat ...
|
||||
gamma_total...
|
||||
P_fa_CROD...
|
||||
P_fa_CAMP...
|
||||
P_fa_SDL...
|
||||
P_fa_ROD...
|
||||
P_d_CROD...
|
||||
P_d_CAMP...
|
||||
P_d_SDL...
|
||||
P_d_ROD;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,82 @@
|
||||
% Input: y,A,lambda,tau,Kit
|
||||
% Output: x_hat_wl,x_hat_d
|
||||
|
||||
% Main structure of cVAMP
|
||||
function [x_hat_wl, x_hat_d] = cVAMPro(y, A, lambda, tau, Kit)
|
||||
|
||||
% Initialization
|
||||
[M, N] = size(A);
|
||||
gamma = M / N;
|
||||
k = 0;
|
||||
p = ctranspose(A) * y;
|
||||
h_1 = p;
|
||||
Q_1 = gamma;
|
||||
tau_d = 1;
|
||||
|
||||
% Iteration
|
||||
while ((k < Kit) && (tau_d > tau))
|
||||
% Factorized Part
|
||||
x_1 = ST(h_1, lambda, Q_1);
|
||||
chi_1 = F1(x_1, lambda, Q_1);
|
||||
% Message Passing
|
||||
h_2 = x_1 / chi_1 - h_1;
|
||||
Q_2 = 1 / chi_1 - Q_1;
|
||||
% Gaussian Part
|
||||
t1 = (p + h_2) / Q_2;
|
||||
t2 = ctranspose(A) * (A * (p + h_2)) / ((Q_2 + 1) * Q_2);
|
||||
x_2 = t1 - t2;
|
||||
chi_2 = gamma / (Q_2 + 1) + (1 - gamma) / Q_2;
|
||||
% Message Passing
|
||||
h_1_next = x_2 ./ chi_2 - h_2;
|
||||
Q_1_next = 1 / chi_2 - Q_2;
|
||||
tau_d = norm(h_1_next - h_1, Inf) / norm(h_1_next, Inf);
|
||||
k = k + 1;
|
||||
% output
|
||||
x_hat_wl = x_1;
|
||||
x_hat_d = h_1_next / Q_1_next;
|
||||
% next
|
||||
h_1 = h_1_next;
|
||||
Q_1 = Q_1_next;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
% SoftThreshold function
|
||||
function x = ST(h_1, lambda, Q_1)
|
||||
[N, M] = size(h_1);
|
||||
x = zeros(N, M);
|
||||
|
||||
for i = 1:N
|
||||
sign = h_1(i) ./ abs(h_1(i));
|
||||
diff = abs(h_1(i)) - lambda(i);
|
||||
x(i) = sign .* (diff ./ Q_1) .* SF(diff);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
% Heaviside's step function
|
||||
function v = SF(a)
|
||||
|
||||
if a > 0
|
||||
v = 1;
|
||||
elseif a == 0
|
||||
v = 0; % at zero points
|
||||
else
|
||||
v = 0;
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
% Calculation of chi_1
|
||||
function chi_1 = F1(x_1, lambda, Q_1)
|
||||
[N, M] = size(x_1);
|
||||
count = 0;
|
||||
|
||||
for i = 1:N
|
||||
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) * (abs(x_1(i)) > 1e-4);
|
||||
end
|
||||
|
||||
chi_1 = count / (2 * N * Q_1);
|
||||
end
|
||||
@@ -0,0 +1,122 @@
|
||||
clear;
|
||||
close all;
|
||||
clc;
|
||||
|
||||
load test_MC_Pfa_Rmf_Dcs.mat;
|
||||
|
||||
Fontsize = 18;
|
||||
plot_width = 800;
|
||||
plot_height = 600;
|
||||
Linewidth = 2;
|
||||
Markersize = 8;
|
||||
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
loglog(P_fa, P_fa_CROD(1,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
loglog(P_fa, P_fa_CROD(3,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_CROD(5,:), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_CROD(7,:), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_CROD(9,:), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('SNR = -30', 'SNR = -20', 'SNR = -10', 'SNR = 0', 'SNR = 10');
|
||||
xlabel('P_{fa} Set');
|
||||
ylabel('Actual P_{fa}');
|
||||
ylim([8e-5,1]);
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('MF P_{fa}')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(3);
|
||||
semilogx(P_fa_CROD(1,:), P_d_CROD(1,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
semilogx(P_fa_CROD(3,:), P_d_CROD(3,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_CROD(5,:), P_d_CROD(5,:), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_CROD(7,:), P_d_CROD(7,:), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_CROD(9,:), P_d_CROD(9,:), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('SNR = -30', 'SNR = -20', 'SNR = -10', 'SNR = 0', 'SNR = 10');
|
||||
xlabel('P_{fa}');
|
||||
ylabel('P_{d}');
|
||||
xlim([8e-5,1]);
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('MF ROC')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(5);
|
||||
semilogy(SNR, P_fa_CROD(:,1), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
semilogy(SNR, P_fa_CROD(:,3), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_CROD(:,5), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_CROD(:,7), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_CROD(:,9), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('Pfa = 1e-4', 'Pfa = 1e-3', 'Pfa = 1e-2', 'Pfa = 1e-1', 'Pfa = 1');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{fa}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('MF')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(6);
|
||||
plot(SNR, P_d_CROD(:,1), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(SNR, P_d_CROD(:,3), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_CROD(:,5), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_CROD(:,7), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_CROD(:,9), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('Pfa = 1e-4', 'Pfa = 1e-3', 'Pfa = 1e-2', 'Pfa = 1e-1', 'Pfa = 1');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{d}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('MF')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
@@ -0,0 +1,119 @@
|
||||
clear;
|
||||
close all;
|
||||
clc;
|
||||
|
||||
load test_MC_Pfa_Rmf_Dmf.mat;
|
||||
|
||||
Fontsize = 18;
|
||||
plot_width = 800;
|
||||
plot_height = 600;
|
||||
Linewidth = 2;
|
||||
Markersize = 8;
|
||||
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
loglog(P_fa, P_fa_MF(1,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
loglog(P_fa, P_fa_MF(3,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_MF(5,:), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_MF(7,:), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('SNR = -30', 'SNR = -20', 'SNR = -10', 'SNR = 0');
|
||||
xlabel('P_{fa} Set');
|
||||
ylabel('Actual P_{fa}');
|
||||
ylim([8e-5,1]);
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('MF P_{fa}')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(3);
|
||||
semilogx(P_fa_MF(1,:), P_d_MF(1,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
semilogx(P_fa_MF(3,:), P_d_MF(3,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_MF(5,:), P_d_MF(5,:), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_MF(7,:), P_d_MF(7,:), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_MF(9,:), P_d_MF(9,:), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('SNR = -30', 'SNR = -20', 'SNR = -10', 'SNR = 0', 'SNR = 0');
|
||||
xlabel('P_{fa}');
|
||||
ylabel('P_{d}');
|
||||
xlim([8e-5,1]);
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('MF ROC')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(5);
|
||||
semilogy(SNR, P_fa_MF(:,1), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
semilogy(SNR, P_fa_MF(:,3), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_MF(:,5), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_MF(:,7), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_MF(:,9), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('Pfa = 1e-4', 'Pfa = 1e-3', 'Pfa = 1e-2', 'Pfa = 1e-1', 'Pfa = 1');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{fa}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('MF')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(6);
|
||||
plot(SNR, P_d_MF(:,1), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(SNR, P_d_MF(:,3), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_MF(:,5), '-d', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_MF(:,7), '-s', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_MF(:,9), '-^', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('Pfa = 1e-4', 'Pfa = 1e-3', 'Pfa = 1e-2', 'Pfa = 1e-1', 'Pfa = 1');
|
||||
xlabel('SNR');
|
||||
ylabel('P_{d}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('MF')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
@@ -0,0 +1,147 @@
|
||||
clear;
|
||||
close all;
|
||||
clc;
|
||||
|
||||
load test_MC_Pfa_Rmf_Dmf.mat;
|
||||
load test_MC_Pfa_Rmf_Dcs.mat;
|
||||
|
||||
Fontsize = 18;
|
||||
plot_width = 800;
|
||||
plot_height = 600;
|
||||
Linewidth = 2;
|
||||
Markersize = 8;
|
||||
% SNR = [0,5,8,10,12,15,20,25];
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
loglog(P_fa, P_fa_MF(9,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
loglog(P_fa, P_fa_CROD(9,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_MF(1,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
loglog(P_fa, P_fa_CROD(1,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('MFMF,SNR = 10dB', 'MFCS,SNR = 10dB','MFMF,SNR = -20dB', 'MFCS,SNR = -20dB');
|
||||
xlabel('P_{fa} Set');
|
||||
ylabel('Actual P_{fa}');
|
||||
ylim([8e-5,1]);
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('P_{fa}')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
figure(3);
|
||||
semilogx(P_fa_CROD(7,:), P_d_CROD(7,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
semilogx(P_fa_MF(7,:), P_d_MF(7,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_CROD(5,:), P_d_CROD(5,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_MF(5,:), P_d_MF(5,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_CROD(1,:), P_d_CROD(1,:), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogx(P_fa_MF(1,:), P_d_MF(1,:), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('MFCS,SNR = 15dB', 'MFMF,SNR = 15dB', 'MFCS,SNR = 10dB', 'MFMF,SNR = 10dB', 'MFCS,SNR = -5dB', 'MFMF,SNR = -5dB');
|
||||
xlabel('P_{fa}');
|
||||
ylabel('P_{d}');
|
||||
xlim([8e-5,1]);
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('ROC')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
% figure(3);
|
||||
% semilogx(P_fa_CROD(9,:), P_d_CROD(9,:), '-o', ...
|
||||
% 'Linewidth', Linewidth, ...
|
||||
% 'MarkerSize', Markersize);
|
||||
% hold on;
|
||||
% grid on;
|
||||
% semilogx(P_fa_MF(9,:), P_d_MF(9,:), '-+', ...
|
||||
% 'Linewidth', Linewidth, ...
|
||||
% 'MarkerSize', Markersize);
|
||||
% legend('MFCS,SNR = 10dB', 'MFMF,SNR = 10dB');
|
||||
% xlabel('P_{fa}');
|
||||
% ylabel('P_{d}');
|
||||
% % xlim([8e-5,1]);
|
||||
% set(gca, 'FontSize', Fontsize);
|
||||
% title('ROC')
|
||||
% set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
% set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
figure(5);
|
||||
semilogy(SNR, P_fa_MF(:,5), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
semilogy(SNR, P_fa_CROD(:,5), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_MF(:,3), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_CROD(:,3), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_MF(:,1), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
semilogy(SNR, P_fa_CROD(:,1), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
legend('MFMF,Pfa = 1e-2', 'MFCS,Pfa = 1e-2', 'MFMF,Pfa = 1e-3', 'MFCS,Pfa = 1e-3', 'MFMF,Pfa = 1e-4', 'MFCS,Pfa = 1e-4');
|
||||
xlabel('SNR/dB');
|
||||
ylabel('P_{fa}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('Actual P_{fa}')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
|
||||
figure(7);
|
||||
plot(SNR, P_d_MF(:,5), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
hold on;
|
||||
grid on;
|
||||
plot(SNR, P_d_CROD(:,5), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_MF(:,3), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_CROD(:,3), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_MF(:,1), '-+', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
plot(SNR, P_d_CROD(:,1), '-o', ...
|
||||
'Linewidth', Linewidth, ...
|
||||
'MarkerSize', Markersize);
|
||||
|
||||
legend('MFMF,Pfa = 1e-2', 'MFCS,Pfa = 1e-2', 'MFMF,Pfa = 1e-3', 'MFCS,Pfa = 1e-3', 'MFMF,Pfa = 1e-4', 'MFCS,Pfa = 1e-4');
|
||||
xlabel('SNR/dB');
|
||||
ylabel('P_{d}');
|
||||
set(gca, 'FontSize', Fontsize);
|
||||
title('Actual P_{d}')
|
||||
set(gcf, 'position', [200, 300, plot_width, plot_height]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
@@ -0,0 +1,11 @@
|
||||
test_Pfa_noT_Rmf_Dmf.m 距离匹配滤波，多普勒匹配滤波【RD图，分布验证，目标场景图】
|
||||
test_Pfa_noT_Rmf_Dcs.m 距离匹配滤波，多普勒压缩感知【RD图，分布验证，目标场景图】
|
||||
|
||||
test_MC_Pfa_Rmf_Dmf.m 距离匹配滤波，多普勒匹配滤波【蒙特卡洛实验】
|
||||
test_MC_Pfa_Rmf_Dmf.mat
|
||||
test_MC_Pfa_Rmf_Dcs.m 距离匹配滤波，多普勒压缩感知【蒙特卡洛实验】
|
||||
test_MC_Pfa_Rmf_Dcs.mat
|
||||
|
||||
plot_MC_Pfa_Rmf_Dmf.m 距离匹配滤波，多普勒匹配滤波【蒙特卡洛结果绘制】
|
||||
plot_MC_Pfa_Rmf_Dcs.m 距离匹配滤波，多普勒压缩感知【蒙特卡洛结果绘制】
|
||||
plot_Pd_Pfa_diff.m 对比图绘制
|
||||
@@ -0,0 +1,323 @@
|
||||
clc
|
||||
clear
|
||||
close all;
|
||||
|
||||
%% 蜿よ焚隶ｾ鄂ｮ
|
||||
B = 5e6; % 菫｡蜿ｷ蟶ｦ�??5MHz
|
||||
Tp = 20e-6; % 閼牙??100us
|
||||
fs = 2 * B; % 驥��ｷ鬚醍紫
|
||||
Ts = 1 / fs; % 驥��ｷ蜻ｨ譛�
|
||||
K = B / Tp; % 郤ｿ諤ｧ隹�｢醍紫
|
||||
fc = 1.25e9; % 霓ｽ豕｢鬚醍紫1.25GHz
|
||||
PRF = 5000; % 閼牙�驥榊?埼?醍紫
|
||||
Tr = 1 / PRF; % 閼牙�驥榊?埼龍髫?
|
||||
|
||||
|
||||
|
||||
numP = 64; % 閼牙�謨ｰ驥�
|
||||
|
||||
t = 0: 1 / fs: Tr - 1 / fs;
|
||||
tP = 0: 1 / fs: Tr * numP - 1 / fs;
|
||||
|
||||
c = 3e8; % 蜈蛾?
|
||||
|
||||
slow_len = 128;
|
||||
|
||||
%%
|
||||
sigma_n = 0.1;
|
||||
alpha_prop = 0.8;
|
||||
|
||||
% SNR = -30;
|
||||
SNR = [-5,0,5,8,10,12,15,20,25];
|
||||
len_SNR = length(SNR);
|
||||
|
||||
P_fa = [1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1];
|
||||
len_P_fa = length(P_fa);
|
||||
|
||||
rep_time = 200;
|
||||
|
||||
% rou >= c / 2B = 15(m)
|
||||
|
||||
%%
|
||||
delta_VAMP = 1e-6;
|
||||
iter_max = 1000;
|
||||
lambda_val = 0.05;
|
||||
lambda = zeros(slow_len,1) + lambda_val;
|
||||
|
||||
%% 莠ｧ逕溷書蟆�ｿ｡蜿ｷ
|
||||
N = Tr * fs;
|
||||
N_high = Tp * fs;
|
||||
|
||||
signal_t = zeros(1, N * numP);
|
||||
for j = 1: numP
|
||||
for i = 1: N_high
|
||||
tp = ((j - 1) * N + i) * (1 / fs);
|
||||
signal_t(1, i + (j-1)*N) = exp(1j * 2 * pi * (fc * tp + 0.5 * K * tp .^ 2));
|
||||
end
|
||||
end
|
||||
|
||||
multiple_r = signal_t(1, 1: N) * signal_t(1, 1: N)';
|
||||
% figure(1);
|
||||
% subplot(311)
|
||||
% plot(tP,real(signal_t));
|
||||
% xlabel('譌ｶ髣ｴ/t');ylabel('蟷�ｺｦ');
|
||||
% title('蜿大ｰ�ｿ｡蜿ｷ');
|
||||
|
||||
%% 隶ｾ鄂ｮ逶?譬? 1
|
||||
distance1 = 18000;
|
||||
v1 = 187.5;
|
||||
|
||||
tau1 = distance1 * 2 / c;
|
||||
n_tau1 = round(tau1 * fs);
|
||||
|
||||
f_d1 = 2 * fc * v1 / c;
|
||||
|
||||
alpha1 = alpha_prop;
|
||||
% alpha1 = 0;
|
||||
|
||||
signal_r1 = zeros(1, N * numP);
|
||||
for j = 1: numP
|
||||
for i = 1: N
|
||||
temp = i - n_tau1;
|
||||
if temp >= 1 && temp <= N_high
|
||||
tp = ((j - 1) * N + i - n_tau1) * (1 / fs);
|
||||
signal_r1(1, i + (j-1)*N) = alpha1 * exp(1j*2*pi*((fc - f_d1) * tp + 0.5 * K * tp .^ 2));
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
%% 逕滓�蝗樊ｳ｢
|
||||
signal_r = signal_r1;
|
||||
% figure(1);
|
||||
% subplot(312)
|
||||
% plot(tP, real(signal_r));
|
||||
% xlabel('譌ｶ髣ｴ/t');
|
||||
% ylabel('蟷�ｺｦ');
|
||||
% title('蝗樊ｳ｢菫｡蜿ｷ');
|
||||
|
||||
%%
|
||||
P_fa_CROD_cnt = zeros(len_SNR, len_P_fa, rep_time);
|
||||
P_d_CROD_cnt = zeros(len_SNR, len_P_fa, rep_time);
|
||||
|
||||
F_ori = dftmtx(slow_len);
|
||||
F = F_ori(1:numP,:);
|
||||
multiple_d = F(:,1)' * F(:,1);
|
||||
F_inv = conj(F)/slow_len;
|
||||
|
||||
A = (sqrt(slow_len) * eye(numP)) * F_inv;
|
||||
|
||||
n = slow_len;
|
||||
m = numP;
|
||||
gamma = numP / slow_len;
|
||||
|
||||
parfor rep = 1: rep_time
|
||||
% for rep = 1: rep_time
|
||||
|
||||
for cnt_SNR = 1: len_SNR
|
||||
|
||||
%% 蝎?螢ｰ�?�炊
|
||||
noise = random('Normal', 0, sigma_n/sqrt(2), 1, N * numP) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, N * numP);
|
||||
alpha_SNR = sqrt(10^(SNR(cnt_SNR)/10) * sigma_n^2 / (multiple_r * multiple_d));
|
||||
signal_r_n = signal_r * alpha_SNR + noise;
|
||||
|
||||
% figure(1);
|
||||
% subplot(313)
|
||||
% plot(tP, real(signal_r_n));
|
||||
% xlabel('譌ｶ髣ｴ/t');
|
||||
% ylabel('蟷�ｺｦ');
|
||||
% title('蝗樊ｳ｢+蝎?螢ｰ菫｡蜿?');
|
||||
|
||||
|
||||
%% 蛹ｹ驟肴ｻ､豕｢窶披疲潔Tr蛻貞�
|
||||
Srange = zeros(N, numP);
|
||||
|
||||
% A = generate_matrix(transpose(signal_t(1, 1: N)), 1);
|
||||
|
||||
count_fa=0;
|
||||
|
||||
for i = 1: numP
|
||||
sr = signal_r_n(1, 1 + (i - 1) * N: i * N);
|
||||
st = signal_t(1, 1 + (i - 1) * N: i * N);
|
||||
|
||||
% A = generate_matrix(transpose(st), 1);
|
||||
% mf = A' * transpose(sr);
|
||||
% filterred_rf_r = mf ./ multiple_r;
|
||||
|
||||
% 蛹ｹ驟肴ｻ､豕｢
|
||||
rf_r_fft_conj = conj(fft(st));
|
||||
filterred_rf_r = ifft(fft(sr) .* rf_r_fft_conj);
|
||||
filterred_rf_r = filterred_rf_r ./ multiple_r;
|
||||
|
||||
Srange(:,i) = filterred_rf_r';
|
||||
|
||||
end
|
||||
|
||||
|
||||
% figure(101)
|
||||
% mesh(abs(Srange))
|
||||
% title('謖欝r霑幄?悟源驟肴ｻ､豕｢扈捺�?')
|
||||
|
||||
%% R蛹ｹ驟肴ｻ､豕｢窶披泌�蟶�ｪ瑚ｯ?
|
||||
% Before MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2)
|
||||
% b~N(0,sigma_n^2 / 2)
|
||||
|
||||
% After MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2 / multiple_r)
|
||||
% b~N(0,sigma_n^2 / 2 / multiple_r)
|
||||
% a^2 + b^2 ~ chi^2(2) * sigma_n^2 / 2 / multiple_r
|
||||
|
||||
% stat_R = abs(Srange).^2;
|
||||
% P_fa = 0.1;
|
||||
% kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2 / multiple;
|
||||
% sum(sum(stat_R>kd))/numel(stat_R)
|
||||
|
||||
%% 螟壽勸蜍呈ｻ､豕?
|
||||
Srd = zeros(N, slow_len);
|
||||
|
||||
count_Pfa_d = zeros(length(P_fa), N);
|
||||
count_Pd_d = zeros(length(P_fa), 1);
|
||||
|
||||
target_node_d = round(v1 / (c / fc * PRF / 2 / slow_len)) + 1;
|
||||
target_node_r = n_tau1 + 1;
|
||||
|
||||
% for i = 1: N
|
||||
for i = target_node_r-N_high: target_node_r+N_high
|
||||
x_slow = Srange(i, :);
|
||||
y = (sqrt(slow_len) * eye(numP)) * transpose(x_slow);
|
||||
|
||||
[x_LASSO,y_d] = cVAMPro(y,A,lambda,delta_VAMP,iter_max);
|
||||
|
||||
% CROD豎ょ悉蛛?
|
||||
rho_active = sum(abs(x_LASSO) > 1e-3)/n;
|
||||
Q_hat = (gamma - rho_active)/(1 - rho_active);
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat = (gamma-Rho)/(1-Rho);
|
||||
x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat;
|
||||
|
||||
% CROD豎る葎髯仙柱譽鬪檎ｻ溯ｮ｡驥�
|
||||
RSS = sum(abs(y - A * x_LASSO).^2)/m;
|
||||
chi = Rho*(1 - Rho)/(gamma - Rho);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_CROD = sqrt(2*chi_hat) / Q_hat;
|
||||
|
||||
sigma_CROD = sigma_CROD*0.895;
|
||||
|
||||
stat_CROD = abs(x_d_CROD / sigma_CROD).^2;
|
||||
|
||||
h_thd = chi2inv(1 - P_fa, 2) / 2;
|
||||
|
||||
% 譽豬?
|
||||
for j = 1:length(P_fa)
|
||||
|
||||
if i ~= target_node_r
|
||||
count_Pfa_d(j,i) = sum(stat_CROD > h_thd(j)) / slow_len;
|
||||
else
|
||||
stat_index = ones(size(stat_CROD));
|
||||
stat_index(target_node_d) = 0;
|
||||
count_Pfa_d(j,i) = sum(stat_CROD(stat_index > 0) > h_thd(j)) / sum(stat_index);
|
||||
count_Pd_d(j) = stat_CROD(target_node_d) > h_thd(j);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
Srd(i, :) = fftshift(transpose(x_d_CROD ./ multiple_d ./ 2));
|
||||
|
||||
end
|
||||
|
||||
target_near_count = count_Pfa_d(:,target_node_r-N_high: target_node_r+N_high);
|
||||
% P_fa_actual = mean(count_Pfa_d,2);
|
||||
P_fa_actual = mean(target_near_count,2);
|
||||
|
||||
for cnt_h_th = 1: len_P_fa
|
||||
P_fa_CROD_cnt(cnt_SNR, cnt_h_th, rep) = P_fa_actual(cnt_h_th);
|
||||
P_d_CROD_cnt(cnt_SNR, cnt_h_th, rep) = count_Pd_d(cnt_h_th);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
fprintf('%d\n', rep);
|
||||
|
||||
end
|
||||
|
||||
P_fa_CROD = mean(P_fa_CROD_cnt, 3);
|
||||
P_d_CROD = mean(P_d_CROD_cnt, 3);
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
loglog(P_fa,P_fa_CROD(1,:), 'linewidth', 2);
|
||||
% hold on;
|
||||
% grid on;
|
||||
% loglog(P_fa,P_fa_MF(2,:), 'linewidth', 2);
|
||||
% loglog(P_fa,P_fa_MF(3,:), 'linewidth', 2);
|
||||
% loglog(P_fa,P_fa_MF(4,:), 'linewidth', 2);
|
||||
% loglog(P_fa,P_fa_MF(5,:), 'linewidth', 2);
|
||||
% loglog(P_fa,P_fa_MF(6,:), 'linewidth', 2);
|
||||
% legend('SNR = 0','SNR = 2','SNR = 4','SNR = 6','SNR = 8','SNR = 10');
|
||||
xlabel('P_fa');
|
||||
ylabel('Actual P_fa');
|
||||
|
||||
figure(2);
|
||||
semilogx(P_fa,P_d_CROD(1,:), 'linewidth', 2);
|
||||
% hold on;
|
||||
% grid on;
|
||||
% semilogx(P_fa,P_d_MF(2,:), 'linewidth', 2);
|
||||
% semilogx(P_fa,P_d_MF(3,:), 'linewidth', 2);
|
||||
% semilogx(P_fa,P_d_MF(4,:), 'linewidth', 2);
|
||||
% semilogx(P_fa,P_d_MF(5,:), 'linewidth', 2);
|
||||
% semilogx(P_fa,P_d_MF(6,:), 'linewidth', 2);
|
||||
% legend('SNR = 0','SNR = 2','SNR = 4','SNR = 6','SNR = 8','SNR = 10');
|
||||
xlabel('P_fa');
|
||||
ylabel('P_d');
|
||||
|
||||
% save test_Rmf_Dmf.mat ...
|
||||
% distance_temp...
|
||||
% speed_temp...
|
||||
% Srd...
|
||||
% Srange...
|
||||
% tP...
|
||||
% signal_t...
|
||||
% signal_r...
|
||||
% signal_r_n...
|
||||
% F...
|
||||
% slow_len;
|
||||
|
||||
save test_MC_Pfa_Rmf_Dcs.mat ...
|
||||
SNR...
|
||||
P_fa...
|
||||
P_fa_CROD...
|
||||
P_d_CROD...
|
||||
sigma_n...
|
||||
B...
|
||||
Tp...
|
||||
fs...
|
||||
Ts...
|
||||
K...
|
||||
fc...
|
||||
PRF...
|
||||
Tr...
|
||||
numP...
|
||||
t...
|
||||
tP...
|
||||
slow_len;
|
||||
@@ -0,0 +1,276 @@
|
||||
clc
|
||||
clear
|
||||
close all;
|
||||
|
||||
%% 参数设置
|
||||
B = 5e6; % 信号带宽5MHz
|
||||
Tp = 20e-6; % 脉宽100us
|
||||
fs = 2 * B; % 采样频率
|
||||
Ts = 1 / fs; % 采样周期
|
||||
K = B / Tp; % 线性调频率
|
||||
fc = 1.25e9; % 载波频率1.25GHz
|
||||
PRF = 5000; % 脉冲重复频率
|
||||
Tr = 1 / PRF; % 脉冲重复间隔
|
||||
|
||||
numP = 64; % 脉冲数量
|
||||
|
||||
t = 0: 1 / fs: Tr - 1 / fs;
|
||||
tP = 0: 1 / fs: Tr * numP - 1 / fs;
|
||||
|
||||
c = 3e8; % 光速
|
||||
|
||||
slow_len = 128;
|
||||
|
||||
%%
|
||||
sigma_n = 0.1;
|
||||
alpha_prop = 0.8;
|
||||
|
||||
% SNR = -30;
|
||||
SNR = [-5,0,5,8,10,12,15,20,25];
|
||||
len_SNR = length(SNR);
|
||||
|
||||
P_fa = [1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1];
|
||||
len_P_fa = length(P_fa);
|
||||
|
||||
rep_time = 500;
|
||||
|
||||
% rou >= c / 2B = 15(m)
|
||||
|
||||
%% 产生发射信号
|
||||
N = Tr * fs;
|
||||
N_high = Tp * fs;
|
||||
|
||||
signal_t = zeros(1, N * numP);
|
||||
for j = 1: numP
|
||||
for i = 1: N_high
|
||||
tp = ((j - 1) * N + i) * (1 / fs);
|
||||
signal_t(1, i + (j-1)*N) = exp(1j * 2 * pi * (fc * tp + 0.5 * K * tp .^ 2));
|
||||
end
|
||||
end
|
||||
|
||||
multiple_r = signal_t(1, 1: N) * signal_t(1, 1: N)';
|
||||
% figure(1);
|
||||
% subplot(311)
|
||||
% plot(tP,real(signal_t));
|
||||
% xlabel('时间/t');ylabel('幅度');
|
||||
% title('发射信号');
|
||||
|
||||
%% 设置目标 1
|
||||
distance1 = 18000;
|
||||
v1 = 187.5;
|
||||
|
||||
tau1 = distance1 * 2 / c;
|
||||
n_tau1 = round(tau1 * fs);
|
||||
|
||||
f_d1 = 2 * fc * v1 / c;
|
||||
|
||||
alpha1 = alpha_prop;
|
||||
% alpha1 = 0;
|
||||
|
||||
signal_r1 = zeros(1, N * numP);
|
||||
for j = 1: numP
|
||||
for i = 1: N
|
||||
temp = i - n_tau1;
|
||||
if temp >= 1 && temp <= N_high
|
||||
tp = ((j - 1) * N + i - n_tau1) * (1 / fs);
|
||||
signal_r1(1, i + (j-1)*N) = alpha1 * exp(1j*2*pi*((fc - f_d1) * tp + 0.5 * K * tp .^ 2));
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
%% 生成回波
|
||||
signal_r = signal_r1;
|
||||
% figure(1);
|
||||
% subplot(312)
|
||||
% plot(tP, real(signal_r));
|
||||
% xlabel('时间/t');
|
||||
% ylabel('幅度');
|
||||
% title('回波信号');
|
||||
|
||||
%%
|
||||
P_fa_MF_cnt = zeros(len_SNR, len_P_fa, rep_time);
|
||||
P_d_MF_cnt = zeros(len_SNR, len_P_fa, rep_time);
|
||||
|
||||
% 多普勒处理
|
||||
F_ori = dftmtx(slow_len);
|
||||
F = F_ori(1:numP,:);
|
||||
multiple_d = F(:,1)' * F(:,1);
|
||||
|
||||
parfor rep = 1: rep_time
|
||||
% for rep = 1: rep_time
|
||||
|
||||
for cnt_SNR = 1: len_SNR
|
||||
|
||||
%% 噪声处理
|
||||
noise = random('Normal', 0, sigma_n/sqrt(2), 1, N * numP) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, N * numP);
|
||||
alpha_SNR = sqrt(10^(SNR(cnt_SNR)/10) * sigma_n^2 / (multiple_r * multiple_d));
|
||||
signal_r_n = signal_r * alpha_SNR + noise;
|
||||
|
||||
% figure(1);
|
||||
% subplot(313)
|
||||
% plot(tP, real(signal_r_n));
|
||||
% xlabel('时间/t');
|
||||
% ylabel('幅度');
|
||||
% title('回波+噪声信号');
|
||||
|
||||
|
||||
%% 匹配滤波——按Tr划分
|
||||
Srange = zeros(N, numP);
|
||||
|
||||
% A = generate_matrix(transpose(signal_t(1, 1: N)), 1);
|
||||
|
||||
count_fa=0;
|
||||
|
||||
for i = 1: numP
|
||||
sr = signal_r_n(1, 1 + (i - 1) * N: i * N);
|
||||
st = signal_t(1, 1 + (i - 1) * N: i * N);
|
||||
|
||||
% A = generate_matrix(transpose(st), 1);
|
||||
% mf = A' * transpose(sr);
|
||||
% filterred_rf_r = mf ./ multiple_r;
|
||||
|
||||
% 匹配滤波
|
||||
rf_r_fft_conj = conj(fft(st));
|
||||
filterred_rf_r = ifft(fft(sr) .* rf_r_fft_conj);
|
||||
filterred_rf_r = filterred_rf_r ./ multiple_r;
|
||||
|
||||
Srange(:,i) = filterred_rf_r';
|
||||
|
||||
end
|
||||
|
||||
|
||||
% figure(101)
|
||||
% mesh(abs(Srange))
|
||||
% title('按Tr进行匹配滤波结果')
|
||||
|
||||
%% R匹配滤波——分布验证
|
||||
% Before MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2)
|
||||
% b~N(0,sigma_n^2 / 2)
|
||||
|
||||
% After MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2 / multiple_r)
|
||||
% b~N(0,sigma_n^2 / 2 / multiple_r)
|
||||
% a^2 + b^2 ~ chi^2(2) * sigma_n^2 / 2 / multiple_r
|
||||
|
||||
% stat_R = abs(Srange).^2;
|
||||
% P_fa = 0.1;
|
||||
% kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2 / multiple;
|
||||
% sum(sum(stat_R>kd))/numel(stat_R)
|
||||
|
||||
%% 多普勒滤波
|
||||
Srd = zeros(N, slow_len);
|
||||
for i = 1: N
|
||||
% for i = 1201: 1201
|
||||
% Srd(i, :) = fftshift(fft([Srange(i, :),zeros(1,slow_len-numP)]));
|
||||
% x_slow = [Srange(i, :),zeros(1,slow_len-numP)];
|
||||
% % y' = F * x'
|
||||
% y_slow2 = F * transpose(x_slow);
|
||||
y_slow2 = transpose(F) * transpose(Srange(i, :))./ multiple_d;
|
||||
% figure(444)
|
||||
% plot(abs(y_slow2))
|
||||
|
||||
Srd(i, :) = fftshift(transpose(y_slow2));
|
||||
end
|
||||
|
||||
%% D匹配滤波——分布验证
|
||||
% Before MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2 / multiple)
|
||||
% b~N(0,sigma_n^2 / 2 / multiple)
|
||||
|
||||
% After MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2 / multiple_r / multiple_d)
|
||||
% b~N(0,sigma_n^2 / 2 / multiple_r / multiple_d)
|
||||
% a^2 + b^2 ~ chi^2(2) * sigma_n^2 / 2 / multiple_r / multiple_d
|
||||
|
||||
% 检测
|
||||
target_node_d = round(v1 / (c / fc * PRF / 2 / slow_len)) + 1 + slow_len / 2;
|
||||
target_node_r = n_tau1 + 1;
|
||||
|
||||
stat_D = abs(Srd).^2;
|
||||
H0_index = zeros(size(Srd));
|
||||
H0_index(target_node_r-N_high:target_node_r+N_high, :) = 1;
|
||||
H0_index(target_node_r, target_node_d) = 0;
|
||||
H1_index = zeros(size(Srd));
|
||||
H1_index(target_node_r, target_node_d) = 1;
|
||||
|
||||
kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2 / multiple_r / multiple_d;
|
||||
|
||||
for cnt_h_th = 1: len_P_fa
|
||||
P_fa_MF_cnt(cnt_SNR, cnt_h_th, rep) = sum(sum(stat_D(H0_index > 0) > kd(cnt_h_th))) / sum(sum(H0_index));
|
||||
P_d_MF_cnt(cnt_SNR, cnt_h_th, rep) = sum(sum(stat_D(H1_index > 0) > kd(cnt_h_th))) / sum(sum(H1_index));
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
fprintf('%d\n', rep);
|
||||
|
||||
end
|
||||
|
||||
P_fa_MF = mean(P_fa_MF_cnt, 3);
|
||||
P_d_MF = mean(P_d_MF_cnt, 3);
|
||||
|
||||
%% plot
|
||||
figure(1);
|
||||
loglog(P_fa,P_fa_MF(1,:), 'linewidth', 2);
|
||||
% hold on;
|
||||
% grid on;
|
||||
% loglog(P_fa,P_fa_MF(2,:), 'linewidth', 2);
|
||||
% loglog(P_fa,P_fa_MF(3,:), 'linewidth', 2);
|
||||
% loglog(P_fa,P_fa_MF(4,:), 'linewidth', 2);
|
||||
% loglog(P_fa,P_fa_MF(5,:), 'linewidth', 2);
|
||||
% loglog(P_fa,P_fa_MF(6,:), 'linewidth', 2);
|
||||
% legend('SNR = 0','SNR = 2','SNR = 4','SNR = 6','SNR = 8','SNR = 10');
|
||||
xlabel('P_fa');
|
||||
ylabel('Actual P_fa');
|
||||
|
||||
figure(2);
|
||||
semilogx(P_fa,P_d_MF(1,:), 'linewidth', 2);
|
||||
% hold on;
|
||||
% grid on;
|
||||
% semilogx(P_fa,P_d_MF(2,:), 'linewidth', 2);
|
||||
% semilogx(P_fa,P_d_MF(3,:), 'linewidth', 2);
|
||||
% semilogx(P_fa,P_d_MF(4,:), 'linewidth', 2);
|
||||
% semilogx(P_fa,P_d_MF(5,:), 'linewidth', 2);
|
||||
% semilogx(P_fa,P_d_MF(6,:), 'linewidth', 2);
|
||||
% legend('SNR = 0','SNR = 2','SNR = 4','SNR = 6','SNR = 8','SNR = 10');
|
||||
xlabel('P_fa');
|
||||
ylabel('P_d');
|
||||
|
||||
distance_temp = 0: (c / fs / 2) : (Tr * c / 2 - c / fs / 2);
|
||||
speed_temp = (-PRF / 2: PRF / slow_len: PRF / 2 - PRF / slow_len) * (c / fc) / 2 ;
|
||||
|
||||
% save test_Rmf_Dmf.mat ...
|
||||
% distance_temp...
|
||||
% speed_temp...
|
||||
% Srd...
|
||||
% Srange...
|
||||
% tP...
|
||||
% signal_t...
|
||||
% signal_r...
|
||||
% signal_r_n...
|
||||
% F...
|
||||
% slow_len;
|
||||
|
||||
save test_MC_Pfa_Rmf_Dmf.mat ...
|
||||
SNR...
|
||||
P_fa...
|
||||
P_fa_MF...
|
||||
P_d_MF...
|
||||
sigma_n...
|
||||
B...
|
||||
Tp...
|
||||
fs...
|
||||
Ts...
|
||||
K...
|
||||
fc...
|
||||
PRF...
|
||||
Tr...
|
||||
numP...
|
||||
t...
|
||||
tP...
|
||||
slow_len;
|
||||
@@ -0,0 +1,380 @@
|
||||
clc
|
||||
clear
|
||||
% rng(1)
|
||||
|
||||
%% 参数设置
|
||||
B = 5e6; % 信号带宽5MHz
|
||||
Tp = 20e-6; % 脉宽100us
|
||||
fs = 2 * B; % 采样频率
|
||||
Ts = 1 / fs; % 采样周期
|
||||
K = B / Tp; % 线性调频率
|
||||
fc = 1.25e9; % 载波频率1.25GHz
|
||||
PRF = 5000; % 脉冲重复频率
|
||||
Tr = 1 / PRF; % 脉冲重复间隔
|
||||
|
||||
numP = 64; % 脉冲数量
|
||||
|
||||
t = 0: 1 / fs: Tr - 1 / fs;
|
||||
tP = 0: 1 / fs: Tr * numP - 1 / fs;
|
||||
|
||||
c = 3e8; % 光速
|
||||
|
||||
%%
|
||||
sigma_n = 0.1; % 噪声标准差
|
||||
alpha_prop = 0.8; % 目标散射点强度
|
||||
|
||||
SNR = 50; % 信噪比（积累后信噪比）
|
||||
|
||||
% rou >= c / 2B = 15(m)
|
||||
|
||||
slow_len = 128; % 多普勒维采样点
|
||||
|
||||
% 是否有目标，0无目标，1有目标
|
||||
has_target = 0;
|
||||
|
||||
%% 产生发射信号
|
||||
% t=-Tp/2:1/fs:Tp/2-1/fs;
|
||||
N = Tr * fs;
|
||||
N_high = Tp * fs;
|
||||
|
||||
signal_t = zeros(1, N * numP);
|
||||
for j = 1: numP
|
||||
for i = 1: N_high
|
||||
tp = ((j - 1) * N + i) * (1 / fs);
|
||||
signal_t(1, i + (j-1)*N) = exp(1j * 2 * pi * (fc * tp + 0.5 * K * tp .^ 2));
|
||||
end
|
||||
end
|
||||
|
||||
% signal_t = zeros(1, N);
|
||||
% for i = 1:N_high
|
||||
% tp = (i - 1) * (1 / fs);
|
||||
% signal_t(1, i) = exp(1j*2*pi*(fc*tp+0.5*K*tp.^2));
|
||||
% end
|
||||
|
||||
figure(1);
|
||||
subplot(311)
|
||||
plot(tP,real(signal_t));
|
||||
xlabel('时间/t');ylabel('幅度');
|
||||
title('发射信号');
|
||||
|
||||
%%
|
||||
% 距离匹配滤波增益
|
||||
multiple_r = signal_t(1, 1: N) * signal_t(1, 1: N)';
|
||||
|
||||
% 多普勒匹配滤波增益
|
||||
F_ori = dftmtx(slow_len);
|
||||
F = F_ori(1:numP,:);
|
||||
multiple_d = F(:,1)' * F(:,1);
|
||||
|
||||
%% 设置目标 1
|
||||
distance1 = 18000;
|
||||
v1 = 187.5;
|
||||
|
||||
tau1 = distance1 * 2 / c; % 时延
|
||||
n_tau1 = round(tau1 * fs); % 对应采样点
|
||||
|
||||
f_d1 = 2 * fc * v1 / c; % 多普勒频率
|
||||
|
||||
% 根据 SNR 设置回波散射强度
|
||||
if has_target
|
||||
alpha1 = alpha_prop * sqrt(10^(SNR/10) * sigma_n^2 / (multiple_r * multiple_d));
|
||||
else
|
||||
alpha1 = 0;
|
||||
end
|
||||
|
||||
% 生成目标 1 回波
|
||||
signal_r1 = zeros(1, N * numP);
|
||||
for j = 1: numP
|
||||
for i = 1: N
|
||||
temp = i - n_tau1;
|
||||
if temp >= 1 && temp <= N_high
|
||||
% tp = (temp - 1) * (1 / fs);
|
||||
tp = ((j - 1) * N + i - n_tau1) * (1 / fs);
|
||||
signal_r1(1, i + (j-1)*N) = alpha1 * exp(1j*2*pi*((fc - f_d1) * tp + 0.5 * K * tp .^ 2));
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
% signal_r1 = zeros(1, N);
|
||||
% for i = 1:N
|
||||
% temp = i - n_tau1;
|
||||
% if temp >= 1 && temp <= N_high
|
||||
% signal_r1(1, i) = alpha1 * signal_t(1, temp);
|
||||
% end
|
||||
% end
|
||||
|
||||
%% 生成回波
|
||||
signal_r = signal_r1;
|
||||
figure(1);
|
||||
subplot(312)
|
||||
plot(tP, real(signal_r));
|
||||
% plot(real(signal_r));
|
||||
xlabel('时间/t');
|
||||
ylabel('幅度');
|
||||
title('回波信号');
|
||||
|
||||
%% 噪声处理
|
||||
noise = random('Normal', 0, sigma_n/sqrt(2), 1, N * numP) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, N * numP);
|
||||
signal_r_n = signal_r + noise;
|
||||
|
||||
figure(1);
|
||||
subplot(313)
|
||||
plot(tP, real(signal_r_n));
|
||||
xlabel('时间/t');
|
||||
ylabel('幅度');
|
||||
title('回波+噪声信号');
|
||||
|
||||
|
||||
%% 对回波求多普勒频移
|
||||
% s_R_slow = zeros(1, slow_len);
|
||||
%
|
||||
% slow_freqs = ((0:slow_len - 1) .* (1 / Tr)) / slow_len;
|
||||
%
|
||||
% init_idx = 1;
|
||||
% while abs(signal_r_n(init_idx)) == 0
|
||||
% init_idx = init_idx + 1;
|
||||
% end
|
||||
%
|
||||
% for i = 0:numP - 1
|
||||
% s_R_slow(i + 1) = signal_r_n(init_idx + i * N + 1);
|
||||
% end
|
||||
%
|
||||
% s_R_fft = fftshift(fft(s_R_slow));
|
||||
|
||||
% figure(3)
|
||||
% plot(real(s_R_fft))
|
||||
|
||||
% [~, max_index_s_R] = max(s_R_fft);
|
||||
% freq_s_R = slow_freqs(max_index_s_R);
|
||||
%
|
||||
% f_d = 1 / Tr - freq_s_R;
|
||||
% v = f_d * c / (2 * (fc - B / 2))
|
||||
|
||||
|
||||
%% 距离匹配滤波——按Tr划分
|
||||
Srange = zeros(N, numP);
|
||||
|
||||
for i = 1: numP
|
||||
sr = signal_r_n(1, 1 + (i - 1) * N: i * N);
|
||||
st = signal_t(1, 1 + (i - 1) * N: i * N);
|
||||
|
||||
% % 匹配滤波
|
||||
% A = generate_matrix(transpose(st), 1);
|
||||
% mf = A' * transpose(sr);
|
||||
% filterred_rf_r = mf ./ multiple;
|
||||
|
||||
% 匹配滤波
|
||||
rf_r_fft_conj = conj(fft(st));
|
||||
filterred_rf_r = ifft(fft(sr) .* rf_r_fft_conj);
|
||||
filterred_rf_r = filterred_rf_r ./ multiple_r;
|
||||
|
||||
% figure(2000)
|
||||
% plot(abs(filterred_rf_r))
|
||||
|
||||
Srange(:,i) = filterred_rf_r';
|
||||
|
||||
end
|
||||
|
||||
|
||||
figure(101)
|
||||
mesh(abs(Srange))
|
||||
title('按Tr进行匹配滤波结果')
|
||||
|
||||
%% R匹配滤波——分布验证
|
||||
% Before MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2)
|
||||
% b~N(0,sigma_n^2 / 2)
|
||||
|
||||
% After MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2 / multiple_r)
|
||||
% b~N(0,sigma_n^2 / 2 / multiple_r)
|
||||
% a^2 + b^2 ~ chi^2(2) * sigma_n^2 / 2 / multiple_r
|
||||
|
||||
% stat_R = abs(Srange).^2;
|
||||
% P_fa = 0.1;
|
||||
% kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2 / multiple;
|
||||
% sum(sum(stat_R>kd))/numel(stat_R)
|
||||
|
||||
%% 多普勒滤波
|
||||
% F_ori = dftmtx(slow_len);
|
||||
% F = F_ori(1:numP,:);
|
||||
% multiple_d = F(:,1)' * F(:,1);
|
||||
|
||||
F_inv = conj(F)/slow_len;
|
||||
% ttttttt = F_inv * transpose(F);
|
||||
|
||||
% VAMP参数
|
||||
delta_VAMP = 1e-6;
|
||||
iter_max = 1000;
|
||||
lambda_val = 0.05;
|
||||
lambda = zeros(slow_len,1) + lambda_val;
|
||||
|
||||
% 放大到行正交
|
||||
A = (sqrt(slow_len) * eye(numP)) * F_inv;
|
||||
|
||||
n = slow_len;
|
||||
m = numP;
|
||||
gamma = numP / slow_len;
|
||||
|
||||
% 设定虚警率
|
||||
P_fa = [1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1];
|
||||
count_Pfa_d = zeros(size(length(P_fa), N));
|
||||
|
||||
Srd = zeros(N, slow_len);
|
||||
for i = 1: N
|
||||
% for i = 1201: 1201
|
||||
x_slow = Srange(i, :); % 慢时间采样
|
||||
y = (sqrt(slow_len) * eye(numP)) * transpose(x_slow); % 进行相应放大
|
||||
|
||||
% LASSO 求解
|
||||
% x_LASSO = FISTA(y, A, lambda_val, delta_VAMP);
|
||||
[x_LASSO,y_d] = cVAMPro(y,A,lambda,delta_VAMP,iter_max);
|
||||
% x_LASSO = x_LASSO ./ multiple_d ./ 2;
|
||||
|
||||
% CROD求去偏
|
||||
rho_active = sum(abs(x_LASSO) > 1e-3)/n;
|
||||
Q_hat = (gamma - rho_active)/(1 - rho_active);
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = 1;
|
||||
while(diff > 1e-4)
|
||||
Rho_pre = Rho;
|
||||
Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
|
||||
diff = abs(Rho - Rho_pre);
|
||||
end
|
||||
Q_hat = (gamma-Rho)/(1-Rho);
|
||||
x_d_CROD = x_LASSO + A'*(y - A*x_LASSO)/Q_hat;
|
||||
|
||||
% CROD求门限和检验统计量
|
||||
RSS = sum(abs(y - A * x_LASSO).^2)/m;
|
||||
chi = Rho*(1 - Rho)/(gamma - Rho);
|
||||
if chi ~= 0
|
||||
chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
|
||||
z = -(1 - chi + chi_temp) / (2*chi);
|
||||
z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
|
||||
G_prime = (z + 1/chi);
|
||||
G_wprime = (z_prime + 1/chi/chi);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
else
|
||||
G_prime = gamma;
|
||||
G_wprime = gamma*(1-gamma);
|
||||
chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
|
||||
+ (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
|
||||
end
|
||||
sigma_CROD = sqrt(2*chi_hat) / Q_hat;
|
||||
|
||||
% 【经验分布、待解决】
|
||||
sigma_CROD = sigma_CROD*0.895;
|
||||
|
||||
stat_CROD = abs(x_d_CROD / sigma_CROD).^2; % 统计量
|
||||
|
||||
h_thd = chi2inv(1 - P_fa, 2) / 2; % 门限
|
||||
|
||||
% 检测
|
||||
% 目标点位
|
||||
target_node_d = round(v1 / (c / fc * PRF / 2 / slow_len)) + 1;
|
||||
target_node_r = n_tau1 + 1;
|
||||
|
||||
for j = 1:length(P_fa)
|
||||
|
||||
if i ~= target_node_r % 非目标点的距离慢采样结果（H0假设）
|
||||
count_Pfa_d(j,i) = sum(stat_CROD > h_thd(j)) / slow_len;
|
||||
else % 要考虑目标点的距离慢采样结果
|
||||
stat_index = ones(size(stat_CROD));
|
||||
stat_index(target_node_d) = 0; % 去掉目标点（H0假设不考虑目标点）
|
||||
count_Pfa_d(j,i) = sum(stat_CROD(stat_index > 0) > h_thd(j)) / sum(stat_index);
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
|
||||
% if i ~= target_node_r
|
||||
% count_Pfa_d(i) = sum(stat_CROD > h_thd) / slow_len;
|
||||
% else
|
||||
% stat_index = ones(size(stat_CROD));
|
||||
% stat_index(target_node_d) = 0;
|
||||
% count_Pfa_d(i) = sum(stat_CROD(stat_index > 0) > h_thd) / sum(stat_index);
|
||||
% end
|
||||
|
||||
% count_Pfa_d(i)
|
||||
|
||||
% figure(444)
|
||||
% % 归一化
|
||||
% plot(abs(x_d_CROD ./ multiple_d ./ 2))
|
||||
% figure(555)
|
||||
% % 统计量
|
||||
% plot(stat_CROD)
|
||||
|
||||
Srd(i, :) = fftshift(transpose(x_d_CROD ./ multiple_d ./ 2));
|
||||
end
|
||||
|
||||
%% DCS——分布验证
|
||||
output = mean(count_Pfa_d,2); % 把所有多普勒维检测结果求平均得到总体虚警率
|
||||
|
||||
figure(4000)
|
||||
loglog(P_fa,output)
|
||||
xlabel('P_fa Set')
|
||||
ylabel('Actual P_fa')
|
||||
title('P_fa')
|
||||
|
||||
%%
|
||||
% Srd_bf_mf = zeros(N, slow_len);
|
||||
% for i = 1: N
|
||||
% Srd_bf_mf(i, :) = fftshift(fft([Srange_bf_mf(i, :),zeros(1,slow_len-numP)]));
|
||||
% end
|
||||
|
||||
% figure(102)
|
||||
% mesh(abs(Srd))
|
||||
|
||||
% figure(105)
|
||||
% mesh(abs(Srd_bf_mf))
|
||||
|
||||
% figure(1)
|
||||
% subplot(414)
|
||||
% t2 = 0: (c / fs / 2) : (Tr * numP * c / 2 - c / fs / 2);
|
||||
% plot(tP, abs(filterred_rf_r) ./ N_high)
|
||||
% % plot(abs(filterred_rf_r))
|
||||
% xlabel('距离')
|
||||
% title('匹配滤波')
|
||||
|
||||
% distance_temp = (0:N - 1) * fs * c / N / 2 / K;
|
||||
|
||||
%% RD 绘制
|
||||
Srd = Srd';
|
||||
|
||||
lambda = c / fc;
|
||||
% 40*lambda*PRF/2/128
|
||||
distance_temp = 0: (c / fs / 2) : (Tr * c / 2 - c / fs / 2);
|
||||
speed_temp = (-PRF / 2: PRF / slow_len: PRF / 2 - PRF / slow_len) * lambda / 2 ;
|
||||
% speed_temp = (-numP / 2: numP / slow_len: numP / 2 - numP / slow_len) * lambda / Tr / numP / 2 ;
|
||||
|
||||
figure(5)
|
||||
[X, Y] = meshgrid(distance_temp, speed_temp);
|
||||
mesh(X, Y, abs(Srd));
|
||||
xlabel('距离(m)');
|
||||
ylabel('速度(m/s)');
|
||||
zlabel('信号幅值');
|
||||
title('2维RD图');
|
||||
|
||||
figure(6)
|
||||
imagesc(distance_temp, speed_temp, abs(Srd));
|
||||
title('Range-Doppler Image');
|
||||
xlabel('Range (m)');
|
||||
ylabel('Speed (m/s)');
|
||||
set(gca, 'FontSize', 18);
|
||||
set(gcf, 'position', [200, 300, 800, 600]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
% save test_Rmf_Dmf.mat ...
|
||||
% distance_temp...
|
||||
% speed_temp...
|
||||
% Srd...
|
||||
% Srange...
|
||||
% tP...
|
||||
% signal_t...
|
||||
% signal_r...
|
||||
% signal_r_n...
|
||||
% F...
|
||||
% slow_len;
|
||||
@@ -0,0 +1,311 @@
|
||||
clc
|
||||
clear
|
||||
% rng(1)
|
||||
|
||||
%% 参数设置
|
||||
B = 5e6; % 信号带宽5MHz
|
||||
Tp = 20e-6; % 脉宽100us
|
||||
fs = 2 * B; % 采样频率
|
||||
Ts = 1 / fs; % 采样周期
|
||||
K = B / Tp; % 线性调频率
|
||||
fc = 1.25e9; % 载波频率1.25GHz
|
||||
PRF = 5000; % 脉冲重复频率
|
||||
Tr = 1 / PRF; % 脉冲重复间隔
|
||||
|
||||
numP = 64; % 脉冲数量
|
||||
|
||||
t = 0: 1 / fs: Tr - 1 / fs;
|
||||
tP = 0: 1 / fs: Tr * numP - 1 / fs;
|
||||
|
||||
c = 3e8; % 光速
|
||||
|
||||
%%
|
||||
sigma_n = 0.1; % 噪声标准差
|
||||
alpha_prop = 0.8; % 目标散射点强度
|
||||
|
||||
SNR = 50; % 信噪比（积累后信噪比）
|
||||
|
||||
% rou >= c / 2B = 15(m)
|
||||
|
||||
slow_len = 128; % 多普勒维采样点
|
||||
|
||||
% 是否有目标，0无目标，1有目标
|
||||
has_target = 1;
|
||||
|
||||
%% 产生发射信号
|
||||
% t=-Tp/2:1/fs:Tp/2-1/fs;
|
||||
N = Tr * fs;
|
||||
N_high = Tp * fs;
|
||||
|
||||
signal_t = zeros(1, N * numP);
|
||||
for j = 1: numP
|
||||
for i = 1: N_high
|
||||
tp = ((j - 1) * N + i) * (1 / fs);
|
||||
signal_t(1, i + (j-1)*N) = exp(1j * 2 * pi * (fc * tp + 0.5 * K * tp .^ 2));
|
||||
end
|
||||
end
|
||||
|
||||
% signal_t = zeros(1, N);
|
||||
% for i = 1:N_high
|
||||
% tp = (i - 1) * (1 / fs);
|
||||
% signal_t(1, i) = exp(1j*2*pi*(fc*tp+0.5*K*tp.^2));
|
||||
% end
|
||||
|
||||
figure(1);
|
||||
subplot(311)
|
||||
plot(tP,real(signal_t));
|
||||
xlabel('时间/t');ylabel('幅度');
|
||||
title('发射信号');
|
||||
|
||||
%%
|
||||
% 距离匹配滤波增益
|
||||
multiple_r = signal_t(1, 1: N) * signal_t(1, 1: N)';
|
||||
|
||||
% 多普勒匹配滤波增益
|
||||
F_ori = dftmtx(slow_len);
|
||||
F = F_ori(1:numP,:);
|
||||
multiple_d = F(:,1)' * F(:,1);
|
||||
|
||||
%% 设置目标 1
|
||||
distance1 = 18000;
|
||||
v1 = 187.5;
|
||||
|
||||
tau1 = distance1 * 2 / c; % 时延
|
||||
n_tau1 = round(tau1 * fs); % 对应采样点
|
||||
|
||||
f_d1 = 2 * fc * v1 / c; % 多普勒频率
|
||||
|
||||
% 根据 SNR 设置回波散射强度
|
||||
if has_target
|
||||
alpha1 = alpha_prop * sqrt(10^(SNR/10) * sigma_n^2 / (multiple_r * multiple_d));
|
||||
else
|
||||
alpha1 = 0;
|
||||
end
|
||||
|
||||
% 生成目标 1 回波
|
||||
signal_r1 = zeros(1, N * numP);
|
||||
for j = 1: numP
|
||||
for i = 1: N
|
||||
temp = i - n_tau1;
|
||||
if temp >= 1 && temp <= N_high
|
||||
% tp = (temp - 1) * (1 / fs);
|
||||
tp = ((j - 1) * N + i - n_tau1) * (1 / fs);
|
||||
signal_r1(1, i + (j-1)*N) = alpha1 * exp(1j*2*pi*((fc - f_d1) * tp + 0.5 * K * tp .^ 2));
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
% signal_r1 = zeros(1, N);
|
||||
% for i = 1:N
|
||||
% temp = i - n_tau1;
|
||||
% if temp >= 1 && temp <= N_high
|
||||
% signal_r1(1, i) = alpha1 * signal_t(1, temp);
|
||||
% end
|
||||
% end
|
||||
|
||||
%% 生成回波
|
||||
signal_r = signal_r1;
|
||||
figure(1);
|
||||
subplot(312)
|
||||
plot(tP, real(signal_r));
|
||||
% plot(real(signal_r));
|
||||
xlabel('时间/t');
|
||||
ylabel('幅度');
|
||||
title('回波信号');
|
||||
|
||||
|
||||
%% 噪声处理
|
||||
noise = random('Normal', 0, sigma_n/sqrt(2), 1, N * numP) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, N * numP);
|
||||
signal_r_n = signal_r + noise;
|
||||
|
||||
figure(1);
|
||||
subplot(313)
|
||||
plot(tP, real(signal_r_n));
|
||||
xlabel('时间/t');
|
||||
ylabel('幅度');
|
||||
title('回波+噪声信号');
|
||||
|
||||
|
||||
%% 对回波求多普勒频移
|
||||
% s_R_slow = zeros(1, slow_len);
|
||||
%
|
||||
% slow_freqs = ((0:slow_len - 1) .* (1 / Tr)) / slow_len;
|
||||
%
|
||||
% init_idx = 1;
|
||||
% while abs(signal_r_n(init_idx)) == 0
|
||||
% init_idx = init_idx + 1;
|
||||
% end
|
||||
%
|
||||
% for i = 0:numP - 1
|
||||
% s_R_slow(i + 1) = signal_r_n(init_idx + i * N + 1);
|
||||
% end
|
||||
%
|
||||
% s_R_fft = fftshift(fft(s_R_slow));
|
||||
|
||||
% figure(3)
|
||||
% plot(real(s_R_fft))
|
||||
|
||||
% [~, max_index_s_R] = max(s_R_fft);
|
||||
% freq_s_R = slow_freqs(max_index_s_R);
|
||||
%
|
||||
% f_d = 1 / Tr - freq_s_R;
|
||||
% v = f_d * c / (2 * (fc - B / 2))
|
||||
|
||||
|
||||
%% 距离匹配滤波——按Tr划分
|
||||
Srange = zeros(N, numP);
|
||||
|
||||
for i = 1: numP
|
||||
sr = signal_r_n(1, 1 + (i - 1) * N: i * N);
|
||||
st = signal_t(1, 1 + (i - 1) * N: i * N);
|
||||
|
||||
% % 匹配滤波
|
||||
% A = generate_matrix(transpose(st), 1);
|
||||
% mf = A' * transpose(sr);
|
||||
% filterred_rf_r = mf ./ multiple;
|
||||
|
||||
% 匹配滤波
|
||||
rf_r_fft_conj = conj(fft(st));
|
||||
filterred_rf_r = ifft(fft(sr) .* rf_r_fft_conj);
|
||||
filterred_rf_r = filterred_rf_r ./ multiple_r;
|
||||
|
||||
% figure(2000)
|
||||
% plot(abs(filterred_rf_r))
|
||||
|
||||
Srange(:,i) = filterred_rf_r';
|
||||
|
||||
end
|
||||
|
||||
figure(101)
|
||||
mesh(abs(Srange))
|
||||
title('按Tr进行匹配滤波结果')
|
||||
|
||||
%% R匹配滤波——分布验证
|
||||
% Before MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2)
|
||||
% b~N(0,sigma_n^2 / 2)
|
||||
|
||||
% After MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2 / multiple_r)
|
||||
% b~N(0,sigma_n^2 / 2 / multiple_r)
|
||||
% a^2 + b^2 ~ chi^2(2) * sigma_n^2 / 2 / multiple_r
|
||||
|
||||
% stat_R = abs(Srange).^2;
|
||||
% P_fa = 0.1;
|
||||
% kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2 / multiple;
|
||||
% sum(sum(stat_R>kd))/numel(stat_R)
|
||||
|
||||
%% 多普勒滤波
|
||||
Srd = zeros(N, slow_len);
|
||||
for i = 1: N
|
||||
% for i = 1201: 1201 % 目标点
|
||||
% Srd(i, :) = fftshift(fft([Srange(i, :),zeros(1,slow_len-numP)]));
|
||||
% x_slow = [Srange(i, :),zeros(1,slow_len-numP)];
|
||||
% % y' = F * x'
|
||||
% y_slow2 = F * transpose(x_slow);
|
||||
y_slow2 = transpose(F) * transpose(Srange(i, :))./ multiple_d;
|
||||
% figure(444)
|
||||
% plot(abs(y_slow2))
|
||||
|
||||
Srd(i, :) = fftshift(transpose(y_slow2));
|
||||
end
|
||||
|
||||
%% D匹配滤波——分布验证
|
||||
% Before MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2 / multiple)
|
||||
% b~N(0,sigma_n^2 / 2 / multiple)
|
||||
|
||||
% After MF
|
||||
% a+bi
|
||||
% a~N(0,sigma_n^2 / 2 / multiple_r / multiple_d)
|
||||
% b~N(0,sigma_n^2 / 2 / multiple_r / multiple_d)
|
||||
% a^2 + b^2 ~ chi^2(2) * sigma_n^2 / 2 / multiple_r / multiple_d
|
||||
|
||||
% 虚警率计算
|
||||
stat_D = abs(Srd).^2;
|
||||
P_fa = [1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1];
|
||||
kd = sigma_n^2 * chi2inv(1 - P_fa, 2) / 2 / multiple_r / multiple_d;
|
||||
count_Pfa = zeros(1,length(P_fa));
|
||||
for i = 1:length(P_fa)
|
||||
count_Pfa(i) = sum(sum(stat_D>kd(i)))/numel(stat_D);
|
||||
end
|
||||
figure(4000)
|
||||
loglog(P_fa,count_Pfa)
|
||||
xlabel('P_fa Set')
|
||||
ylabel('Actual P_fa')
|
||||
title('P_fa')
|
||||
%%
|
||||
% Srd_bf_mf = zeros(N, slow_len);
|
||||
% for i = 1: N
|
||||
% Srd_bf_mf(i, :) = fftshift(fft([Srange_bf_mf(i, :),zeros(1,slow_len-numP)]));
|
||||
% end
|
||||
|
||||
% figure(102)
|
||||
% mesh(abs(Srd))
|
||||
|
||||
% figure(105)
|
||||
% mesh(abs(Srd_bf_mf))
|
||||
|
||||
% figure(1)
|
||||
% subplot(414)
|
||||
% t2 = 0: (c / fs / 2) : (Tr * numP * c / 2 - c / fs / 2);
|
||||
% plot(tP, abs(filterred_rf_r) ./ N_high)
|
||||
% % plot(abs(filterred_rf_r))
|
||||
% xlabel('距离')
|
||||
% title('匹配滤波')
|
||||
|
||||
% distance_temp = (0:N - 1) * fs * c / N / 2 / K;
|
||||
|
||||
%% RD 绘制
|
||||
Srd = Srd';
|
||||
|
||||
lambda = c / fc;
|
||||
|
||||
distance_temp = 0: (c / fs / 2) : (Tr * c / 2 - c / fs / 2);
|
||||
speed_temp = (-PRF / 2: PRF / slow_len: PRF / 2 - PRF / slow_len) * lambda / 2 ;
|
||||
% speed_temp = (-numP / 2: numP / slow_len: numP / 2 - numP / slow_len) * lambda / Tr / numP / 2 ;
|
||||
|
||||
figure(5)
|
||||
[X, Y] = meshgrid(distance_temp, speed_temp);
|
||||
mesh(X, Y, abs(Srd));
|
||||
xlabel('距离(m)');
|
||||
ylabel('速度(m/s)');
|
||||
zlabel('信号幅值');
|
||||
title('2维RD图');
|
||||
|
||||
figure(6)
|
||||
imagesc(distance_temp, speed_temp, abs(Srd));
|
||||
title('Range-Doppler Image');
|
||||
xlabel('Range (m)');
|
||||
ylabel('Speed (m/s)');
|
||||
set(gca, 'FontSize', 18);
|
||||
set(gcf, 'position', [200, 300, 800, 600]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
figure(7)
|
||||
sence = zeros(1,N);
|
||||
sence(n_tau1) = 0.8;
|
||||
plot(distance_temp,sence,'Linewidth', 2);
|
||||
legend('Target')
|
||||
title('Target Scene');
|
||||
xlabel('Range (m)');
|
||||
ylabel('Target Scattering Intensity');
|
||||
ylim([0,1]);
|
||||
set(gca, 'FontSize', 18);
|
||||
set(gcf, 'position', [200, 300, 800, 600]);
|
||||
set(gca,'fontsize',20,'fontname','Times');
|
||||
|
||||
% save test_Rmf_Dmf.mat ...
|
||||
% distance_temp...
|
||||
% speed_temp...
|
||||
% Srd...
|
||||
% Srange...
|
||||
% tP...
|
||||
% signal_t...
|
||||
% signal_r...
|
||||
% signal_r_n...
|
||||
% F...
|
||||
% slow_len;
|
||||
@@ -0,0 +1,52 @@
|
||||
clc; clear;
|
||||
|
||||
phi = 0.5 * pi;
|
||||
N = 5000;
|
||||
t = linspace(-2 * pi, 2 * pi, N)';
|
||||
x = sin(2 * pi * t) + 3 * sin(3 * 2 * pi * t);
|
||||
x_2 = sin(2 * pi * t + phi) + 3 * sin(3 * 2 * pi * t + phi);
|
||||
|
||||
f = fftshift(fft(x));
|
||||
F = dftmtx(N);
|
||||
f_hat = fftshift(F * x);
|
||||
|
||||
f_2 = fftshift(fft(x_2));
|
||||
|
||||
J = circshift(eye(N), 1, 2);
|
||||
[J_eigenvector, J_eigenvalue] = eig(J);
|
||||
J_eigenvalue = diag(J_eigenvalue);
|
||||
|
||||
CMatrix = toeplitz(x);
|
||||
[CMatrix_eigenvector, CMatrix_eigenvalue] = eig(CMatrix);
|
||||
CMatrix_eigenvalue = diag(CMatrix_eigenvalue);
|
||||
|
||||
CMatrix_eigenvalue_hat = zeros(size(CMatrix_eigenvalue));
|
||||
for i = 1: N
|
||||
tmp = f(i) * J_eigenvalue(i) ^ (i - 1);
|
||||
CMatrix_eigenvalue_hat(i) = CMatrix_eigenvalue_hat(i) + tmp;
|
||||
end
|
||||
|
||||
|
||||
if 1 == 1
|
||||
subplot(311);
|
||||
plot(x);
|
||||
|
||||
subplot(312);
|
||||
plot(abs(f));
|
||||
|
||||
subplot(313);
|
||||
% plot(abs(f_hat));
|
||||
plot(abs(f_2));
|
||||
|
||||
% polarplot(angle(J_eigenvalue),abs(J_eigenvalue),"o")
|
||||
else
|
||||
subplot(311);
|
||||
plot(abs(f_hat));
|
||||
|
||||
subplot(312);
|
||||
plot(sort(abs(CMatrix_eigenvalue)));
|
||||
|
||||
subplot(313);
|
||||
plot(sort(abs(CMatrix_eigenvalue_hat)));
|
||||
|
||||
end
|
||||
Reference in New Issue
Block a user