Add lyh-非满秩相变

.gitignore

@@ -37,3 +37,5 @@ Manifest.toml
 *.fig
 *.tif
 *.bmp
+*.jpg
+*.jpeg

0 - example/lyh - 相变/捷变频相变/FAR_block_sparsity_phase1.m

@@ -0,0 +1,48 @@
+close all;
+clear all;
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 20;
+epi = 0.02;
+result = zeros(N,25);
+for col = 4:N
+    for block_sparsity = 10:18
+        success_count = 0;
+        for loop = 1:trial
+            FAR_model = zeros(N,M*N);
+            %Cn = randperm(M)-1
+            for n = 0 : N-1
+                Cn = floor(rand()*M);
+                for q = 0 : N-1
+                    for p = 0:M-1
+                        FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                    end
+                end
+            end
+            col_choose = randperm(N,col);
+            FAR_model = FAR_model(col_choose,:);
+            sparse_signal = zeros(M,N);
+            block = randperm(N,block_sparsity);
+            sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+            y = FAR_model * sparse_signal(:);
+            cvx_begin
+            variable x(M,N) complex
+            norm21 = 0;
+            for i = 1:N
+                norm21 = norm21 + norm(x(:,i));
+            end
+            minimize(norm21)
+            subject to
+            FAR_model * x(:) == y
+            cvx_end
+            if norm(x(:)-sparse_signal(:))<tol
+                success_count = success_count+1;
+            end
+        end
+        result(col,block_sparsity) = success_count/trial;
+    end
+end
+save('FARblockepsilon2.mat');

0 - example/lyh - 相变/捷变频相变/FAR_sparsity_phase.m

@@ -0,0 +1,47 @@
+close all;
+clear all;
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 50;
+epi = 0.02;
+result = zeros(N,25);
+for col = 4:4:128
+    for block_sparsity = 1:25
+        success_count = 0;
+        for loop = 1:trial
+            FAR_model = zeros(N,M*N);
+            %Cn = randperm(M)-1
+            for n = 0 : N-1
+                Cn = floor(rand()*M);
+                for q = 0 : N-1
+                    for p = 0:M-1
+                        FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                    end
+                end
+            end
+            
+            col_choose = randperm(N,col);
+            FAR_model = FAR_model(col_choose,:);
+            sparse_signal = zeros(M,N);
+
+            block = randperm(N,block_sparsity);
+            sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+            
+            y = FAR_model * sparse_signal(:);
+            cvx_begin
+            variable x(M*N) complex
+            minimize(norm(x,1))
+            subject to
+            FAR_model * x == y
+            cvx_end
+            if norm(x-sparse_signal(:))<tol
+                success_count = success_count+1;
+            end
+        end
+        result(col,block_sparsity) = success_count/trial;
+    end
+end
+save('FARepsilon.mat');

0 - example/lyh - 相变/捷变频相变/theoretic.m

@@ -4,5 +4,4 @@ function n = theoretic(m,s,d)
     f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
     g = diff(f,t);
     t1 = solve(g);
-    n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
+    n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
-end

0 - example/lyh - 相变/高斯相变/gauss_phase_real.m

@@ -0,0 +1,27 @@
+N = 100;
+gauss_phase_res = zeros(100,100);
+for col = 1:100
+    %¾ØÕóÉú³É
+    for p =1:100
+        suc = 0;
+        for loop = 1:50
+            x1 = zeros(N,1);
+            q = randperm(N,p);
+            x1(q) = randn(p,1);
+            fai = randn(col,N);
+            b = fai*x1;
+            cvx_begin quiet
+                variable x(N)
+                minimize( norm( x, 1 ) )
+                subject to
+                    fai * x == b
+            cvx_end
+            %disp((norm(x-x1,1)))
+            if (norm(x-x1,1))<10e-5
+                 suc = suc+1;
+             end
+        end
+        gauss_phase_res(col,p)=suc/50;
+    end
+end
+save gauss_phase_real;

0 - example/lyh - 非满秩相变/FAR_block_sparsity_phase1.m

@@ -0,0 +1,229 @@
+close all;
+clear all;
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 50;
+epi = 0.02;
+result = zeros(N,25);
+for col = 4:4:128
+    for block_sparsity = 10:18
+        success_count = 0;
+        for loop = 1:trial
+            FAR_model = zeros(N,M*N);
+            %Cn = randperm(M)-1
+            for n = 0 : N-1
+                Cn = floor(rand()*M);
+                for q = 0 : N-1
+                    for p = 0:M-1
+                        FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                    end
+                end
+            end
+            col_choose = randperm(N,col);
+            FAR_model = FAR_model(col_choose,:);
+            sparse_signal = zeros(M,N);
+            block = randperm(N,block_sparsity);
+            sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+            y = FAR_model * sparse_signal(:);
+            cvx_begin
+            variable x(M,N) complex
+            norm21 = 0;
+            for i = 1:N
+                norm21 = norm21 + norm(x(:,i));
+            end
+            minimize(norm21)
+            subject to
+            FAR_model * x(:) == y
+            cvx_end
+            if norm(x(:)-sparse_signal(:))<tol
+                success_count = success_count+1;
+            end
+        end
+        result(col,block_sparsity) = success_count/trial;
+    end
+end
+%save('FARblockepsilon2.mat');
+
+%%
+close all;
+clear all;
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 30;
+epi = 0.02;
+result = zeros(1,25);
+col = 128;
+Prob = [1/3,1/3,1/6,1/6];
+for block_sparsity = 4:18
+    success_count = 0;
+    for loop = 1:trial
+        FAR_model = zeros(N,M*N);
+        %Cn = randperm(M)-1
+        for n = 0 : N-1
+            Cn = randsrc(1,1,[[0,1,2,3];Prob]);
+            for q = 0 : N-1
+                for p = 0:M-1
+                    FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                end
+            end
+        end
+        %col_choose = randperm(N,col);
+        %FAR_model = FAR_model(col_choose,:);
+        sparse_signal = zeros(M,N);
+        block = randperm(N,block_sparsity);
+        sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+        y = FAR_model * sparse_signal(:);
+        cvx_begin
+        variable x(M,N) complex
+        norm21 = 0;
+        for i = 1:N
+            norm21 = norm21 + norm(x(:,i));
+        end
+        minimize(norm21)
+        subject to
+        FAR_model * x(:) == y
+        cvx_end
+        if norm(x(:)-sparse_signal(:))<tol
+            success_count = success_count+1;
+        end
+    end
+    result(block_sparsity) = success_count/trial;
+end
+%%
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 30;
+epi = 0.02;
+result2 = zeros(2,25);
+col = 128;
+Prob = [1/3,1/3,1/6,1/6];
+for block_sparsity = 4:18
+    success_count = 0;
+    success_count2 = 0;
+    for loop = 1:trial
+        FAR_model = zeros(N,M*N);
+        %Cn = randperm(M)-1
+        for n = 0 : N-1
+            Cn = randsrc(1,1,[[0,1,2,3];Prob]);
+            for q = 0 : N-1
+                for p = 0:M-1
+                    FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                end
+            end
+        end
+        %col_choose = randperm(N,col);
+        %FAR_model = FAR_model(col_choose,:);
+        sparse_signal = zeros(M,N);
+        block = randperm(N,block_sparsity);
+        sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+        y = FAR_model * sparse_signal(:);
+        cvx_begin
+        variable x(M,N) complex
+        norm21 = 0;
+        for i = 1:N
+            norm21 = norm21 + norm(FAR_model(:,M*i-3:M*i)*x(:,i));
+        end
+        minimize(norm21)
+        subject to
+        FAR_model * x(:) == y
+        cvx_end
+        err1 = 0;
+        for i =1:N
+            err1 = err1 + norm(FAR_model(:,M*i-3:M*i)*(x(:,i)-sparse_signal(:,i)));
+        end
+        if err1<tol
+            success_count2 = success_count2 +1;
+        end
+        if norm(x(:)-sparse_signal(:))<tol
+            success_count = success_count+1;
+        end
+    end
+    result2(1,block_sparsity) = success_count/trial;
+    result2(2,block_sparsity) = success_count2/trial;
+end
+
+
+%%
+
+
+close all;
+clear all;
+clc;
+M = 4;
+N = 64;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 50;
+epi = 0.02;
+%result = zeros(N,25);
+%for col = 4:4:128
+col = 64;
+block_sparsity = 6;
+err = 0.1;
+sigma = 0.5*sqrt(block_sparsity);
+
+        
+        
+FAR_model = zeros(N,M*N);
+%Cn = randperm(M)-1
+for n = 0 : N-1
+    Cn = floor(rand()*M);
+    for q = 0 : N-1
+        for p = 0:M-1
+            FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn-1i*2*pi*q/N*n*(1+Cn*epi));
+        end
+    end
+end
+%col_choose = randperm(N,col);
+%FAR_model = FAR_model(col_choose,:);
+sparse_signal = zeros(M,N);
+block = randperm(N,block_sparsity);
+sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+y = FAR_model * sparse_signal(:) + sigma*(randn(N,1)) + 1i*sigma*(randn(N,1));
+cvx_begin
+variable x(M,N) complex
+norm21 = 0;
+for i = 1:N
+    norm21 = norm21 + norm(x(:,i));
+end
+minimize(norm21)
+subject to
+norm(FAR_model * x(:) - y)<err
+cvx_end
+
+
+x_block_norm = zeros(N,1);
+x_norm = zeros(N,1);
+for i = [1:N]
+    x_block_norm(i) = norm(FAR_model(:,M*i-3:M*i)*x(:,i)); 
+    x_norm(i) = norm(FAR_model(:,M*i-3:M*i)*sparse_signal(:,i));
+end
+
+hold on 
+plot(x_block_norm,'--o');
+plot(x_norm,'-.s');
+lgh = legend("Estimated","Ground Truth", ...
+    "MF");
+set(lgh,'interpreter','latex','FontName','Times New Roman')
+%set(gcf,'interpreter','latex','FontName','Times New Roman')
+xlabel("\fontname{Times New Roman}Velocity Cell Index");
+ylabel("\fontname{Times New Roman}Test Statistics \it{T_i}");
+xlim([0 75]);
+
+
+
+ 
+%end
+
+
+
+

0 - example/lyh - 非满秩相变/FAR_noise.m

