diff --git a/.gitignore b/.gitignore index b7b90c0..51094f4 100644 --- a/.gitignore +++ b/.gitignore @@ -33,7 +33,6 @@ Manifest.toml *.asv *.pdf -0 - example *.mat *.fig *.tif diff --git a/0 - example/nsq - 恒虚警代码/cVAMPro.m b/0 - example/nsq - 恒虚警代码/cVAMPro.m new file mode 100644 index 0000000..ee21f8d --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/cVAMPro.m @@ -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 diff --git a/0 - example/nsq - 恒虚警代码/plot_Pd_Pfa_vamp_cal_stat.m b/0 - example/nsq - 恒虚警代码/plot_Pd_Pfa_vamp_cal_stat.m new file mode 100644 index 0000000..ffb6471 --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/plot_Pd_Pfa_vamp_cal_stat.m @@ -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'); \ No newline at end of file diff --git a/0 - example/nsq - 恒虚警代码/stat_window.m b/0 - example/nsq - 恒虚警代码/stat_window.m new file mode 100644 index 0000000..356cde7 --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/stat_window.m @@ -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 \ No newline at end of file diff --git a/0 - example/nsq - 恒虚警代码/test_Pd_Pfa_vamp_cal_stat.m b/0 - example/nsq - 恒虚警代码/test_Pd_Pfa_vamp_cal_stat.m new file mode 100644 index 0000000..2d35ff1 --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/test_Pd_Pfa_vamp_cal_stat.m @@ -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; + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/FISTA.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/FISTA.m new file mode 100644 index 0000000..24f7223 --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/FISTA.m @@ -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 \ No newline at end of file diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m new file mode 100644 index 0000000..dc7b77f --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m @@ -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]); + + + + + + + + + + diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m new file mode 100644 index 0000000..e018e4d --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 2/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_SNR.m @@ -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; + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/FISTA.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/FISTA.m new file mode 100644 index 0000000..24f7223 --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/FISTA.m @@ -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 \ No newline at end of file diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m new file mode 100644 index 0000000..40ef17b --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m @@ -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]); + + + + + + + + + + diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m new file mode 100644 index 0000000..b436281 --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 3/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_p0.m @@ -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; + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/FISTA.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/FISTA.m new file mode 100644 index 0000000..24f7223 --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/FISTA.m @@ -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 \ No newline at end of file diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m new file mode 100644 index 0000000..baf4bff --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/plot_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m @@ -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]); + + + + + + + + + + diff --git a/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m new file mode 100644 index 0000000..212a723 --- /dev/null +++ b/0 - example/nsq - 恒虚警代码/那/20221025_那斯琦_ICASSP2022/Fig. 4/test_Pfa_Pd_SNR_PF_CROD_CAMP_SDL_ROD_gamma.m @@ -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; + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/0 - example/wzk - RD图相关/cVAMPro.m b/0 - example/wzk - RD图相关/cVAMPro.m new file mode 100644 index 0000000..ee21f8d --- /dev/null +++ b/0 - example/wzk - RD图相关/cVAMPro.m @@ -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 diff --git a/0 - example/wzk - RD图相关/plot_MC_Pfa_Rmf_Dcs.m b/0 - example/wzk - RD图相关/plot_MC_Pfa_Rmf_Dcs.m new file mode 100644 index 0000000..68733b5 --- /dev/null +++ b/0 - example/wzk - RD图相关/plot_MC_Pfa_Rmf_Dcs.m @@ -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'); \ No newline at end of file diff --git a/0 - example/wzk - RD图相关/plot_MC_Pfa_Rmf_Dmf.m b/0 - example/wzk - RD图相关/plot_MC_Pfa_Rmf_Dmf.m new file mode 100644 index 0000000..6a32cb9 --- /dev/null +++ b/0 - example/wzk - RD图相关/plot_MC_Pfa_Rmf_Dmf.m @@ -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'); \ No newline at end of file diff --git a/0 - example/wzk - RD图相关/plot_Pd_Pfa_diff.m b/0 - example/wzk - RD图相关/plot_Pd_Pfa_diff.m new file mode 100644 index 0000000..93b10e3 --- /dev/null +++ b/0 - example/wzk - RD图相关/plot_Pd_Pfa_diff.m @@ -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'); + diff --git a/0 - example/wzk - RD图相关/readme.txt b/0 - example/wzk - RD图相关/readme.txt new file mode 100644 index 0000000..e19ce53 --- /dev/null +++ b/0 - example/wzk - RD图相关/readme.txt @@ -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 对比图绘制 \ No newline at end of file diff --git a/0 - example/wzk - RD图相关/test_MC_Pfa_Rmf_Dcs.m b/0 - example/wzk - RD图相关/test_MC_Pfa_Rmf_Dcs.m new file mode 100644 index 0000000..1a8bcc4 --- /dev/null +++ b/0 - example/wzk - RD图相关/test_MC_Pfa_Rmf_Dcs.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('按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); + + 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; \ No newline at end of file diff --git a/0 - example/wzk - RD图相关/test_MC_Pfa_Rmf_Dmf.m b/0 - example/wzk - RD图相关/test_MC_Pfa_Rmf_Dmf.m new file mode 100644 index 0000000..302653e --- /dev/null +++ b/0 - example/wzk - RD图相关/test_MC_Pfa_Rmf_Dmf.m @@ -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; \ No newline at end of file diff --git a/0 - example/wzk - RD图相关/test_Pfa_noT_Rmf_Dcs.m b/0 - example/wzk - RD图相关/test_Pfa_noT_Rmf_Dcs.m new file mode 100644 index 0000000..9800d74 --- /dev/null +++ b/0 - example/wzk - RD图相关/test_Pfa_noT_Rmf_Dcs.m @@ -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; \ No newline at end of file diff --git a/0 - example/wzk - RD图相关/test_Pfa_noT_Rmf_Dmf.m b/0 - example/wzk - RD图相关/test_Pfa_noT_Rmf_Dmf.m new file mode 100644 index 0000000..9279965 --- /dev/null +++ b/0 - example/wzk - RD图相关/test_Pfa_noT_Rmf_Dmf.m @@ -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; \ No newline at end of file diff --git a/Test.m b/Test.m new file mode 100644 index 0000000..a9637cc --- /dev/null +++ b/Test.m @@ -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 \ No newline at end of file