ASSP算法实现

fdd566

1. function[DD] = CS_ASSP2(Y,P,pth)
   
2. load data1;
   
3. %step1
   
4.     s = 1;
   
5.     k = 1;
   
6.     Is_k1 = [];
   
7.     R_k1 = Y;
   
8.     R_s = +inf;
   
9.     D_k1 = 0;
   
10.     DD = zeros(M*L,R);
    
11. %step2
    
12.     while 1
    
13.          %(1)correlation
    
14.          Z = P'*R_k1;
    
15.          %(2)support estimate
    
16.          for i = 1:L
    
17.              F_Z(i) = norm(Z((i-1)*M+1:i*M,:));
    
18.          end
    
19.          [val,pos] = sort(F_Z,'descend');
    
20.          Is_bpk = union(Is_k1,pos(1:s));
    
21.          %(3)support pruning
    
22.          %D = zeros(M*L,R);
    
23.          P_Is_bpk = [];
    
24.          for i = 1:size(Is_bpk,2)
    
25.             P_Is_bpk(:,(i-1)*M+1:i*M) = P(:,(Is_bpk(i)-1)*M+1:Is_bpk(i)*M);
    
26.          end
    
27.          D_Is_bpk = pinv(P_Is_bpk)*Y;%theta_Is
    
28.          F_D = [];
    
29.          for i = 1:size(Is_bpk,2)
    
30.              F_D(i) = norm(D_Is_bpk((i-1)*M+1:i*M,:));
    
31.          end
    
32.          [val,pos] = sort(F_D,'descend');
    
33.          Is_bk = Is_bpk(pos(1:s));
    
34.          %(4)matrix estimate
    
35.          P_Is_bk = [];
    
36.          for i = 1:size(Is_bk,2)
    
37.             P_Is_bk(:,(i-1)*M+1:i*M) = P(:,(Is_bk(i)-1)*M+1:Is_bk(i)*M);
    
38.          end
    
39.          D_Is_bk = pinv(P_Is_bk)*Y;
    
40.          %(5)residue update
    
41.          R_k = Y - P_Is_bk*D_Is_bk;
    
42.          %(6)matrix update
    
43.         if norm(R_k1) > norm(R_k)
    
44.             %(7)Iteration with fixed sparsity level
    
45.             Is_k1 = Is_bk;
    
46.             R_k1 = R_k;
    
47.             D_k1 = D_Is_bk;
    
48.             k = k+1;
    
49.         else
    
50.             %(8)update sparsity level
    
51.             D_s = D_k1;
    
52.             R_s = R_k1;
    
53.             Is_s = Is_k1;
    
54.             s = s+1;
    
55.         end
    
56.         %condition 1
    
57.         if (norm(R_k) > norm(R_s)) && (norm(R_s) < 1)
    
58.             fprintf('condition 1\n');
    
59.             break;
    
60.         end
    
61.         %condition 1 end
    
62.         %condition 2
    
63.         F_D_s = [];
    
64.         for i = 1:size(Is_bk,2)
    
65.             F_D_s(i) = norm(D_Is_bk((i-1)*M+1:i*M,:));
    
66.         end
    
67.         [val,pos] = sort(F_D_s,'ascend');
    
68.         if val(1) <= sqrt(M*R)*pth
    
69.             fprintf('condition 2\n');
    
70.             break;
    
71.         end
    
72.         %condition 2 end
    
73.     end
    
74.     for i = 1:size(Is_s,2)
    
75.         DD((Is_s(i)-1)*M+1:Is_s(i)*M,:) = D_s((i-1)*M+1:i*M,:);
    
76.     end
    
77. end
    
78. 
    
79. 
    
80. 
    

复制代码

利用ofdm符合间的H具有相同的支撑集，以及不同天线具有共同支撑集实现的sp算法的改进 下面是测试程序。

1. clear all;close all;clc;
   
2. M = 32;%观测值个数（长度）
   
3. N = 4096;
   
4. K = 6;%稀疏度
   
5. L = 64;
   
6. Np = round(N*0.17);
   
7. R = 1;
   
8. %Np = round(Eita_p*N)
   
9. save data1 M N K L N R;
   
10. SNR_val = 30;
    
11. pth = 0.1;
    
12. for i = 1:16
    
13.     MSE_singel_val = [];
    
14.     Eita_p_singel_val = [];
    
15.     for j = 1:20
    
16.         D = CS_D();
    
17.         [MSE_singel_val(j),Eita_p_singel_val(j)] = MSE_Eita(D,Np,SNR_val,pth);
    
18.     end
    
19.     MSE_val(i) = sum(MSE_singel_val)/20
    
20.     Eita_p_val(i) = Eita_p_singel_val(1);
    
21.     Np = Np + floor((N*0.21 - N*0.17)/16);
    
22. end
    
23. 
    
24. SNR_val = 20;
    
25. pth = 0.05;
    
26. Np =  round(N*0.17);
    
27. for i = 1:16
    
28.     MSE_singel_val = [];
    
29.     Eita_p_singel_val = [];
    
30.     for j = 1:20
    
31.         D = CS_D();
    
32.         [MSE_singel_val(j),Eita_p_singel_val(j)] = MSE_Eita(D,Np,SNR_val,pth);
    
33.     end
    
34.     MSE_val_2(i) = sum(MSE_singel_val)/20
    
35.     Eita_p_val_2(i) = Eita_p_singel_val(1);
    
36.     Np = Np + floor((N*0.21 - N*0.17)/16);
    
37. end
    
38. 
    
39. SNR_val = 10;
    
40. pth = 0.05;
    
41. Np =  round(N*0.17);
    
42. for i = 1:16
    
43.     MSE_singel_val = [];
    
44.     Eita_p_singel_val = [];
    
45.     for j = 1:20
    
46.         D = CS_D();
    
47.         [MSE_singel_val(j),Eita_p_singel_val(j)] = MSE_Eita(D,Np,SNR_val,pth);
    
48.     end
    
49.     MSE_val_3(i) = sum(MSE_singel_val)/20
    
50.     Eita_p_val_3(i) = Eita_p_singel_val(1);
    
51.     Np = Np + floor((N*0.21 - N*0.17)/16);
    
52. end
    
53. 
    
54. 
    
55. figure;
    
56. semilogy(Eita_p_val,MSE_val,'-ro');%绘制
    
57. axis([0.17 0.20 1.0e-04 1.0e01]);
    
58. grid on;
    
59. hold on;
    
60. semilogy(Eita_p_val_2,MSE_val_2,'-b+');%绘制
    
61. hold on;
    
62. semilogy(Eita_p_val_3,MSE_val_3,'-k*');%绘制
    
63. hold off;
    
64. legend('SNR_30','SNR_20','SNR_10');

复制代码