@@ -0,0 +1,80 @@
+N = 32;
+M =4;
+R =2;
+epi = 0.02;
+noise = [0.1,0.1^(0.5)];
+%contri1 = zeros(M,N);
+%contri2 = zeros(M,N);
+con1 = zeros(20,2);
+con2 = zeros(20,2);
+re1 = zeros(20,2);
+re2 = zeros(20,2);
+
+for noise_k = 1:2
+for k = 1:20
+    for loop = 1:50
+
+    
+        FAR_model = zeros(N,M*N);
+        order = randperm(M)-1;
+        for n = 0 : N-1
+            Cn = order(ceil(rand()*R));
+            for q = 0 : N-1
+                for p = 0:M-1
+                    FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                end
+            end
+        end
+        
+        x = zeros(M,N);
+        col = randperm(N,k);
+        x(:,col) = randn(M,k) + 1i*randn(M,k);
+        y = FAR_model * x(:) + (randn(N,1)+1i*randn(N,1))*noise(noise_k);
+        
+        cvx_begin
+        variable x_e(M,N) complex
+        norm21 = 0;
+        for i = 1:N
+            norm21 = norm21 + norm(FAR_model(:,(i-1)*M+1:i*M)*x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        norm(FAR_model*x_e(:) - y) <= sqrt(N)*noise(noise_k);
+        cvx_end
+        contri2 = zeros(N,1);
+        re_err = zeros(N,1);
+        for i = 1:N
+            if ismember(i,col)
+                re_err = re_err+FAR_model(:,(i-1)*M+1:i*M)*(x_e(:,i));
+            end
+            contri2(i) = norm(FAR_model(:,(i-1)*M+1:i*M)*(x_e(:,i)));
+       	end
+		re1(k,noise_k) = re1(k,noise_k) +norm(FAR_model*x(:)-re_err)/norm(FAR_model*x(:));
+        con1(k,noise_k) = con1(k,noise_k) + 1-sum(contri2(col))/sum(contri2);      
+
+
+        cvx_begin
+        variable x_e(M,N) complex
+        norm21 = 0;
+        for i = 1:N
+            norm21 = norm21 + norm(x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        norm(FAR_model*x_e(:) - y) <= sqrt(N)*noise(noise_k);
+        cvx_end
+        re_err = zeros(N,1);
+        contri2 = zeros(N,1);
+        for i = 1:N
+            if ismember(i,col)
+                re_err = re_err+FAR_model(:,(i-1)*M+1:i*M)*(x_e(:,i));
+            end
+            contri2(i) = norm(FAR_model(:,(i-1)*M+1:i*M)*(x_e(:,i)));
+       	end
+		re2(k,noise_k) = re2(k,noise_k) +norm(FAR_model*x(:)-re_err)/norm(FAR_model*x(:));
+        con2(k,noise_k) = con2(k,noise_k) + 1-sum(contri2(col))/sum(contri2);   
+    end
+end
+end
+
+save("FAR_noise_0206.mat")

0 - example/lyh - 非满秩相变/gauss_n32m2d2base_class.m

@@ -0,0 +1,53 @@
+%初值设定
+N = 64; %最大行数
+D = 32; %块数
+m = 4;  %块的列数
+d = 2;  %块内rank/维度
+result = zeros(N,D);    %保留结果
+%按行数循环
+for n = 1:N
+A = zeros(n,m*D);
+B = zeros(n,d*D);
+%按稀疏度循环
+for k = 1:D
+    %50次重复实验
+    for j = 1:50
+	re_err = 0;
+        for i = 0:D-1
+        %生成高斯基集
+            set = randn(n,d);
+            
+        %生成随机的模为1的向量，与基集相乘得到块不满秩的高斯矩阵
+	    theta = rand(1,4)*2*pi;
+	    tmp = [sin(theta);cos(theta)];
+            A(:,m*i+1:m*i+m) = set*tmp;
+	   
+	    B(:,d*i+1:d*i+d) = set;
+        end
+        x = zeros(m,D);
+        col = randperm(D,k);
+        x(:,col) = randn(m,k);
+        y = A * x(:);
+        
+        %凸优化利用l21范数求解基集下的恢复问题
+        cvx_begin 
+        variable x_e(d,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        B*x_e(:) == y
+        cvx_end
+        
+        for i = 1:D
+            re_err = re_err+norm(B(:,(i-1)*d+1:i*d)*x_e(:,i)-A(:,(i-1)*m+1:i*m)*x(:,i));
+       	end
+	if re_err<10e-4
+		result(n,k) = result(n,k)+1;
+	end 
+    end
+end
+end
+save('gaussD32m2d2base.mat');

0 - example/lyh - 非满秩相变/gauss_n32m2d2rand.m

@@ -0,0 +1,55 @@
+%初值设定
+N = 64; %最大行数
+D = 32; %块数
+m = 4;  %块的列数
+d = 2;  %块内rank/维度
+result = zeros(N,D);    %保留结果
+
+%按行数循环
+for n = 1:N             
+A = zeros(n,m*D);       
+B = zeros(n,d*D);
+%按稀疏度循环
+for k = 1:D
+    %50次重复实验
+    for j = 1:50
+	re_err = 0;
+        for i = 0:D-1
+            %生成高斯基集
+            set = randn(n,d);
+        %生成随机的模为1的向量，与基集相乘得到块不满秩的高斯矩阵
+	    theta = rand(1,4)*2*pi;
+	    tmp = [sin(theta);cos(theta)];
+            A(:,m*i+1:m*i+m) = set*tmp;
+        %随机挑选块列向量作为降维后矩阵
+	    q = randperm(m,d);
+	    B(:,d*i+1:d*i+d) = A(:,m*i+q);
+        end
+        %生成稀疏块信号
+        x = zeros(m,D);
+        col = randperm(D,k);
+        x(:,col) = randn(m,k);
+        y = A * x(:);
+        %凸优化利用特殊l21范数求解降维后的恢复问题
+        cvx_begin 
+        variable x_e(d,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(B(:,(i-1)*d+1:i*d)*x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        B*x_e(:) == y
+        cvx_end
+        
+        %得到重建效果
+        for i = 1:D
+            re_err = re_err+norm(B(:,(i-1)*d+1:i*d)*x_e(:,i)-A(:,(i-1)*m+1:i*m)*x(:,i));
+       	end
+	if re_err<10e-4
+		result(n,k) = result(n,k)+1;
+	end 
+    end
+end
+end
+save('gaussD32m2d2randalter.mat');

0 - example/lyh - 非满秩相变/gauss_n32m2d2rand_class.m

@@ -0,0 +1,53 @@
+%初值设定
+N = 64; %最大行数
+D = 32; %块数
+m = 4;  %块的列数
+d = 2;  %块内rank/维度
+result = zeros(N,D);    %保留结果
+%按行数循环
+for n = 1:N
+A = zeros(n,m*D);
+B = zeros(n,d*D);
+%按稀疏度循环
+for k = 1:D
+    %50次重复实验
+    for j = 1:50
+	re_err = 0;
+        for i = 0:D-1
+        %生成高斯基集
+            set = randn(n,d);
+        %生成随机的模为1的向量，与基集相乘得到块不满秩的高斯矩阵
+	    theta = rand(1,4)*2*pi;
+	    tmp = [sin(theta);cos(theta)];
+            A(:,m*i+1:m*i+m) = set*tmp;
+	    q = randperm(m,d);
+	    B(:,d*i+1:d*i+d) = A(:,m*i+q);
+        end
+        %随机挑选块列向量作为降维后矩阵
+        x = zeros(m,D);
+        col = randperm(D,k);
+        x(:,col) = randn(m,k);
+        y = A * x(:);
+        %凸优化利用l21范数求解降维后的恢复问题
+        cvx_begin 
+        variable x_e(d,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        B*x_e(:) == y
+        cvx_end
+        
+         %得到重建效果
+        for i = 1:D
+            re_err = re_err+norm(B(:,(i-1)*d+1:i*d)*x_e(:,i)-A(:,(i-1)*m+1:i*m)*x(:,i));
+       	end
+	if re_err<10e-4
+		result(n,k) = result(n,k)+1;
+	end 
+    end
+end
+end
+save('gaussD32m2d2rand.mat');

0 - example/lyh - 非满秩相变/gauss_n32m4d2.m

@@ -0,0 +1,52 @@
+%初值设定
+N = 64; %最大行数
+D = 32; %块数
+m = 4;  %块的列数
+d = 2;  %块内rank/维度
+result = zeros(N,D);    %保留结果
+%按行数循环
+for n = 1:N
+A = zeros(n,m*D);
+B = zeros(n,d*D);
+%按稀疏度循环
+for k = 1:D
+    %50次重复实验
+    for j = 1:50
+	re_err = 0;
+        for i = 0:D-1
+        %生成高斯基集
+	    set = randn(n,d);
+        
+        
+        %生成随机的模为1的向量，与基集相乘得到块不满秩的高斯矩阵
+	    theta = rand(1,4)*2*pi;
+	    temp = [sin(theta);cos(theta)];
+            A(:,m*i+1:m*i+m) = set*temp;
+        end
+        x = zeros(m,D);
+        col = randperm(D,k);
+        x(:,col) = randn(m,k);
+        y = A * x(:);
+        
+        %凸优化利用特殊l21范数求解恢复问题
+        cvx_begin
+        variable x_e(m,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(A(:,(i-1)*m+1:i*m)*x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        A*x_e(:) == y
+        cvx_end
+        
+        for i = 1:D
+            re_err = re_err+norm(A(:,(i-1)*m+1:i*m)*(x_e(:,i)-x(:,i)));
+       	end
+	if re_err<10e-4
+		result(n,k) = result(n,k)+1;
+	end 
+    end
+end
+end
+save('gaussD32m4d2alter.mat');

0 - example/lyh - 非满秩相变/gauss_n32m4d2class.m

@@ -0,0 +1,51 @@
+%初值设定
+N = 64; %最大行数
+D = 32; %块数
+m = 4;  %块的列数
+d = 2;  %块内rank/维度
+result = zeros(N,D);    %保留结果
+%按行数循环
+for n = 1:N
+A = zeros(n,m*D);
+B = zeros(n,d*D);
+%按稀疏度循环
+for k = 1:D
+    %50次重复实验
+    for j = 1:50
+	re_err = 0;
+        for i = 0:D-1
+        %生成高斯基集
+	    set = randn(n,d);
+        %生成随机的模为1的向量，与基集相乘得到块不满秩的高斯矩阵
+	    theta = rand(1,4)*2*pi;
+	    temp = [sin(theta);cos(theta)];
+            A(:,m*i+1:m*i+m) = set*temp;
+        end
+        x = zeros(m,D);
+        col = randperm(D,k);
+        x(:,col) = randn(m,k);
+        y = A * x(:);
+        %凸优化利用l21范数求解恢复问题
+        cvx_begin
+        variable x_e(m,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        A*x_e(:) == y
+        cvx_end
+        
+        
+         %得到重建效果
+        for i = 1:D
+            re_err = re_err+norm(A(:,(i-1)*m+1:i*m)*(x_e(:,i)-x(:,i)));
+       	end
+	if re_err<10e-4
+		result(n,k) = result(n,k)+1;
+	end 
+    end
+end
+end
+save('gaussD32m4d2class.mat');

0 - example/lyh - 非满秩相变/gauss_noise.m

@@ -0,0 +1,24 @@
+N = 40; %行数
+D = 32; %块数
+m = 4;  %块的列数
+d = 2;  %块内rank/维度
+resulta = zeros(5000,1);
+resultb = zeros(5000,1);
+%resulta = zeros(5,20,50);    %保留结果
+%resultb = zeros(5,20,50);
+%按行数循环
+A = zeros(N,m*D);
+B = zeros(N,d*D);
+%按稀疏度循环
+n = N;
+parfor i = 0:199
+    t = floor(i/1000) + 1;
+    k = floor(mod(i,1000)/50)+1;
+    j = mod(i,50)+1;
+    [a,b] = test(t,k,j);
+    resulta(i+1) = a;
+    resultb(i+1) = b;
+end
+resulta = reshape(resulta,5,20,50);
+resultb = reshape(resultb,5,20,50);
+save('gaussnoise.mat');

0 - example/lyh - 非满秩相变/phaseplot.m

@@ -0,0 +1,48 @@
+load("result1.mat");
+%ambient_dim = 100;
+%sparsities = 1:100;
+d = 32;%¿é¸öÊý
+m = 2;%¿éÄÚÔªËØÊý
+l1recovery = zeros(100,1);
+for i = 1:32
+    l1recovery(i) = theoretic(m,i,d);
+end
+m = 1:32;
+n = 1:64;
+figure;
+contourf( m, n, result, 20, 'LineColor', 'none' );
+hold on
+cbar = colorbar;
+myplot2 = plot( n, l1recovery, 'Color', 'w', 'Linewidth', 3 );
+mylegend2 = legend(myplot2, '   Theory', 'Location', 'SouthEast' );
+mylegend2.Box = 'off';
+mylegend2.TextColor = [1,1,1];
+mylegend2.FontSize = 12;
+
+figure
+plot(con1(:,1)/400,'*-')
+xlabel("\fontname{Times New Roman} Block Sparsity \it K");
+ylabel("\fontname{Times New Roman} Block Contribution Error");
+savefig("contri_noise001.fig");
+saveas(gca,"contri_noise001.eps");
+
+figure
+plot(con1(:,2)/400,'*-')
+xlabel("\fontname{Times New Roman} Block Sparsity \it K");
+ylabel("\fontname{Times New Roman} Block Contribution Error");
+savefig("contri_noise01.fig");
+saveas(gca,"contri_noise01.eps");
+
+figure
+plot(re1(:,1)/400,'*-')
+xlabel("\fontname{Times New Roman} Block Sparsity \it K");
+ylabel("\fontname{Times New Roman} Reconstruction Error");
+savefig("re_noise001.fig");
+saveas(gca,"re_noise001.eps");
+
+figure
+plot(re1(:,2)/400,'*-')
+xlabel("\fontname{Times New Roman} Block Sparsity \it K");
+ylabel("\fontname{Times New Roman} Reconstruction Error");
+savefig("re_noise01.fig");
+saveas(gca,"re_noise01.eps");

0 - example/lyh - 非满秩相变/test.m

@@ -0,0 +1,228 @@
+%function[out1,out2,out3,out4,out5,out6] = test(t,k,~)
+
+%sigma = [0,0.01, 0.5, 1, 1.5, 2];
+n = 40;
+d = 2;
+D = 32;
+m = 4;
+%t = 1;
+k = 13;
+
+err = zeros(20,1);
+
+A = zeros(n,D*m);
+contri1 = zeros(n,D);
+contri2 = zeros(n,D);
+        for i = 0:32-1
+        %生成高斯基集
+	    set = randn(n,d);            
+        %生成随机的模为1的向量，与基集相乘得到块不满秩的高斯矩阵
+	    theta = rand(1,4)*2*pi;
+	    temp = [sin(theta);cos(theta)];
+            A(:,m*i+1:m*i+m) = set*temp*(i+1);
+        end
+        x = zeros(m,D);
+        col = randperm(D,k);
+        x(:,col) = randn(m,k);
+        y = A * x(:);
+        
+        weight =  diag(rand(D*m,1)*30);
+        %A2 = zeros(n,D*m);
+        %for i= 0:31
+        %    A2(:,m*i+1:m*i+m) = A(:,m*i+1:m*i+m) * weight(i+1);
+        %end
+        A2 = A*weight;
+        
+        %凸优化利用特殊l21范数求解恢复问题
+        cvx_begin
+        variable x_e(m,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(A(:,(i-1)*m+1:i*m)*x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        A*x_e(:) == y;
+        cvx_end
+        re_err = 0;
+        for i = 1:D
+            re_err = re_err+norm(A(:,(i-1)*m+1:i*m)*(x_e(:,i)-x(:,i)));
+            contri1(:,i) = A(:,(i-1)*m+1:i*m)*(x_e(:,i));
+       	end
+		out1 = re_err;
+        x3 = x_e;
+
+
+        
+        cvx_begin
+        variable x_e(m,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(A2(:,(i-1)*m+1:i*m)*x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        A2*x_e(:) == y;
+        cvx_end
+        re_err = 0;
+        for i = 1:D
+            re_err = re_err+norm(A2(:,(i-1)*m+1:i*m)*(x_e(:,i)-inv(weight((i-1)*m+1:i*m,(i-1)*m+1:i*m))*x(:,i)));
+            contri2(:,i) = A2(:,(i-1)*m+1:i*m)*(x_e(:,i));
+       	end
+		out2 = re_err;
+%err(iter) =norm(contri1-contri2,'fro');
+
+
+a1 = 0;
+a2 = 0;
+for i = 1:32
+    a1 = a1 +norm(contri1(:,i));
+    a2 = a2 +norm(contri1(:,i));
+end
+
+
+%         cvx_begin
+%         variable x_e(m,D) 
+%         norm21 = 0;
+%         for i = 1:D
+%             norm21 = norm21 + norm(A2(:,(i-1)*m+1:i*m)*x_e(:,i));
+%         end
+%         minimize(norm21)
+%         subject to
+%         A2*x_e(:) == y;
+%         cvx_end
+%         re_err = 0;
+%         for i = 1:D
+%             re_err = re_err+norm(A2(:,(i-1)*m+1:i*m)*(x_e(:,i)-inv(weight((i-1)*m+1:i*m,(i-1)*m+1:i*m))*x(:,i)));
+%             contri2(:,i) = A2(:,(i-1)*m+1:i*m)*(x_e(:,i));
+%        	end
+% 		out2 = re_err;
+%%
+n = 40;
+d = 2;
+D = 32;
+m = 4;
+%t = 1;
+k = 4;
+A = zeros(n,D*m);
+contri1 = zeros(n,D);
+contri2 = zeros(n,D);
+        for i = 0:32-1
+        %生成高斯基集
+	    set = randn(n,d);            
+        %生成随机的模为1的向量，与基集相乘得到块不满秩的高斯矩阵
+	    theta = rand(1,4)*2*pi;
+	    temp = [sin(theta);cos(theta)];
+            A(:,m*i+1:m*i+m) = set*temp*(i+1);
+        end
+        x = zeros(m,D);
+        col = randperm(D,k);
+        x(:,col) = randn(m,k);
+        y = A * x(:);
+
+        A2 = A * diag(randn(D*m,1)*30);
+        
+        %凸优化利用特殊l21范数求解恢复问题
+        cvx_begin
+        variable x_e(m,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        A*x_e(:) == y;
+        cvx_end
+        re_err = 0;
+        for i = 1:D
+            re_err = re_err+norm(A(:,(i-1)*m+1:i*m)*(x_e(:,i)-x(:,i)));
+            contri1(:,i) = A(:,(i-1)*m+1:i*m)*(x_e(:,i));
+       	end
+		out1 = re_err;
+
+
+
+        cvx_begin
+        variable x_e(m,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        A2*x_e(:) == y;
+        cvx_end
+        re_err = 0;
+        for i = 1:D
+            re_err = re_err+norm(A2(:,(i-1)*m+1:i*m)*(x_e(:,i)-inv(weight((i-1)*m+1:i*m,(i-1)*m+1:i*m))*x(:,i)));
+            contri2(:,i) = A2(:,(i-1)*m+1:i*m)*(x_e(:,i));
+       	end
+		out2 = re_err;
+%%
+n = 40;
+d = 2;
+D = 32;
+m = 4;
+%t = 1;
+k = 6;
+A = zeros(n,D*m);
+noise = 0.01;
+contri1 = zeros(n,D);
+contri2 = zeros(n,D);
+        for i = 0:32-1
+        %生成高斯基集
+	    set = randn(n,d);            
+        %生成随机的模为1的向量，与基集相乘得到块不满秩的高斯矩阵
+	    theta = rand(1,4)*2*pi;
+	    temp = [sin(theta);cos(theta)];
+            A(:,m*i+1:m*i+m) = set*temp*(i+1);
+        end
+        x = zeros(m,D);
+        col = randperm(D,k);
+        x(:,col) = randn(m,k);
+        y = A * x(:) + randn(n,1)*noise;
+        
+        weight =  diag(rand(D*m,1)*30);
+        %A2 = zeros(n,D*m);
+        %for i= 0:31
+        %    A2(:,m*i+1:m*i+m) = A(:,m*i+1:m*i+m) * weight(i+1);
+        %end
+        A2 = A*weight;
+        
+        %凸优化利用特殊l21范数求解恢复问题
+        cvx_begin
+        variable x_e(m,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(A(:,(i-1)*m+1:i*m)*x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+        norm(A*x_e(:) - y) <= sqrt(40)*0.01;
+        cvx_end
+        re_err = 0;
+        for i = 1:D
+            re_err = re_err+norm(A(:,(i-1)*m+1:i*m)*(x_e(:,i)-x(:,i)));
+            contri1(:,i) = A(:,(i-1)*m+1:i*m)*(x_e(:,i));
+       	end
+		out1 = re_err;
+
+
+        
+        cvx_begin
+        variable x_e(m,D) 
+        norm21 = 0;
+        for i = 1:D
+            norm21 = norm21 + norm(A2(:,(i-1)*m+1:i*m)*x_e(:,i));
+        end
+        minimize(norm21)
+        subject to
+       norm(A2*x_e(:) - y) <= sqrt(40)*0.01;
+        cvx_end
+        re_err = 0;
+        for i = 1:D
+            re_err = re_err+norm(A2(:,(i-1)*m+1:i*m)*(x_e(:,i)-inv(weight((i-1)*m+1:i*m,(i-1)*m+1:i*m))*x(:,i)));
+            contri2(:,i) = A2(:,(i-1)*m+1:i*m)*(x_e(:,i));
+       	end
+		out2 = re_err;
+

0 - example/lyh - 非满秩相变/theoretic.m

@@ -0,0 +1,7 @@
+function n = theoretic(m,s,d)
+    syms t;
+    syms u;
+    f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
+    g = diff(f,t);
+    t1 = solve(g);
+    n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);

0 - example/lyh - 非满秩相变/untitled.m

@@ -0,0 +1,147 @@
+
+len = 256;
+w = 0.05;
+m = [0:len-1];
+n = m;
+B = m - n';
+B = (sin(2*pi*w*B)./(pi*B));
+for i = 1:len
+    B(i,i) = 2*w;
+end
+t = [0:0.1:25.5];
+r = 1/sqrt(2*pi*1)*exp(-t.^2/2*1);
+gau = toeplitz(r);
+
+
+
+
+%%
+len = 64;
+w = 0.2;
+m = [0:len-1];
+n = m;
+B = m - n';
+B_new = sin(2*pi*w*B)./(pi*B);
+for i = 1:len
+    B_new(i,i) = 2*w;
+end
+
+cor = 0.236;
+
+final = zeros(len^2,len^2);
+for j = 1:len
+    for k = 1:len
+        tmp = sin(2*pi*w*(cor*B+(j-k)))./(pi*(cor*B+(j-k)));
+        tmp(cor*B+(j-k) == 0) = 2*w;
+        final((j-1)*len+1:j*len,(k-1)*len+1:k*len) = tmp.*B_new;
+    end
+end
+
+[a,b] = eig(final);
+plot(diag(b));
+c = diag(b);
+M = c'*c;
+plot(sort(M(:)));
+
+
+
+%%
+
+K = kron(B,B);
+
+[a,b] = eig(B);
+plot(diag(b));
+c = diag(b);
+M = c*c';
+plot(sort(M(:)));
+
+
+f = zeros(1,1200);
+f(1:120) = randn(120,1);
+f(1081:1200) = randn(120,1);
+t = ifft(f);
+t = t(1:64);
+
+
+
+
+
+%B = B + diag(exp(1i*2*pi*0.3*[0:len-1]))*B*diag(exp(1i*2*pi*0.3*[0:len-1]))';
+
+[a,b] = eig(B);
+% plot(sort(abs(diag(b))))
+
+D = B(1:8,:);
+E = B(1:2:16,:);
+[U2,S2,V2] = svd(D);
+[U3,S3,V3] = svd(E);
+
+
+
+
+
+
+
+
+%sel = (1:2:1023);
+o = randn(1024,1024);
+%o = orth(o);
+sel = randperm(1024,512);
+% C = o(sel,:)*b;
+% D = B(sel,:);
+% E = o(sel,randperm(1024,204));
+D = B(1:512,:);
+C = B(sel,:);
+E = B(1:2:1023,:);
+[U,S,V] = svd(C);
+[U2,S2,V2] = svd(D);
+[U3,S3,V3] = svd(E);
+s = diag(S);
+s2 = diag(S2);
+s3 = diag(S3);
+figure;
+subplot(1,4,1);
+plot(s);
+subplot(1,4,2);
+plot(s2);
+subplot(1,4,3);
+plot(s3);
+q = sum(s);
+subplot(1,4,4);
+plot(flipud(abs(diag(b))));
+
+
+%t1 = squeeze(resulta(:,:,1));
+%t2 = squeeze(resultb(:,:,1));
+% t1 = squeeze(resulta(:,:,2));
+% t2 = squeeze(resultb(:,:,2));
+% figure;
+
+result_cona = reshape(result_cona,200,20,2);
+result_conb = reshape(result_conb,200,20,2);
+for i = 1:2
+t1 = squeeze(result_cona(:,:,i)+t3(:,:,i))/2;
+t2 = squeeze(result_conb(:,:,i)+t4(:,:,i))/2;
+figure
+hold on
+plot(mean(t1),'-r.');
+plot(mean(t2),'-bo');
+h = legend("$P_{\ell_{2,1}}'$","$P_{\ell_{2,1}}$","Location","Southeast","Fontsize",15);
+set(h,'Interpreter','latex');
+xlabel("\fontname{Times New Roman} Block Sparsity \it s_B");
+ylabel("\fontname{Times New Roman} Block Contribution Error");
+end
+
+
+xlabel("\fontname{Times New Roman} Block Sparsity \it K");
+ylabel("\fontname{Times New Roman} Block Contribution Error");
+
+
+t3 = result_cona;
+t4 = result_conb;
+
+
+
+
+
+

0 - example/nsq - 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

0 - example/nsq - 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]);
+
+
+
+
+
+
+
+
+
+

0 - example/nsq - 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;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

0 - example/nsq - 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

0 - example/nsq - 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]);
+
+
+
+
+
+
+
+
+
+

0 - example/nsq - 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;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

0 - example/nsq - 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

0 - example/nsq - 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]);
+
+
+
+
+
+
+
+
+
+

0 - example/nsq - 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;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

0 - example/nsq - 恢复算法/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 = sum(abs(z_pre - z))/sum(abs(z)); 
+    x_pre = x;
+    z_pre = z;
+    t_pre = t;
+    k = k + 1;
+end
+
+end

0 - example/nsq - 恢复算法/cVAMPa_dampling.m

@@ -0,0 +1,79 @@
+function[x, x_d, hat_Q1, sigma_d, ifcvg] = cVAMPa_dampling(y, A, lambda, alpha, delta, iter_max, sigma)
+
+[M, N] = size(A);
+
+gamma = M/N;
+p = A'*y;
+h1 = p/gamma;
+hat_Q1 = gamma;
+
+%Eigenvalue Decomposition
+[V, D] = eig(A'*A);
+d = diag(D); 
+
+t = 0;
+diff = 1;
+
+while((diff > delta) && (t < iter_max))
+    
+    h1_pre = h1;
+    
+    % Factorized
+    x1 = sft_thd(h1, lambda/hat_Q1);
+    chi1 = sum((abs(x1) > 1e-4).* (2 - lambda./((hat_Q1*abs(x1) + lambda)))) / 2 / N / hat_Q1;
+    % Message F to G
+    hat_Q2 = 1/chi1 - hat_Q1;
+    h2 = (x1/chi1 - h1*hat_Q1)/hat_Q2;
+    % Gaussian
+    tmp = V'*(p + h2*hat_Q2);
+    tmp = tmp./(d+hat_Q2);
+    x2 = V*tmp;
+    chi2 = sum(1./(d+hat_Q2))/N;
+    % Message G to F
+    hat_Q1 = 1/chi2 - hat_Q2;
+    h1 = alpha*(x2/chi2 - h2*hat_Q2)/hat_Q1+(1-alpha)*h1;
+    
+    diff = sum(abs(h1_pre - h1))/sum(abs(h1));
+    t = t+1;
+    
+end
+
+ifcvg = diff <= delta;
+
+x = x1;
+x_d = h1;
+
+chi = chi1;
+t = -hat_Q2;
+t_prime = -1/mean((1./(d+hat_Q2)).^2);
+G_prime = t + 1/chi;
+G_wprime = t_prime + 1/chi/chi;
+RSS = sum(abs(y - A*x).^2)/M;
+hat_chi = gamma*G_wprime/(2*G_prime-2*chi*G_wprime)*RSS +...
+    (-G_wprime*gamma+G_prime*G_prime)/(2*G_prime-2*chi*G_wprime)*sigma^2;
+sigma_d = sqrt(2*hat_chi)/hat_Q1;
+
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

0 - example/nsq - 恢复算法/cal_debiased_LASSO.m

@@ -0,0 +1,37 @@
+function [x_d, hat_Q1, sigma_d] = cal_debiased_LASSO(x, A, y, lambda, sigma)
+
+[M, N] = size(A);
+gamma = M/N;
+hat_Q1 = gamma;
+[~, D] = eig(A'*A);
+d = diag(D);
+
+diff = 1;
+T = 1000;
+t = 0;
+
+while (t < T) && (diff > 1e-6)
+    
+    Q1_pre = hat_Q1;
+    rho = mean((2 - lambda./(hat_Q1*abs(x) + lambda)).*(abs(x) > 1e-4))/2;
+    hat_Q1 = rho/mean(1./(d + (1-rho)*hat_Q1/rho));
+    diff = abs(Q1_pre - hat_Q1);
+    t = t+1;
+
+end
+
+x_d = x + 1/hat_Q1*A'*(y - A*x);
+
+chi = rho/hat_Q1;
+hat_Q2 = 1/chi - hat_Q1;
+
+t = -hat_Q2;
+t_prime = -1/mean((1./(d+hat_Q2)).^2);
+G_prime = t + 1/chi;
+G_wprime = t_prime + 1/chi/chi;
+RSS = sum(abs(y - A*x).^2)/M;
+hat_chi = gamma*G_wprime/(2*G_prime-2*chi*G_wprime)*RSS +...
+    (-G_wprime*gamma+G_prime*G_prime)/(2*G_prime-2*chi*G_wprime)*sigma^2;
+sigma_d = sqrt(2*hat_chi)/hat_Q1;
+
+end

0 - example/nsq - 恢复算法/generate_matrix_new.m

@@ -0,0 +1,18 @@
+function [ A ] = generate_matrix_new(y1, y2)
+    A = [];
+    l = length(y1);
+    temp1 = y1;
+    temp2 = y2;
+    for i = 1:l
+        t1 = circshift(temp1, i-1);
+        t2 = circshift(temp2, i-1);
+        if i - 1 > 0
+            t1(1:i - 1,1) = 0;
+        end
+        if i - 1 > 0
+            t2(1:i - 1,1) = 0;
+        end
+        A = [A,t1,t2];
+    end
+end
+

0 - example/nsq - 恢复算法/model.m

@@ -0,0 +1,146 @@
+clc;
+clear;
+close all;
+
+%% 参数设置
+sigma_n = 0.1;
+gamma = 0.5;
+
+% 信号参数
+B = 5e5;        %信号带宽
+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;
+
+target_scattering = [0.8, 1, 0.9]; %扩展目标各点散射强度
+
+
+%% 生成矩阵 A
+% 生成发射信号 signal_t 及 
+N = Tr * fs;
+N_high = Tp * fs;
+
+signal_t  = zeros(1, N);
+signal_td  = 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));
+    tp2 = (i - 0.5) * (1 / fs);
+    signal_td(1, i + 1) = exp(1j*2*pi*(fc*tp2+0.5*K*tp2.^2));
+end
+
+A = generate_matrix_new(signal_t.', signal_td.');
+
+J = A'*A;
+lambda_J=eig(J);
+
+figure(1)
+spy(A)
+
+
+%% 生成回波 y
+% 设置目标 - 扩展目标，由三个点组成
+distance1 = 52000;
+tau1 = distance1 * 2 / c;
+n_tau1 = round(tau1 * fs);
+
+alpha1 = 0.6; % 扩展目标整体散射强度
+signal_r1_t1  = zeros(1, N);
+signal_r2_t1  = zeros(1, N);
+signal_r3_t1  = zeros(1, N);
+
+for i = 1:N
+    temp = i - n_tau1;
+    if temp >= 1 && temp <= N_high
+        signal_r1_t1(1, i) = alpha1 * target_scattering(1) * signal_t(1, temp);
+        if i + 1 <= N
+            signal_r2_t1(1, i + 1) = alpha1 * target_scattering(2) * signal_t(1, temp);
+        end
+        if i + 1 <= N
+            signal_r3_t1(1, i + 2) = alpha1 * target_scattering(3) * signal_t(1, temp);
+        end
+    end
+end
+signal_r1 = signal_r1_t1 + signal_r2_t1 + signal_r3_t1;
+
+% 回波
+signal_r = signal_r1;
+
+% 加入噪声
+noise = random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r)) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r));
+signal_r_n = signal_r + noise;
+
+y = signal_r_n.';
+figure(2)
+subplot(211);
+plot(t, real(signal_t));
+title('发射信号')
+xlabel('时间')
+ylabel('幅度')
+
+subplot(212);
+plot(t, real(y))
+title('接收信号(y)')
+xlabel('时间')
+ylabel('幅度')
+
+%% 理论 x
+distance_node = round((distance1 * 2 / c) * fs);
+x_t = zeros(1, 2 * N);
+for i = 1: length(target_scattering)
+    x_t((distance_node + i) * 2 - 1) = alpha1 * target_scattering(i);
+end
+x = x_t.';
+
+figure(3)
+plot(t2, x);
+title('目标散射点(x)')
+xlabel('时间')
+ylabel('幅度')
+
+
+
+%% 验证
+y_t = A * x;
+y_r = signal_r.';
+
+figure(4)
+subplot(311);
+plot(real(y_r))
+title('实际回波')
+
+subplot(312);
+plot(real(y_t));
+title('计算结果')
+
+subplot(313);
+plot(abs(y_r - y_t));
+title('差异')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

0 - example/nsq - 恢复算法/sft_thd.m

@@ -0,0 +1,20 @@
+function y = sft_thd(x, thd)
+
+if isequal(size(x), size(thd))
+    
+    tmp = abs(x);
+    y = x;
+    y(tmp <= thd) = 0;
+    y(tmp > thd) = (tmp(tmp > thd) - thd(tmp > thd)) .* x(tmp > thd) ./ tmp(tmp > thd);
+    
+else
+    
+    tmp = abs(x);
+    y = x;
+    y(tmp <= thd) = 0;
+    y(tmp > thd) = (tmp(tmp > thd) - thd) .* x(tmp > thd) ./ tmp(tmp > thd);
+    
+end
+
+
+end

0 - example/nsq - 恢复算法/test_model.m

@@ -0,0 +1,175 @@
+clc;
+clear;
+%close all;
+
+%% 参数设置
+sigma_n = 0.1;
+gamma = 0.5;
+
+% 信号参数
+B = 5e5;        %信号带宽
+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;
+
+target_scattering = [0.8, 1, 0.9]; %扩展目标各点散射强度
+
+
+%% 生成矩阵 A
+% 生成发射信号 signal_t 及 
+N = Tr * fs;
+N_high = Tp * fs;
+
+signal_t  = zeros(1, N);
+signal_td  = 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));
+    tp2 = (i - 0.5) * (1 / fs);
+    signal_td(1, i + 1) = exp(1j*2*pi*(fc*tp2+0.5*K*tp2.^2));
+end
+
+A = generate_matrix_new(signal_t.', signal_td.');
+
+temp = 0;
+for i = 1: size(A, 1)
+
+    for j = 1: size(A, 2)
+
+        temp = temp + abs(A(mod(i, size(A, 1))+1, mod(j+1, size(A, 2))+1) - A(i, j));
+
+    end
+
+end
+
+%% 生成回波 y
+% 设置目标 - 扩展目标，由三个点组成
+distance1 = 52000;
+tau1 = distance1 * 2 / c;
+n_tau1 = round(tau1 * fs);
+
+alpha1 = 0.6; % 扩展目标整体散射强度
+signal_r1_t1  = zeros(1, N);
+signal_r2_t1  = zeros(1, N);
+signal_r3_t1  = zeros(1, N);
+
+for i = 1:N
+    temp = i - n_tau1;
+    if temp >= 1 && temp <= N_high
+        signal_r1_t1(1, i) = alpha1 * target_scattering(1) * signal_t(1, temp);
+        if i + 1 <= N
+            signal_r2_t1(1, i + 1) = alpha1 * target_scattering(2) * signal_t(1, temp);
+        end
+        if i + 1 <= N
+            signal_r3_t1(1, i + 2) = alpha1 * target_scattering(3) * signal_t(1, temp);
+        end
+    end
+end
+signal_r1 = signal_r1_t1 + signal_r2_t1 + signal_r3_t1;
+
+% 回波
+signal_r = signal_r1;
+
+% 加入噪声
+noise = random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r)) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r));
+signal_r_n = signal_r + noise;
+
+y = signal_r_n.';
+
+%% 理论 x
+distance_node = round((distance1 * 2 / c) * fs);
+x_t = zeros(1, 2 * N);
+for i = 1: length(target_scattering)
+    x_t((distance_node + i) * 2 - 1) = alpha1 * target_scattering(i);
+end
+x = x_t.';
+
+%% Experiment
+
+%% Parameters setting
+lambda = 0.002;
+alpha = 1/4;
+delta = 1e-8*alpha;
+iter_max = round(2000/alpha);
+n = size(A, 2);
+
+% J = A'*A;
+J1 = A*A';
+lambda_J=eig(J1);
+% histogram(lambda_J, 100);
+
+%% Normalized
+A = A / sqrt(lambda_J(end));
+y = y / sqrt(lambda_J(end));
+sigma_n = sigma_n / sqrt(lambda_J(end));
+
+%% cVAMP
+tic;
+[x_VAMP, x_d, hat_Q1, sigma_d, ifcvg] = cVAMPa_dampling(y, A, lambda, alpha, delta, iter_max, sigma_n);
+toc;
+
+%% cvx
+tic;
+cvx_begin quiet
+    variable x_cvx(n, 1) complex
+    z = lambda*sum(abs(x_cvx)) + 0.5*sum(pow_abs((y - A * x_cvx), 2));
+    minimize(z)
+cvx_end
+
+[x_d_cal, hat_Q1_cal, sigma_d_cal] = cal_debiased_LASSO(x_cvx, A, y, lambda, sigma_n);
+toc;
+
+sigma_ex = std(x_d_cal - x, 1);
+tmp = (x_d_cal - x)/sigma_ex;
+
+[h_r, p_r, k_r, c_r] = kstest(real(tmp)*sqrt(2));
+[h_i, p_i, k_i, c_i] = kstest(imag(tmp)*sqrt(2));
+
+%% results
+% whether cVAMP algorithm converges
+ifcvg 
+
+% whether the output of cVAMP converges to the LASSO solution
+sum(abs(x_cvx - x_VAMP))
+
+% check the results from cVAMP and "calculation"
+abs(hat_Q1 - hat_Q1_cal)
+abs(sigma_d - sigma_d_cal)
+sum(abs(x_d - x_d_cal))
+
+% accuracy of estimating the variance
+abs(sigma_d_cal - sigma_ex)/abs(sigma_ex)
+
+% p-value of KS-test 
+% the larger, the higher probability it is drawn from Gaussian distribution
+p_r
+p_i
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

