script_demo.m 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155
  1. % This matlab script provides a demonstration of the informed and
  2. % constrained NMF algorithms proposed in the article:
  3. % [1] C. DORFFER, M. Puigt, G. Delmaire and G. Roussel, "Informed Nonnegative
  4. % Matrix Factorization Methods for Mobile Sensor Network Calibration", in
  5. % IEEE Transactions on Signal and Information Processing over Networks,
  6. % vol. 4, no. 4, pp. 667-682, Dec. 2018.
  7. %
  8. % Author: Clément DORFFER
  9. % Date: 26/11/2018
  10. % @: clement.dorffer@ensta-bretagne.fr
  11. %%
  12. clear
  13. close all
  14. clc
  15. addpath('functions/')
  16. rng(3)
  17. %% Simulation parameters
  18. N_Ref = 3; % Nb. of reference measurements
  19. N_Cpt = 25; % Nb. of mobile sensors
  20. Mu_beta = .9; % Mean sensors gain
  21. Mu_alpha = 5; % Mean sensors offset
  22. Bound_beta = [.01;1.5]; % Gain boundaries
  23. Bound_alpha = [3.5;6.5]; % Offset boundaries
  24. MV = .9; % Missing Value prop.
  25. RV = 0; % RendezVous prop.
  26. var_n = 0; % Noise variance
  27. %% Scene simulation
  28. [xx,yy] = meshgrid((-1:2/9:1),(-1:2/9:1));
  29. sig = [1;1];
  30. y=@(xx,yy) exp(-[xx;yy]'*diag(sig)*[xx;yy]);
  31. for i = 1:100
  32. g(i) = y(xx(i),yy(i));
  33. end
  34. g = g-min(g);
  35. g = .5*(g/max(g))+1e-5;
  36. G_theo = [ones(100,1),g']; % Theoretical matrix G (see eq.(3) of [1])
  37. figure
  38. subplot(221)
  39. imagesc(reshape(g,10,10))
  40. %% Sensors simulation
  41. F_theo = [max(Bound_beta(1),min(Bound_beta(2),Mu_beta+.5*randn(1,N_Cpt)));...
  42. max(Bound_alpha(1),min(Bound_alpha(2),Mu_alpha+.5*randn(1,N_Cpt)))];
  43. F_theo = [F_theo,[0;1]]; % Theoretical matrix F (see eq.(3) of [1])
  44. %% Data simulation
  45. X_theo = G_theo*F_theo; % Theoretical matrix X (see eq.(2) of [1])
  46. W = zeros(100,N_Cpt+1);
  47. idx_Ref = randperm(100);
  48. idx_Ref = idx_Ref(1:N_Ref); % Reference measurement locations
  49. W(idx_Ref,end) = 1;
  50. N_RV = round(N_Cpt*RV); % Nb. of sensors having a RendezVous
  51. idx_CptRV = randperm(N_Cpt);
  52. idx_CptRV = idx_CptRV(1:N_RV); % Selection of sensors having a RendezVous
  53. idx_RefRV = randi(N_Ref,1,N_Cpt);
  54. idx_RefRV = idx_Ref(idx_RefRV(1:N_RV)); % Selection of the references for each RendezVous
  55. for i = 1 : N_RV
  56. W(idx_RefRV(i),idx_CptRV(i)) = 1;
  57. end
  58. N_data = round((1-MV)*(N_Cpt)*(100-N_Ref)); % Nb. of measurements in data matrix X
  59. xCpt = 1 : 100;
  60. xCpt(idx_Ref) = []; % Reference free locations
  61. [xx,yy] = meshgrid(xCpt,1:N_Cpt); % Possibly sensed locations
  62. idx_data = randperm((100-N_Ref)*N_Cpt);
  63. for i = 1 : N_data
  64. W(xx(idx_data(i)),yy(idx_data(i))) = 1; % Sensor measurement placement
  65. end
  66. N = var_n*randn(100,N_Cpt+1); % Noise simulation
  67. N(:,end) = 0;
  68. N = max(N,-X_theo);
  69. SNR = snr(X_theo(W~=0),N(W~=0));
  70. X = W.*(X_theo+N); % Data matrix X
  71. subplot(222)
  72. imagesc(X)
  73. %% Calibration parameters
  74. % % Common parameters
  75. N_iter = 1.e5; % Maximum nb. of iterations
  76. Omega_G = [ones(100,1),W(:,end)]; % Mask on known values in G (see eq.(14) of [1])
  77. Omega_F = [zeros(2,N_Cpt),[1;1]]; % Mask on known values in F (see eq.(15) of [1])
  78. Phi_G = [ones(100,1),X(:,end)]; % Known values in G (see eq.(14) of [1])
  79. Phi_F = [zeros(2,N_Cpt),[0;1]]; % Known values in F (see eq.(15) of [1])
  80. Ginit = abs(randn(100,2)); % Initial matrix G
  81. Ginit = (1-Omega_G).*Ginit+Phi_G;
  82. Finit = [max(Bound_beta(1),min(Bound_beta(2),Mu_beta+.5*randn(1,N_Cpt)));...
  83. max(Bound_alpha(1),min(Bound_alpha(2),Mu_alpha+.5*randn(1,N_Cpt)))];
  84. Finit = [Finit,[0;1]]; % Initial matrix F
  85. Finit = (1-Omega_F).*Finit+Phi_F;
  86. % % Parameters for the "average constrained" approach (called ACIN-Cal in [1])
  87. Mean_F = mean(F_theo(:,1:N_Cpt),2);
  88. mu = 10; % Regularization weight
  89. % % Parameters for the Sparsity based regularization (called SpIN-Cal in [1])
  90. % Dictionary construction
  91. D = real(ifft(diag(ones(100,1))));
  92. D = [D(:,1:15),g',D(:,1:15),g'];
  93. D(51:end,1:16) = 0;
  94. D(1:50,17:end) = 0;
  95. D = D*diag(1./sqrt(sum(D.^2))); % Atoms normalisation
  96. k = 2; % Nb. of atoms to be choosen by the OMP
  97. lambda = 10; % Regularization weight
  98. %% Calibration
  99. % % IN-Cal
  100. [ G_IN_Cal , F_IN_Cal , RMSE_IN_Cal ] = IN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , F_theo , N_iter );
  101. % % ACIN-Cal
  102. [ G_ACIN_Cal , F_ACIN_Cal , RMSE_ACIN_Cal ] = ACIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , Mean_F , mu , F_theo , N_iter );
  103. % % SpIN-Cal
  104. [ G_SpIN_Cal , F_SpIN_Cal , RMSE_SpIN_Cal ] = SpIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , F_theo , N_iter );
  105. % % SpAIN-Cal
  106. [ G_SpAIN_Cal , F_SpAIN_Cal , RMSE_SpAIN_Cal ] = SpAIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , Mean_F , mu , F_theo , N_iter );
  107. subplot(223)
  108. semilogy(RMSE_IN_Cal)
  109. axis([0 N_iter 1.e-16 1])
  110. hold all
  111. semilogy(RMSE_ACIN_Cal)
  112. semilogy(RMSE_SpIN_Cal)
  113. semilogy(RMSE_SpAIN_Cal)