基于 压缩感知的尽头: 原子范数最小化 中的原子范数最小化算法， 笔者做了一些matlab的仿真， 作为简单的例程，希望帮助大家进一步理解算法和自定义的拓展。
由于凸问题的求解需要使用 CVX， 因此需要读者先自行安装好 matlab 的 CVX包。
1-D 无噪场景假设接收天线有 64 64 64根， 有 3 3 3个单天线的目标源同时发射信号， 接收信号可以表示为： z = [ a ( f 1 ) , a ( f 2 ) , a ( f 3 ) ] [ s 1 s 2 s 3 ] T = A ( f ) s z=[a(f_1), a(f_2), a(f_3)][s_1\; s_2\; s_3]^T=A(f)s z=[a(f1),a(f2),a(f3)][s1s2s3]T=A(f)s 其中， a ( f ) = [ 1 , e i 2 π f , … , e i 2 π ( M − 1 ) f ] T \boldsymbol{a}(f)=\left[1, e^{i 2 \pi f}, \ldots, e^{i 2 \pi(M-1) f}\right]^{T} a(f)=[1,ei2πf,…,ei2π(M−1)f]T 代表天线响应矢量。 这里假设没有噪声。 那么 从 z z z 中估计 f f f， 我们使用原子范数最小化算法。 代码如下：
N = 64; % number of antennas
f = [0.1, 0.4, 0.3]; % f = cos(theta)
A = exp(1j * 2 * pi * (0 : N-1)' * f); % A = [a(f1), ..., a(fK)]
s = ones(3, 1); % 3 sources
z = A * s;
cvx_begin sdp
cvx_solver sdpt3
variable T(N, N) hermitian toeplitz
variable x
minimize(0.5 * x + 0.5 * T(1,1))
[ x z'; z T] >= 0;
cvx_end
[Phi,P] = rootmusic(T, 3, 'corr');
Phi / 2 / pi
中间 cvx_begin 和 cvx_end 包围的部分， 我们就是在求解 和 原子范数最小化等价的 SDP问题。 CVX可以返回解得的 T。 最后， 我们使用 rootmusic函数， 可以将
T
T
T 对应的
f
f
f 得到， 即为 Phi。 （本来是做范德蒙德分解，然后得到， 但通过rootmusic提取出
f
f
f 更加便捷。）
输出结果如下：
>> ANM_1D
Calling SDPT3 4.0: 4225 variables, 128 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------
num. of constraints = 128
dim. of sdp var = 130, num. of sdp blk = 1
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
HKM 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
-------------------------------------------------------------------
0|0.000|0.000|4.3e+03|1.1e+01|2.2e+05| 0.000000e+00 0.000000e+00| 0:0:00| chol 1 1
1|0.893|1.000|4.6e+02|2.3e-01|2.9e+04|-5.933532e+00 -3.889838e+01| 0:0:00| chol 1 1
2|0.988|1.000|5.4e+00|2.3e-02|3.8e+02|-2.814220e+00 -3.796329e+01| 0:0:00| chol 1 1
3|0.943|1.000|3.1e-01|2.3e-03|3.3e+01|-3.966408e+00 -1.985735e+01| 0:0:00| chol 1 1
4|0.888|0.790|3.5e-02|6.6e-04|1.3e+01|-2.602549e+00 -1.421772e+01| 0:0:00| chol 1 1
5|1.000|1.000|1.2e-09|7.0e-03|2.8e+00|-2.914921e+00 -5.690264e+00| 0:0:00| chol 1 1
6|0.980|0.987|1.4e-10|9.2e-05|3.6e-02|-2.998482e+00 -3.034630e+00| 0:0:00| chol 1 1
7|0.979|0.987|4.9e-11|1.4e-06|5.0e-04|-2.999973e+00 -3.000466e+00| 0:0:00| chol 1 1
8|0.978|0.986|9.6e-11|2.1e-08|7.3e-06|-2.999999e+00 -3.000007e+00| 0:0:00| chol 1 1
9|0.976|0.990|2.1e-10|2.2e-10|2.1e-07|-3.000000e+00 -3.000000e+00| 0:0:00| chol 1 1
10|0.950|0.988|1.3e-11|2.5e-11|9.0e-09|-3.000000e+00 -3.000000e+00| 0:0:01|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations = 10
primal objective value = -3.00000000e+00
dual objective value = -3.00000001e+00
gap := trace(XZ) = 8.95e-09
relative gap = 1.28e-09
actual relative gap = 1.31e-09
rel. primal infeas (scaled problem) = 1.26e-11
rel. dual " " " = 2.48e-11
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual " " " = 0.00e+00
norm(X), norm(y), norm(Z) = 7.4e-01, 1.4e+01, 8.0e+01
norm(A), norm(b), norm(C) = 6.5e+01, 1.7e+00, 1.5e+01
Total CPU time (secs) = 0.50
CPU time per iteration = 0.05
termination code = 0
DIMACS: 1.4e-11 0.0e+00 1.5e-10 0.0e+00 1.3e-09 1.3e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +3
ans =
0.4000
0.1000
0.3000
中间是CVX的求解过程。 可以通过在 cvx_begin 后面 输入 quiet参数使其不输出。 最后的结果就是根据原子范数最小化求得的
f
f
f 的值， 可以看到。 分毫不差。
把 z= A * s替换为 z = A * s + 0.5 * (randn(N, 1) + 1j * randn(N, 1));， 可以人为的加入噪声，当噪声大小增大时， 可以看到估计结果会出现一些误差。
s=ones(3,1)可以替换为 s=randn(3,1)之类的，加入随机性。
能看到这里， 说明已经对原子范数有了一定的了解。 那这里就不再多废话， 直接上代码：
N = 64; % number of antennas
L = 5;
f = [0.1, 0.4, 0.3]; % f = cos(theta)
A = exp(1j * 2 * pi * (0 : N-1)' * f); % A = [a(f1), ..., a(fK)]
S = randn(3, L) + 1j * randn(3, L); % 3 sources
Z = A * s ;
cvx_begin sdp quiet
cvx_solver sdpt3
variable T(N, N) hermitian toeplitz
variable X(L, L) hermitian
minimize(trace(X) + trace(T))
[X Z'; Z T] >= 0;
cvx_end
[Phi,P] = rootmusic(T, 3, 'corr');
Phi / 2 / pi
N = 64; % number of antennas
L = 5;
sigma = 0.5;
f = [0.1, 0.4, 0.3]; % f = cos(theta)
A = exp(1j * 2 * pi * (0 : N-1)' * f); % A = [a(f1), ..., a(fK)]
S = randn(3, L) + 1j * randn(3, L); % 3 sources
Y = A * s + sigma * (randn(N, 1) + 1j * randn(N, 1));
lambda = sqrt(N * (L + log(N) + sqrt( 2 * L * log(N)))*sigma);
cvx_begin sdp quiet
cvx_solver sdpt3
variable T(N, N) hermitian toeplitz
variable X(L, L) hermitian
variable Z(N, L) complex
minimize(lambda * (trace(X) + trace(T)) + 0.5*sum_square_abs(vec(Y - Z)))
[X Z'; Z T] >= 0;
cvx_end
[Phi,P] = rootmusic(T, 3, 'corr');
Phi / 2 / pi