0 - example/wzk - PSK_ROC/test_function.m

@@ -0,0 +1,69 @@
+
+% load st.mat
+% [ A ] = generate_matrix_long(transpose(s_T));
+
+
+N = 512;
+A = dftmtx(N) / sqrt(N);
+A = A(1: 256, :);
+
+Target_index = 321;
+SNR = 10;
+lambda = 0.0005;
+
+%%
+[ x_mf, x_cs ] = yAxn_recovery( A, SNR, Target_index, lambda );
+
+figure(1111)
+subplot(211)
+plot(abs(x_mf))
+subplot(212)
+plot(abs(x_cs))
+
+
+%%
+rep_time = 500;
+[ res_MF, res_CS, P_fa, H1_MF_cnt, H0_MF_cnt, H1_CS_cnt, H0_CS_cnt ] = yAxn_MC( A, SNR, Target_index, lambda, rep_time );
+
+P_fa_MF = res_MF(:,1);
+P_d_MF = res_MF(:,2);
+P_fa_CS = res_CS(:,1);
+P_d_CS = res_CS(:,2);
+
+figure;
+loglog(P_fa,P_fa_MF, 'linewidth', 2);
+hold on;
+loglog(P_fa,P_fa_CS, 'linewidth', 2);
+legend('MF','CS')
+xlabel('P_{fa}');
+ylabel('Actual P_{fa}');
+
+figure;
+semilogx(P_fa,P_d_MF, 'linewidth', 2);
+hold on;
+semilogx(P_fa,P_d_CS, 'linewidth', 2);
+legend('MF','CS')
+xlabel('P_{fa} set');
+ylabel('P_{d}');
+
+figure;
+semilogx(P_fa_MF,P_d_MF, 'linewidth', 2);
+hold on;
+grid on;
+semilogx(P_fa_CS,P_d_CS, 'linewidth', 2);
+xlim([1e-4,1])
+legend('MF','CS')
+xlabel('Actual P_{fa}');
+ylabel('P_{d}');
+title('ROC');
+
+
+figure
+subplot(211)
+histogram(real(H0_MF_cnt(5,:)))
+hold on;
+histogram(real(H1_MF_cnt))
+subplot(212)
+histogram(real(H0_CS_cnt(5,:)))
+hold on;
+histogram(real(H1_CS_cnt))

0 - example/wzk - PSK_ROC/test_signal_t.m

@@ -0,0 +1,57 @@
+clc
+clear
+close all
+
+rng(6)
+
+PRF = 5000;         % 脉冲重复频率
+Tr = 1 / PRF;       % 脉冲重复间隔
+Tp = 2e-5;
+P = 10;
+B = P / Tp;
+fs = 10 * B;
+Ts = 1/fs;
+fc = 1.25e9;
+
+t = 0:1/fs:1-1/fs;
+
+tsin = Tp / P * 2;
+fsin = 1 / tsin;
+
+N = round(Tr * fs);
+N_high = round(Tp * fs / P);
+
+signal_t  = zeros(1, N);
+signal_td  = zeros(1, N);
+
+sign_p = sign(randn(1, P));
+
+for p = 1: P
+    for i = 1: N_high
+            tp = (i - 1) * (1 / fs);
+            signal_t(1, (p-1)*N_high + i) = sign_p(p) * exp(1j * 2 * pi * fsin * tp);
+            tp2 = (i - 0.5) * (1 / fs);
+            signal_td(1, (p-1)*N_high + i) = sign_p(p) * exp(1j * 2 * pi * fsin * tp2);
+    end
+end
+
+A = zeros(N, 2 * N);
+A(:, 1) = transpose(signal_t);
+A(:, 2) = transpose([0, signal_td(1: end-1)]);
+for k = 3: 2 * N
+    A(:, k) = [0; A(1: end - 1, k - 2)];
+end
+
+
+
+figure(1)
+subplot(211)
+plot(real(A(1:N_high*P+20, 1)))
+hold on;
+plot(real(A(1:N_high*P+20, 2)))
+plot(real(A(1:N_high*P+20, 3)))
+subplot(212)
+plot(imag(A(1:N_high*P+20, 1)))
+hold on;
+plot(imag(A(1:N_high*P+20, 2)))
+plot(imag(A(1:N_high*P+20, 3)))

0 - example/wzk - PSK_ROC/yAxn_MC.m

@@ -0,0 +1,106 @@
+function [ res_MF, res_CS, P_fa, H1_MF_cnt, H0_MF_cnt, H1_CS_cnt, H0_CS_cnt ] = yAxn_MC( A, SNR, Target_index, lambda, rep_time )
+
+    % matrix parameters
+    [M, N] = size(A);
+    mul = A(:,1)'*A(:,1);
+    
+    % H0 sample
+    L = round(Target_index + N * 0.1);
+    R = round(N * 0.9);
+    Lambda_C = L: R;
+    
+    % parameters
+    sigma_n = 0.1;
+    Target_amplitude = sqrt(10^(SNR/10) * sigma_n^2 / mul);
+
+    
+    %% CS setting
+    % CS-parameters
+%     lambda = 0.0005;
+    alpha = 1/4;
+    delta = 1e-8*alpha;
+    
+    % CS-normalization
+    J1 = A*A';
+    lambda_J=eig(J1);
+    A_norm = A / sqrt(lambda_J(end));
+    Hp = A_norm'*A_norm;
+    [~, D] = eig(Hp);
+    d = diag(D);
+    sigma_n_norm = sigma_n / sqrt(lambda_J(end));
+    
+    
+    %% generate x
+    x = zeros(N, 1);
+    x(Target_index) = Target_amplitude;
+    
+    
+    %% MC parameters
+    P_fa = [1e-6, 1e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1];
+    len_P_fa = length(P_fa);
+    
+    P_fa_MF_cnt = zeros(len_P_fa, rep_time);
+    P_d_MF_cnt = zeros(len_P_fa, rep_time);
+    P_fa_CS_cnt = zeros(len_P_fa, rep_time);
+    P_d_CS_cnt = zeros(len_P_fa, rep_time);  
+    
+    H1_index = zeros(N, 1);
+    H1_index(Target_index) = 1;
+    H0_index = ones(size(x));
+    H0_index(Target_index) = 0;
+    H1_MF_cnt = zeros(1, rep_time);
+    H0_MF_cnt = zeros(N - 1, rep_time);
+    H1_CS_cnt = zeros(1, rep_time);
+    H0_CS_cnt = zeros(N - 1, rep_time);
+    
+    %% MC
+    parfor rep = 1: rep_time
+%     for rep = 1: rep_time
+        %% generate 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;
+
+        %% recover
+        % MF
+        x_mf = A' * y ./ mul;
+        % sigma_MF = sqrt(var(x_mf(Lambda_C)));
+        sigma_MF = sigma_n / sqrt(mul);
+        x_mf_norm = x_mf ./ sigma_MF;
+        stat_MF = abs(x_mf ./ sigma_MF).^2;
+
+        % CS
+        y_norm = y / sqrt(lambda_J(end));
+        x_FISTA = FISTA_v1(y_norm, A_norm, lambda, delta, Hp);
+        [x_d_cal_f, hat_Q1_cal_f, sigma_d_cal_f] = cal_debiased_LASSO_v1(x_FISTA, A_norm, y_norm, lambda, sigma_n_norm, d);
+        % sigma_CS = sqrt(var(x_d_cal_f(Lambda_C)));
+        sigma_CS = sigma_d_cal_f;
+        x_cs_norm = x_d_cal_f ./ sigma_CS;
+        stat_CS = abs(x_d_cal_f ./ sigma_CS).^2;
+
+        %% detect
+        H1_MF_cnt(rep) = x_mf(Target_index);
+        H0_MF_cnt(:,rep) = [x_mf(1:Target_index-1);x_mf(Target_index+1:end)];
+        H1_CS_cnt(rep) = x_d_cal_f(Target_index);
+        H0_CS_cnt(:,rep) = [x_d_cal_f(1:Target_index-1);x_d_cal_f(Target_index+1:end)];
+        
+        kd = chi2inv(1 - P_fa, 2) / 2;
+        for cnt_h_th = 1: len_P_fa
+            P_fa_MF_cnt(cnt_h_th, rep) = sum(stat_MF(H0_index > 0) > kd(cnt_h_th)) / sum(H0_index);
+            P_d_MF_cnt(cnt_h_th, rep) = sum(stat_MF(H1_index > 0) > kd(cnt_h_th)) / sum(H1_index);
+            P_fa_CS_cnt(cnt_h_th, rep) = sum(stat_CS(H0_index > 0) > kd(cnt_h_th)) / sum(H0_index);
+            P_d_CS_cnt(cnt_h_th, rep) = sum(stat_CS(H1_index > 0) > kd(cnt_h_th)) / sum(H1_index);
+        end
+        
+        fprintf('%d\n', rep);
+    end
+    
+    P_fa_MF = mean(P_fa_MF_cnt, 2);
+    P_d_MF = mean(P_d_MF_cnt, 2);
+    P_fa_CS = mean(P_fa_CS_cnt, 2);
+    P_d_CS = mean(P_d_CS_cnt, 2);
+    
+    res_MF = [P_fa_MF, P_d_MF];
+    res_CS = [P_fa_CS, P_d_CS];
+    
+end
+

0 - example/wzk - PSK_ROC/yAxn_recovery.m

@@ -0,0 +1,57 @@
+function [ x_mf_norm, x_cs_norm ] = yAxn_recovery( A, SNR, Target_index, lambda )
+
+    % matrix parameters
+    [M, N] = size(A);
+    mul = A(:,1)'*A(:,1);
+    
+    % H0 sample
+    L = round(Target_index + N * 0.1);
+    R = round(N * 0.9);
+    Lambda_C = L: R;
+    
+    % parameters
+    sigma_n = 0.1;
+    Target_amplitude = sqrt(10^(SNR/10) * sigma_n^2 / mul);
+
+    
+    %% CS setting
+    % CS-parameters
+%     lambda = 0.0005;
+    alpha = 1/4;
+    delta = 1e-8*alpha;
+    
+    % CS-normalization
+    J1 = A*A';
+    lambda_J=eig(J1);
+    A_norm = A / sqrt(lambda_J(end));
+    Hp = A_norm'*A_norm;
+    [~, D] = eig(Hp);
+    d = diag(D);
+    sigma_n_norm = sigma_n / sqrt(lambda_J(end));
+    
+    
+    %% generate x
+    x = zeros(N, 1);
+    x(Target_index) = Target_amplitude;
+    
+    
+    %% generate 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;
+    
+    
+    %% recover
+    % MF
+    x_mf = A' * y ./ mul;
+    sigma_MF = sqrt(var(x_mf(Lambda_C)));
+    x_mf_norm = x_mf ./ sigma_MF;
+    
+    % CS
+    y_norm = y / sqrt(lambda_J(end));
+    x_FISTA = FISTA_v1(y_norm, A_norm, lambda, delta, Hp);
+    [x_d_cal_f, hat_Q1_cal_f, sigma_d_cal_f] = cal_debiased_LASSO_v1(x_FISTA, A_norm, y_norm, lambda, sigma_n_norm, d);
+    sigma_CS = sqrt(var(x_d_cal_f(Lambda_C)));
+    x_cs_norm = x_d_cal_f ./ sigma_CS;
+    
+end
+

0 - example/wzk - RDA/angle_estimation.m

@@ -0,0 +1,60 @@
+function [ target_list_RDA ] = angle_estimation(target_list_RDA_temp, echo_mtx, A, d, lambda, acc)
+% This function estimates all targets angle
+%
+% Usage: 
+% [target_list_RDA] = angle_estimation(target_list_RDA_temp, echo_mtx, A, d, lambda, acc)
+%
+% Inputs:
+%   target_list_RDA_temp:   Temporary list of targets before angle estimation
+%   echo_mtx:               Original echo data matrix
+%   A:                      Chirp measurement matrix
+%   d:                      Antenna spacing
+%   lambda:                 Wave length
+%   acc:                    Accuracy for angle estimation
+%
+% Outputs:
+%   target_list_RDA:        Target list include range, doppler and angle information
+
+
+    % parameters
+    if nargin < 6
+        acc = 181;
+    end
+    
+    % calculate all targets angle
+    target_list_RDA = zeros(size(target_list_RDA_temp));
+    for i = 1: size(target_list_RDA_temp, 1)
+        t_rda = target_list_RDA_temp(i, :);
+        target_list_RDA(i, 1:2) = t_rda(1:2);
+        target_list_RDA(i, 3) = cal_angle(t_rda(1), t_rda(2), echo_mtx, A, d, lambda, acc);
+    end
+
+end
+
+
+%% sub-functions
+% calculate target actual angle with range and doppler index
+function [ angle ] = cal_angle(r_idx, d_idx, echo_mtx, A, d, lambda, acc)
+
+    % range doppler processing
+    echo_r = zeros([size(echo_mtx, 1), 1, size(echo_mtx, 3)]);
+    for i = 1 : size(echo_mtx, 3)
+        echo_r(:, 1, i) = sum(bsxfun(@times, echo_mtx(:, :, i), A(:, r_idx)'), 2);
+    end
+    echo_r = squeeze(echo_r);
+    echo_r_ex = [echo_r, zeros(size(echo_r))];
+    echo_rd = zeros(size(echo_r_ex, 1), 1);
+    for i = 1 : size(echo_r_ex, 1)
+         temp_v = fftshift(fft(echo_r_ex(i,:)));
+         echo_rd(i) = temp_v(d_idx);
+    end
+    
+    % CBF
+    THETA = linspace(-90, 90, round(acc));
+    fai = exp((0: size(echo_r_ex, 1) - 1)' *...
+        -1j * 2 * pi / lambda * d * sin(THETA / 180 * pi));
+    angle_mf = abs(fai'* echo_rd);
+    [v, idx] = max(angle_mf);
+    angle = THETA(idx);
+
+end

0 - example/wzk - RDA/doppler_process_CS.m

@@ -0,0 +1,200 @@
+function [ echo_rd_mtx, stat_RD, sigma_n_o ] = doppler_process_CS(echo_r_mtx, sigma_n, lambda, gamma, delta)
+% Doppler processing for echo data using compressed sensing
+%
+% Usage: 
+% [echo_rd_mtx, stat_RD, sigma_n_o] = doppler_process_CS(echo_r_mtx, sigma_n, lambda, gamma, delta)
+%
+% Inputs:
+%   echo_r_mtx:     Range processing echo data
+%   sigma_n:        Input noise standard deviation
+%   lambda:         LASSO weight
+%   gamma:          Compressed ratio(0.5)
+%   delta:          Convergence normalized difference(1e-6)
+%
+% Outputs:
+%   echo_rd_mtx:    Range-doppler processing echo data
+%   stat_RD:        Range-doppler statistics
+%   sigma_n_o:      Output noise standard deviation
+
+
+    % parameters
+    if nargin < 5
+        delta = 1e-6;
+    end
+    if nargin < 4
+        gamma = 0.5;
+    end
+    [lenA, lenR, M] = size(echo_r_mtx);
+    N = round(M / gamma);
+    iter_max_VAMP = 1000;  
+    lambda_v = zeros(N, 1) + lambda;
+    
+    % generate mtx
+    F_ori = dftmtx(N);
+    F = F_ori(1:M,:);
+    F_inv = conj(F) / N;
+    
+    % normalization
+    A = (sqrt(N) * eye(M)) * F_inv;
+    echo_r_mtx = sqrt(N) .* echo_r_mtx;
+    sigma_n = sqrt(N) * sigma_n;
+    
+    % doppler processing
+    echo_rd_mtx = zeros(lenA, lenR, N);
+    stat_RD = zeros(lenA, lenR, N);
+    sigma_n_o_cnt = zeros(lenA, lenR);
+    for numA = 1: lenA
+        parfor numR = 1: lenR
+            sample = squeeze(echo_r_mtx(numA, numR, :));
+            y = sample;
+            x_LASSO = cVAMPro(y, A, lambda_v, delta, iter_max_VAMP);
+            [x_d_CROD, sigma_CROD] = CROD(y, A, x_LASSO, lambda, sigma_n);
+            sigma_n_o_cnt(numA, numR) = abs(sigma_CROD);
+            stat_RD(numA, numR, :) = abs(fftshift(x_d_CROD) / sigma_CROD).^2;
+            echo_rd_mtx(numA, numR, :) = fftshift(x_d_CROD);
+        end
+    end
+%     sigma_n_o = mean(sigma_n_o_cnt, 2);
+    sigma_n_o = mean(mean(sigma_n_o_cnt));
+    
+end
+
+
+%% sub-functions
+% algorithm for LASSO
+% y: measurements
+% A: measurement matrix
+% lambda: LASSO weight
+% tau: convergence normalized difference
+% Kit: maximum number of iterations
+% LASSO estimator
+function x_hat_wl = 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;
+        % next
+        h_1 = h_1_next;
+        Q_1 = Q_1_next;
+    end
+
+end
+
+
+% soft threshold function
+% x: processing object
+% thd: threshold
+% y: result
+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
+
+% calculate debiased LASSO estimator
+% y: measurements
+% A: measurement matrix
+% x_LASSO: LASSO estimator
+% lambda: LASSO weight
+% sigma_n: input noise standard deviation
+% x_d_CROD: debiased LASSO estimator
+% sigma_CROD: equivalent noise standard deviation estimator
+function [ x_d_CROD, sigma_CROD ] = CROD(y, A, x_LASSO, lambda, sigma_n)
+
+    [m, n] = size(A);
+    gamma = m / n;
+
+    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;
+
+end

0 - example/wzk - RDA/doppler_process_MF.m

@@ -0,0 +1,45 @@
+function [ echo_rd_mtx, stat_RD, sigma_n_o ] = doppler_process_MF(echo_r_mtx, sigma_n, gamma)
+% Doppler processing for echo data using matching filter
+%
+% Usage: 
+% [echo_rd_mtx, stat_RD, sigma_n_o] = doppler_process_MF(echo_r_mtx, sigma_n, gamma)
+%
+% Inputs:
+%   echo_r_mtx:     Range processing echo data
+%   sigma_n:        Input noise standard deviation
+%   gamma:          Compressed ratio(0.5)
+%
+% Outputs:
+%   echo_rd_mtx:    Range-doppler processing echo data
+%   stat_RD:        Range-doppler statistics
+%   sigma_n_o:      Output noise standard deviation
+
+
+    % parameters
+    if nargin < 3
+        gamma = 0.5;
+    end
+    [lenA, lenR, M] = size(echo_r_mtx);
+    N = round(M / gamma);
+    
+    % generate mtx
+    F_ori = dftmtx(N);
+    F = F_ori(1:M,:);
+    multiple_d = F(:,1)' * F(:,1);
+    
+    % doppler matched filtering
+    echo_rd_mtx = zeros(lenA, lenR, N);
+    for numA = 1: lenA
+        for numR = 1: lenR
+            sample = squeeze(echo_r_mtx(numA, numR, :));
+            dpl_temp = transpose(F) * sample;
+            dpl_norm = fftshift(dpl_temp ./ multiple_d);
+            echo_rd_mtx(numA, numR, :) = dpl_norm;
+        end
+    end
+    
+    % calculate output noise
+    sigma_n_o = sqrt(sigma_n^2 / multiple_d);
+    stat_RD = abs(echo_rd_mtx ./ sigma_n_o).^2;
+    
+end

0 - example/wzk - RDA/echo_processing_CS.m

@@ -0,0 +1,56 @@
+function [ echo_rd_mtx, target_list_RDA, sigma_n_rd ] = ...
+    echo_processing_CS( echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d )
+% This function processes echo data using compressed sensing method.
+%
+% Usage: 
+% [echo_rd_mtx, target_list_RDA, sigma_n_rd] = echo_processing_MF(echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d)
+% 
+% Inputs:
+%   echo_mtx:       Matrix containing the echo data
+%   PRF:            Pulse Repetition Frequency
+%   fs:             Sampling frequency
+%   fc:             Carrier frequency
+%   B:              Bandwidth
+%   D:              Duty ratio
+%   d:              Antenna spacing
+%   sigma_n:        Input noise level of the echo data
+%   P_fa:           False alarm probability threshold
+%   lambda_r:       LASSO weight for range processing
+%   lambda_d:       LASSO weight for doppler processing
+%
+% Outputs:
+%   echo_rd_mtx:        Echo matrix after range-Doppler processing
+%   target_list_RDA:    List of detected targets
+%   sigma_n_rd:         Estimated noise level after range-Doppler processing
+
+
+    % parameters
+    c = 3e8;
+    lambda = c / fc;
+    [num_antenna, N, num_pluse] = size(echo_mtx);
+    
+    % data processing
+    [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pluse);
+    
+    [echo_sFFT_mtx, angle_v] = spatial_FFT(echo_mtx, d, lambda);
+    fprintf('   Spatial_FFT done.\n');
+    
+    [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D);
+    fprintf('   Generate_chirp_mtx done.\n');
+    
+    [echo_r_mtx, sigma_n_o] = range_process_CS(echo_sFFT_mtx, A, sigma_n, lambda_r);
+    fprintf('   Range_processing done.\n');
+    
+    [echo_rd_mtx, stat_RD, sigma_n_rd] = doppler_process_CS(echo_r_mtx, sigma_n_o, lambda_d);
+    fprintf('   Doppler_processing done.\n');
+    
+    [target_list_RDA_temp, target_map] = rda_detection(stat_RD, P_fa);
+    fprintf('   Target_detection done.\n');
+
+    [target_list_RDA] = angle_estimation(target_list_RDA_temp, echo_mtx, A, d, lambda);
+    fprintf('   Angle_estimation done.\n');
+    target_list_RDA(:, 1) = distance_v(target_list_RDA(:, 1));
+    target_list_RDA(:, 2) = speed_v(target_list_RDA(:, 2));
+
+end
+

0 - example/wzk - RDA/echo_processing_CS_v0.m

@@ -0,0 +1,28 @@
+function [ echo_rd_mtx, target_list_RDA, sigma_n_rd ] = ...
+    echo_processing_CS_v0( echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d )
+
+    % parameters
+    c = 3e8;
+    lambda = c / fc;
+    [num_antenna, N, num_pluse] = size(echo_mtx);
+    
+    
+    % data processing
+    [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pluse);
+    
+    [echo_sFFT_mtx, angle_v] = spatial_FFT(echo_mtx, d, lambda);
+    
+    [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D);
+    
+    [echo_r_mtx, sigma_n_o] = range_process_CS(echo_sFFT_mtx, A, sigma_n, lambda_r);
+    
+    [echo_rd_mtx, stat_RD, sigma_n_rd] = doppler_process_CS(echo_r_mtx, sigma_n_o, lambda_d); %差个方差
+    
+    [target_list_RDA_temp, target_map] = rda_detection(stat_RD, P_fa);
+    target_list_RDA = target_list_RDA_temp;
+    target_list_RDA(:, 1) = distance_v(target_list_RDA(:, 1));
+    target_list_RDA(:, 2) = speed_v(target_list_RDA(:, 2));
+    target_list_RDA(:, 3) = angle_v(target_list_RDA(:, 3));
+
+end
+

0 - example/wzk - RDA/echo_processing_MF.m

@@ -0,0 +1,54 @@
+function [ echo_rd_mtx, target_list_RDA, sigma_n_rd ] =...
+    echo_processing_MF( echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa )
+% This function processes echo data using matching filter method.
+%
+% Usage: 
+% [echo_rd_mtx, target_list_RDA, sigma_n_rd] = echo_processing_MF(echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa)
+% Inputs:
+%   echo_mtx:       Matrix containing the echo data
+%   PRF:            Pulse Repetition Frequency
+%   fs:             Sampling frequency
+%   fc:             Carrier frequency
+%   B:              Bandwidth
+%   D:              Duty ratio
+%   d:              Antenna spacing
+%   sigma_n:        Input noise level of the echo data
+%   P_fa:           False alarm probability threshold
+%
+% Outputs:
+%   echo_rd_mtx:        Echo matrix after range-Doppler processing
+%   target_list_RDA:    List of detected targets
+%   sigma_n_rd:         Estimated noise level after range-Doppler processing
+
+
+    % parameters
+    c = 3e8;
+    lambda = c / fc;
+    [num_antenna, N, num_pluse] = size(echo_mtx);
+    
+    
+    % data processing
+    [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pluse);
+    
+    [echo_sFFT_mtx, angle_v] = spatial_FFT(echo_mtx, d, lambda);
+    fprintf('   Spatial_FFT done.\n');
+    
+    [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D);
+    fprintf('   Generate_chirp_mtx done.\n');
+    
+    [echo_r_mtx, sigma_n_o] = range_process_MF(echo_sFFT_mtx, A, sigma_n);
+    fprintf('   Range_processing done.\n');
+    
+    [echo_rd_mtx, stat_RD, sigma_n_rd] = doppler_process_MF(echo_r_mtx, sigma_n_o);
+    fprintf('   Doppler_processing done.\n');
+    
+    [target_list_RDA_temp, target_map] = rda_detection(stat_RD, P_fa);
+    fprintf('   Target_detection done.\n');
+    
+    [target_list_RDA] = angle_estimation(target_list_RDA_temp, echo_mtx, A, d, lambda);
+    fprintf('   Angle_estimation done.\n');
+    target_list_RDA(:, 1) = distance_v(target_list_RDA(:, 1));
+    target_list_RDA(:, 2) = speed_v(target_list_RDA(:, 2));
+
+end
+

0 - example/wzk - RDA/echo_processing_MF_v0.m

@@ -0,0 +1,28 @@
+function [ echo_rd_mtx, target_list_RDA, sigma_n_rd ] =...
+    echo_processing_MF_v0( echo_mtx, PRF, fs, fc, B, D, d, sigma_n, P_fa )
+
+    % parameters
+    c = 3e8;
+    lambda = c / fc;
+    [num_antenna, N, num_pluse] = size(echo_mtx);
+    
+    
+    % data processing
+    [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, num_pluse);
+    
+    [echo_sFFT_mtx, angle_v] = spatial_FFT(echo_mtx, d, lambda);
+    
+    [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D);
+    
+    [echo_r_mtx, sigma_n_o] = range_process_MF(echo_sFFT_mtx, A, sigma_n);
+    
+    [echo_rd_mtx, stat_RD, sigma_n_rd] = doppler_process_MF(echo_r_mtx, sigma_n_o); %差个方差
+    
+    [target_list_RDA_temp, target_map] = rda_detection(stat_RD, P_fa);
+    target_list_RDA = target_list_RDA_temp;
+    target_list_RDA(:, 1) = distance_v(target_list_RDA(:, 1));
+    target_list_RDA(:, 2) = speed_v(target_list_RDA(:, 2));
+    target_list_RDA(:, 3) = angle_v(target_list_RDA(:, 3));
+
+end
+

0 - example/wzk - RDA/generate_chirp_mtx.m

@@ -0,0 +1,81 @@
+function [ A, signal_t ] = generate_chirp_mtx(PRF, B, fs, D, gamma, sign_mid)
+% Generates a chirp measurement matrix
+%
+% Usage: 
+% [A, signal_t] = generate_chirp_mtx(PRF, B, fs, D, gamma, sign_mid)
+%
+% Inputs:
+%   PRF:         Pulse Repetition Frequency
+%   B:           Bandwidth
+%   fs:          Sampling frequency
+%   D:           Duty ratio
+%   gamma:       Compression ratio
+%   sign_mid:    if sign_mid is 0 means the initial frequency is 0,
+%                if sign_mid is 1 means the initial frequency is -B/2,
+%
+% Outputs:
+%   A:           Chirp measurement matrix
+%   signal_t:    Transmitting signal
+
+
+    % parameters
+    if nargin < 5
+        gamma = 0.5;
+    end
+    if nargin < 6
+        sign_mid = 0;
+    end
+    Tr = 1 / PRF;
+    Tp = Tr * D;
+    K = B / Tp;
+    
+    % generate_signal
+    N = Tr * fs;
+    N_high = Tp * fs;
+    N_mtx = round(N / gamma);
+    
+    signal_t  = zeros(1, N);
+    for i = 1: N_high
+        tp = i * (1 / fs) - sign_mid * N_high / fs / 2;
+        signal_t(1, i) = exp(1j * 2 * pi * 0.5 * K * tp .^ 2);
+    end
+
+    signal_t_2fs  = zeros(1, N_mtx);
+    for i = 1: round(N_high / gamma)
+        tp = (i+1) * (1 / fs / 2) - sign_mid * N_high / fs / 2;
+        signal_t_2fs(1, i) = exp(1j * pi * K * tp .^ 2);
+    end
+    
+    % generate chirp matrix
+    A = generate_matrix_by_signal2fs(transpose(signal_t_2fs));
+
+end
+
+
+%% sub-functions
+% generate chirp matrix when gamma=0.5
+function [ mtx ] = generate_matrix_by_signal2fs(signal)
+
+    mtx = [];
+    l = round(length(signal) / 2);
+    temp1 = signal(1:2:end);
+    temp2 = circshift(signal(2:2:end), 1);
+    if length(temp2) ~= length(temp1)
+        temp2 = [temp2; 0];
+    end
+    
+    for i = 1:l
+        t1 = circshift(temp1, i-1);
+        t2 = circshift(temp2, i-1);
+        if i - 1 > 0
+            t1(1:i - 1,1) = 0;
+        end
+        if i - 1 > 0
+            t2(1:i - 1,1) = 0;
+        end
+        mtx = [mtx,t1,t2];
+    end
+    
+end
+
+

0 - example/wzk - RDA/get_range_speed_val.m

@@ -0,0 +1,35 @@
+function [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, numP, gamma_r, gamma_d)
+% Calculates range and speed values for node index
+%
+% Usage: 
+% [distance_v, speed_v] = get_range_speed_val(PRF, fs, fc, numP, gamma_r, gamma_d)
+%
+% Inputs:
+%   PRF:         Pulse Repetition frequency
+%   fs:          Sampling frequency
+%   fc:          Carrier frequency
+%   numP:        Number of pulses
+%   gamma_r:     Range compression ratio(0.5)
+%   gamma_d:     Doppler compression ratio(0.5)
+%
+% Outputs:
+%   distance_v:  Real distance value
+%   speed_v:     Real speed value
+
+
+    % parameters
+    if nargin < 6
+        gamma_d = 0.5;
+    end
+    if nargin < 5
+        gamma_r = 0.5;
+    end
+    c = 3e8;
+    Tr = 1 / PRF;
+    lambda = c / fc;
+    
+    % calculate
+    distance_v = 0: (c / fs / 2 * gamma_r) : (Tr * c / 2 - c / fs / 2 * gamma_r);
+    speed_v = -(-PRF / 2: PRF / numP * gamma_d: PRF / 2 - PRF / numP * gamma_d) * lambda / 2 ;
+    
+end

0 - example/wzk - RDA/main_CS.m

@@ -0,0 +1,81 @@
+clc
+clear
+close all
+
+%% load echo data
+filename = 'Raw_Echo_5dB';
+load(['./data/', filename, '.mat']);
+
+
+%% parameters
+% ladar
+PRF = 5000;
+B = 5e6;
+D = 0.1;
+Tp = 2e-5;
+fs = 5e6;
+fc = 1.25e9;
+
+% antenna
+num_antenna = 18;
+d = 0.12;
+
+% echo
+num_pulse = 64;
+sigma_n = 0.1;
+
+% detect
+P_fa = 1e-6;
+
+% algorithm
+gamma_r = 0.5;
+gamma_d = 0.5;
+lambda_r = 0.005;
+lambda_d = 0.15;
+
+% accurate angle
+sign_AA = 1;
+
+
+%% data processing
+fprintf('Using CS\n');
+fprintf(['File: ', filename, '\n']);
+
+tic;
+if sign_AA == 0
+    [ echo_rd_mtx, target_list_RDA, sigma_n_out ] = ...
+    echo_processing_CS_v0( Raw_Echo, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d );
+else
+    [ echo_rd_mtx, target_list_RDA, sigma_n_out ] = ...
+    echo_processing_CS( Raw_Echo, PRF, fs, fc, B, D, d, sigma_n, P_fa, lambda_r, lambda_d );
+end
+toc;
+
+
+%% plot
+Fontsize = 18;
+plot_width = 800;
+plot_height = 600;
+Linewidth = 2;
+Markersize = 8;
+
+figure(1)
+scatter3(target_list_RDA(:,1), target_list_RDA(:,2),...
+    target_list_RDA(:,3), 'filled', 'o')
+xlabel('Range(m)');
+ylabel('Speed(m/s)');
+zlabel('Angle(°)');
+zlim([-90, 90]);
+ylim([80, 150]);
+xlim([17000, 19000]);
+set(gca, 'FontSize', Fontsize);
+title('node')
+set(gcf, 'position', [100, 200, plot_width+100, plot_height+50]);
+set(gca,'fontsize',18,'fontname','Times');
+
+
+%% save
+save(['./output/', filename, '_CS_Pfa', num2str(P_fa), '.mat'],...
+    'target_list_RDA', 'echo_rd_mtx', 'sigma_n_out', 'P_fa');
+fprintf(['[', filename, '_CS_Pfa', num2str(P_fa), '.mat]', ' saved.\n\n']);
+

0 - example/wzk - RDA/main_MF.m

@@ -0,0 +1,79 @@
+clc
+clear
+close all
+
+%% load echo data
+filename = 'Raw_Echo_5dB';
+load(['./data/', filename, '.mat']);
+
+
+%% parameters
+% ladar
+PRF = 5000;
+B = 5e6;
+D = 0.1;
+Tp = 2e-5;
+fs = 5e6;
+fc = 1.25e9;
+
+% antenna
+num_antenna = 18;
+d = 0.12;
+
+% echo
+num_pulse = 64;
+sigma_n = 0.1;
+
+% detect
+P_fa = 1e-6;
+
+% algorithm
+gamma_r = 0.5;
+gamma_d = 0.5;
+
+% accurate angle
+sign_AA = 1;
+
+
+%% data processing
+fprintf('Using MF\n');
+fprintf(['File: ', filename, '\n']);
+
+tic;
+if sign_AA == 0
+    [ echo_rd_mtx, target_list_RDA, sigma_n_out ] = ...
+    echo_processing_MF_v0( Raw_Echo, PRF, fs, fc, B, D, d, sigma_n, P_fa );
+else
+    [ echo_rd_mtx, target_list_RDA, sigma_n_out ] = ...
+    echo_processing_MF( Raw_Echo, PRF, fs, fc, B, D, d, sigma_n, P_fa );
+end
+toc;
+
+
+%% plot
+Fontsize = 18;
+plot_width = 800;
+plot_height = 600;
+Linewidth = 2;
+Markersize = 8;
+
+figure(1)
+scatter3(target_list_RDA(:,1), target_list_RDA(:,2),...
+    target_list_RDA(:,3), 'filled', 'o')
+xlabel('Range(m)');
+ylabel('Speed(m/s)');
+zlabel('Angle(°)');
+zlim([-90, 90]);
+ylim([80, 150]);
+xlim([17000, 19000]);
+set(gca, 'FontSize', Fontsize);
+title('node')
+set(gcf, 'position', [100, 200, plot_width+100, plot_height+50]);
+set(gca,'fontsize',18,'fontname','Times');
+
+
+%% save
+save(['./output/', filename, '_MF_Pfa', num2str(P_fa), '.mat'],...
+    'target_list_RDA', 'echo_rd_mtx', 'sigma_n_out', 'P_fa');
+fprintf(['[', filename, '_MF_Pfa', num2str(P_fa), '.mat]', ' saved.\n\n']);
+

0 - example/wzk - RDA/plot_result.m

@@ -0,0 +1,38 @@
+clc
+clear
+close all
+
+
+Fontsize = 18;
+plot_width = 800;
+plot_height = 600;
+Linewidth = 2;
+Markersize = 8;
+
+
+Pfa_set = '1e-06';
+method = 'CS';
+SNR = 5;
+
+filename = ['Raw_Echo_', num2str(SNR), 'dB_', method, '_Pfa', Pfa_set];
+
+load(['./output/', filename, '.mat']);
+
+%%
+figure(1)
+scatter3(target_list_RDA(:,1), target_list_RDA(:,2),...
+    target_list_RDA(:,3), 'filled', 'o')
+xlabel('Range(m)');
+ylabel('Speed(m/s)');
+zlabel('Angle(°)');
+zlim([-90, 90]);
+% ylim([80, 150]);
+% xlim([17000, 19000]);
+set(gca, 'FontSize', Fontsize);
+title('node')
+set(gcf, 'position', [100, 200, plot_width+100, plot_height+50]);
+set(gca,'fontsize',18,'fontname','Times');
+
+
+
+

0 - example/wzk - RDA/range_process_CS.m

@@ -0,0 +1,155 @@
+function [ echo_r_mtx, sigma_n_o ] = range_process_CS(echo_mtx, A, sigma_n, lambda, delta)
+% Range processing for echo data using compressed sensing
+%
+% Usage: 
+% [echo_r_mtx, sigma_n_o] = range_process_CS(echo_mtx, A, sigma_n, lambda, delta)
+%
+% Inputs:
+%   echo_mtx:     Original echo data
+%   A:            Chirp measurement matrix
+%   sigma_n:      Input noise standard deviation
+%   lambda:       LASSO weight
+%   delta:        Convergence normalized difference(2e-7)
+%
+% Outputs:
+%   echo_r_mtx:   Range processing echo data
+%   sigma_n_o:    Output noise standard deviation
+
+
+    % parameters
+    if nargin < 5
+        delta = 2e-7;
+    end
+    [M, N] = size(A);
+    [lenA, M, lenP] = size(echo_mtx);
+    
+    % normalization
+    J1 = A*A';
+    lambda_J=eig(J1);
+    A = A / sqrt(lambda_J(end));
+    echo_mtx = echo_mtx ./ sqrt(lambda_J(end));
+    sigma_n = sigma_n / sqrt(lambda_J(end));
+    
+    % compressed sensing
+    echo_r_mtx = zeros(lenA, N, lenP);
+    sigma_n_o_cnt = zeros(lenA, lenP);
+    for numA = 1: lenA
+        parfor numP = 1: lenP
+            sr = echo_mtx(numA, :, numP);
+            y = transpose(sr);
+            % LASSO
+            x_FISTA = FISTA(y, A, lambda, delta);
+            % debiased LASSO
+            [x_d, sigma_w] = cal_debiased_LASSO(x_FISTA, A, y, lambda, sigma_n);
+            sigma_n_o_cnt(numA, numP) = abs(sigma_w);
+            echo_r_mtx(numA, :, numP) = x_d;
+        end
+    end
+    
+    % calculate output noise
+    sigma_n_o = mean(mean(sigma_n_o_cnt));
+    
+end
+
+
+%% sub-functions
+% algorithm for LASSO
+% y: measurements
+% A: measurement matrix
+% lambda: LASSO weight
+% delta: convergence normalized difference
+% z: LASSO estimator
+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
+
+
+% soft threshold function
+% x: processing object
+% thd: threshold
+% y: result
+function y = sft_thd(x, thd)
+
+    if isequal(size(x), size(thd))
+        tmp = abs(x);
+        y = x;
+        y(tmp <= thd) = 0;
+        y(tmp > thd) = (tmp(tmp > thd) - thd(tmp > thd)) .* x(tmp > thd) ./ tmp(tmp > thd);
+    else
+        tmp = abs(x);
+        y = x;
+        y(tmp <= thd) = 0;
+        y(tmp > thd) = (tmp(tmp > thd) - thd) .* x(tmp > thd) ./ tmp(tmp > thd);
+    end
+
+end
+
+
+% calculate debiased LASSO estimator
+% x: LASSO estimator
+% A: measurement matrix
+% y: measurements
+% lambda: LASSO weight
+% sigma_n: input noise standard deviation
+% x_d: debiased LASSO estimator
+% sigma_d: equivalent noise standard deviation
+function [x_d, sigma_d] = cal_debiased_LASSO(x, A, y, lambda, sigma)
+
+    [M, N] = size(A);
+    gamma = M/N;
+    hat_Q1 = gamma;
+    [~, D] = eig(A'*A);
+    d = diag(D);
+
+    diff = 1;
+    T = 1000;
+    t = 0;
+
+    while (t < T) && (diff > 1e-6)
+        Q1_pre = hat_Q1;
+        rho = mean((2 - lambda./(hat_Q1*abs(x) + lambda)).*(abs(x) > 1e-4))/2;
+        hat_Q1 = rho/mean(1./(d + (1-rho)*hat_Q1/rho));
+        diff = abs(Q1_pre - hat_Q1);
+        t = t+1;
+    end
+
+    x_d = x + 1/hat_Q1*A'*(y - A*x);
+
+    chi = rho/hat_Q1;
+    hat_Q2 = 1/chi - hat_Q1;
+
+    t = -hat_Q2;
+    t_prime = -1/mean((1./(d+hat_Q2)).^2);
+    G_prime = t + 1/chi;
+    G_wprime = t_prime + 1/chi/chi;
+    RSS = sum(abs(y - A*x).^2)/M;
+    hat_chi = gamma*G_wprime/(2*G_prime-2*chi*G_wprime)*RSS +...
+        (-G_wprime*gamma+G_prime*G_prime)/(2*G_prime-2*chi*G_wprime)*sigma^2;
+    sigma_d = sqrt(2*hat_chi)/hat_Q1;
+
+end

0 - example/wzk - RDA/range_process_MF.m

@@ -0,0 +1,36 @@
+function [ echo_r_mtx, sigma_n_o ] = range_process_MF(echo_mtx, A, sigma_n)
+% Range processing for echo data using matching filter
+%
+% Usage: 
+% [echo_r_mtx, sigma_n_o] = range_process_MF(echo_mtx, A, sigma_n)
+%
+% Inputs:
+%   echo_mtx:     Original echo data
+%   A:            Chirp measurement matrix
+%   sigma_n:      Input noise standard deviation
+%
+% Outputs:
+%   echo_r_mtx:   Range processing echo data
+%   sigma_n_o:    Output noise standard deviation
+
+
+    % parameters
+    multiple_r = A(:, 1)' * A(:, 1);
+    [M, N] = size(A);
+    [lenA, M, lenP] = size(echo_mtx);
+    
+    % matched filtering
+    echo_r_mtx = zeros(lenA, N, lenP);
+    for numA = 1: lenA
+        for numP = 1: lenP
+            sr = echo_mtx(numA, :, numP);
+            mf_res = A' * transpose(sr);
+            mf_norm = mf_res ./ multiple_r;
+            echo_r_mtx(numA, :, numP) = mf_norm;
+        end
+    end
+    
+    % calculate output noise
+    sigma_n_o = sqrt(sigma_n^2 / multiple_r);
+    
+end

0 - example/wzk - RDA/rda_detection.m

@@ -0,0 +1,47 @@
+function [ target_list_RDA, target_map ] = rda_detection(stat_RD, P_fa)
+% Target detection, estimating distance, velocity, and angle intervals
+%
+% Usage: 
+% [target_list_RDA, target_map] = rda_detection(stat_RD, P_fa)
+%
+% Inputs:
+%   stat_RD:            Range-Doppler statistics
+%   P_fa:               Probability of false alarm
+%
+% Outputs:
+%   target_list_RDA:    Target list include range, doppler and angle information
+%   target_map:         Show targets in RD map
+
+
+    % parameters
+    [lenA, lenR, lenP] = size(stat_RD);
+    target_map = zeros(lenR, lenP);
+    
+    % detection threshold
+    d_thd = chi2inv(1 - P_fa, 2) / 2;
+    
+    % find target
+    target_list_RD = [];
+    for i = 1: lenA
+        RD_map = squeeze(stat_RD(i, :, :));
+        detect_map = zeros(size(RD_map));
+        detect_map(RD_map > d_thd) = 1;
+%         target_map(i, :, :) = detect_map;
+        [r, c] = find(detect_map);
+        temp_node = [r, c];
+        target_list_RD = [target_list_RD; temp_node];
+    end
+    target_list_RD_new = unique(target_list_RD, 'rows');
+    
+    % estimate target angle interval
+    target_list_angle = [];
+    for i = 1: size(target_list_RD_new, 1)
+        tgt = target_list_RD_new(i, :);
+        stat_vct = stat_RD(:, tgt(1), tgt(2));
+        [v, idx] = max(stat_vct);
+        target_map(tgt(1), tgt(2)) = idx;
+        target_list_angle = [target_list_angle; idx];
+    end
+    target_list_RDA = [target_list_RD_new, target_list_angle];
+    
+end

0 - example/wzk - RDA/spatial_FFT.m

@@ -0,0 +1,31 @@
+function [ echo_sFFT_mtx, angle_index ] = spatial_FFT(echo_mtx, d, lambda)
+% FFT for echo space domain
+%
+% Usage: 
+% [echo_sFFT_mtx, angle_index] = spatial_FFT(echo_mtx, d, lambda)
+%
+% Inputs:
+%   echo_mtx:     Matrix containing the echo data
+%   d:            Antenna spacing
+%   lambda:       Wave length
+%
+% Outputs:
+%   echo_sFFT_mtx:  Processed echo matrix
+%   angle_index:    True value of the angle
+
+
+    % parameters
+    [lenA, lenR, lenP] = size(echo_mtx);
+    
+    % calculate angle
+    angle_index = asin(-1 * (-1/2: 1/lenA: 1/2 - 1/lenA) * lambda / d) / pi * 180;
+    
+    % spatial FFT
+    echo_mtx = reshape(echo_mtx, [lenA, lenR * lenP]);
+    echo_sFFT_mtx = zeros(size(echo_mtx));
+    for i = 1: lenR * lenP
+        echo_sFFT_mtx(:, i) = fftshift(fft(echo_mtx(:, i))) ./ sqrt(lenA);
+    end
+    echo_sFFT_mtx = reshape(echo_sFFT_mtx, [lenA, lenR, lenP]);
+    
+end

0 - example/wzk - main_process/angle_est.m

@@ -0,0 +1,43 @@
+% echo_r_mtx: range processing echo data
+% target_list_RD: target list include range and doppler information
+% d: antenna spacing
+% fc: carrier frequency
+% acc: angle estimation accuracy
+
+% target_list_RDA: target list include range, doppler and angle information
+% RA_map: range-angle map
+% THETA: scale of angle
+
+
+function [ target_list_RDA, RA_map, THETA ] = angle_est(echo_r_mtx, target_list_RD, d, fc, acc)
+
+    % parameters
+    if nargin < 5
+        acc = 181;
+    end
+    [lenA, lenR, lenP] = size(echo_r_mtx);
+    c = 3e8;      
+    lambda = c / fc;
+    
+    % initialization
+    echo_mtx = squeeze(echo_r_mtx(:, :, 1));
+    
+    % CBF
+    THETA = linspace(-90, 90, round(acc));
+    RA_map = zeros(length(THETA), lenR);
+    for i = 1: length(THETA)
+        a = exp((0: lenA - 1)' * -1j * 2 * pi / lambda * d * sin(THETA(i) / 180 * pi)); 
+        RA_map(i, :) = a'* echo_mtx;  
+    end
+    RA_map = transpose(RA_map);
+    
+    % calculate angle
+    range_angle = zeros(1, lenR);
+    for j = 1: lenR
+        [v, num] = max(abs(RA_map(j, :)));
+        range_angle(j) = THETA(num);
+    end
+    targets_angle = range_angle(target_list_RD(:, 2));
+    target_list_RDA = [target_list_RD, transpose(targets_angle)];
+
+end

0 - example/wzk - main_process/doppler_process_CS.m

@@ -0,0 +1,187 @@
+% echo_r_mtx: range processing echo data
+% sigma_n: input noise standard deviation
+% lambda: LASSO weight
+% gamma: compressed ratio
+% delta: convergence normalized difference
+
+% echo_rd_mtx: range-doppler processing echo data
+% stat_RD: range-doppler statistics
+
+
+function [ echo_rd_mtx, stat_RD ] = doppler_process_CS(echo_r_mtx, sigma_n, lambda, gamma, delta)
+
+    % parameters
+    if nargin < 5
+        delta = 1e-6;
+    end
+    if nargin < 4
+        gamma = 0.5;
+    end
+    [lenA, lenR, M] = size(echo_r_mtx);
+    N = round(M / gamma);
+    iter_max_VAMP = 1000;  
+    lambda_v = zeros(N, 1) + lambda;
+    
+    % generate mtx
+    F_ori = dftmtx(N);
+    F = F_ori(1:M,:);
+    F_inv = conj(F) / N;
+    
+    % normalization
+    A = (sqrt(N) * eye(M)) * F_inv;
+    echo_r_mtx = sqrt(N) .* echo_r_mtx;
+    sigma_n = sqrt(N) * sigma_n;
+    
+    % doppler matched filtering
+    echo_rd_mtx = zeros(lenA, lenR, N);
+    stat_RD = zeros(lenA, lenR, N);
+    for numA = 1: lenA
+        for numR = 1: lenR
+            sample = squeeze(echo_r_mtx(numA, numR, :));
+            y = sample;
+            x_LASSO = cVAMPro(y, A, lambda_v, delta, iter_max_VAMP);
+            [x_d_CROD, sigma_CROD] = CROD(y, A, x_LASSO, lambda, sigma_n);
+            stat_RD(numA, numR, :) = abs(fftshift(x_d_CROD) / sigma_CROD).^2;
+            echo_rd_mtx(numA, numR, :) = fftshift(x_d_CROD);
+        end
+    end
+    
+end
+
+% algorithm for LASSO
+% y: measurements
+% A: measurement matrix
+% lambda: LASSO weight
+% tau: convergence normalized difference
+% Kit: maximum number of iterations
+% LASSO estimator
+function x_hat_wl = 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;
+        % next
+        h_1 = h_1_next;
+        Q_1 = Q_1_next;
+    end
+
+end
+
+
+% soft threshold function
+% x: processing object
+% thd: threshold
+% y: result
+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
+
+% calculate debiased LASSO estimator
+% y: measurements
+% A: measurement matrix
+% x_LASSO: LASSO estimator
+% lambda: LASSO weight
+% sigma_n: input noise standard deviation
+% x_d_CROD: debiased LASSO estimator
+% sigma_CROD: equivalent noise standard deviation estimator
+function [ x_d_CROD, sigma_CROD ] = CROD(y, A, x_LASSO, lambda, sigma_n)
+
+    [m, n] = size(A);
+    gamma = m / n;
+
+    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;
+
+end

0 - example/wzk - main_process/doppler_process_MF.m

@@ -0,0 +1,38 @@
+% echo_r_mtx: range processing echo data
+% sigma_n: input noise standard deviation
+% gamma: compressed ratio
+
+% echo_rd_mtx: range-doppler processing echo data
+% stat_RD: range-doppler statistics
+
+
+function [ echo_rd_mtx, stat_RD ] = doppler_process_MF(echo_r_mtx, sigma_n, gamma)
+
+    % parameters
+    if nargin < 3
+        gamma = 0.5;
+    end
+    [lenA, lenR, M] = size(echo_r_mtx);
+    N = round(M / gamma);
+    
+    % generate mtx
+    F_ori = dftmtx(N);
+    F = F_ori(1:M,:);
+    multiple_d = F(:,1)' * F(:,1);
+    
+    % doppler matched filtering
+    echo_rd_mtx = zeros(lenA, lenR, N);
+    for numA = 1: lenA
+        for numR = 1: lenR
+            sample = squeeze(echo_r_mtx(numA, numR, :));
+            dpl_temp = transpose(F) * sample;
+            dpl_norm = fftshift(dpl_temp ./ multiple_d);
+            echo_rd_mtx(numA, numR, :) = dpl_norm;
+        end
+    end
+    
+    % calculate output noise
+    sigma_n_o = sqrt(sigma_n^2 / multiple_d);
+    stat_RD = abs(echo_rd_mtx ./ sigma_n_o).^2;
+    
+end

0 - example/wzk - main_process/draw_target_map.m

@@ -0,0 +1,21 @@
+% target_list_RDA: target list, including range, doppler and angle information
+% stat_RD: range-doppler statistics
+
+% target_map: show targets information
+
+
+function [ target_map ] = draw_target_map(target_list_RDA, stat_RD)
+
+    % initialization
+    target_map = ones(size(stat_RD)) .* -300;
+    t_map = target_list_RDA(:, 1);
+    t_R = target_list_RDA(:, 2);
+    t_D = target_list_RDA(:, 3);
+    t_A = target_list_RDA(:, 4);
+    
+    % draw_map
+    for i = 1: size(target_list_RDA, 1)
+        target_map(t_map(i), t_R(i), t_D(i)) = t_A(i);
+    end
+
+end

0 - example/wzk - main_process/generate_chirp_mtx.m

@@ -0,0 +1,72 @@
+% PRF: pulse repetition frequency
+% B: bandwidth
+% fs: sampling rate
+% D: duty ratio
+% gamma: compression ratio
+
+% A: chirp matrix
+% signal_t: transmitting beam
+
+
+function [ A, signal_t ] = generate_chirp_mtx(PRF, B, fs, D, gamma, sign_mid)
+
+    % parameters
+    if nargin < 5
+        gamma = 0.5;
+    end
+    if nargin < 6
+        sign_mid = 0;
+    end
+    Tr = 1 / PRF;
+    Tp = Tr * D;
+    K = B / Tp;
+    
+    % generate_signal
+    N = Tr * fs;
+    N_high = Tp * fs;
+    N_mtx = round(N / gamma);
+    
+    signal_t  = zeros(1, N);
+    for i = 1: N_high
+        tp = i * (1 / fs) - sign_mid * N_high / fs / 2;
+        signal_t(1, i) = exp(1j * 2 * pi * 0.5 * K * tp .^ 2);
+    end
+
+    signal_t_2fs  = zeros(1, N_mtx);
+    for i = 1: round(N_high / gamma)
+        tp = (i+1) * (1 / fs / 2) - sign_mid * N_high / fs / 2;
+        signal_t_2fs(1, i) = exp(1j * pi * K * tp .^ 2);
+    end
+    
+    % generate chirp matrix
+    A = generate_matrix_by_signal2fs(transpose(signal_t_2fs));
+
+end
+
+
+% generate chirp matrix when gamma=0.5
+function [ mtx ] = generate_matrix_by_signal2fs(signal)
+
+    mtx = [];
+    l = round(length(signal) / 2);
+    temp1 = signal(1:2:end);
+    temp2 = circshift(signal(2:2:end), 1);
+    if length(temp2) ~= length(temp1)
+        temp2 = [temp2; 0];
+    end
+    
+    for i = 1:l
+        t1 = circshift(temp1, i-1);
+        t2 = circshift(temp2, i-1);
+        if i - 1 > 0
+            t1(1:i - 1,1) = 0;
+        end
+        if i - 1 > 0
+            t2(1:i - 1,1) = 0;
+        end
+        mtx = [mtx,t1,t2];
+    end
+    
+end
+
+