写点题外话,最近,一直在研究Malioutouv的
这是他毕业论文摘出来的一篇论文,他的毕业论文下载地址:D. M. Malioutov. (2003, Jul.) A Sparse Signal Reconstruction Perspective for Source Localization With Sensor Arrays. M.S. thesis, Mass. Inst. Technol., Cambridge, MA. [Online]. Available: http://ssg.mit.edu/~dmm/publications/malioutov_MS_thesis.pdf
虽然毕业论文还是有一些错误,但是对这篇论文的解释比较详尽,有看不懂的地方可以回本科论文里面找,基本都能解释清楚。
其实我是想说,人和人的天资是天差地别的,最近对压缩感知的学习,让我深深地为自己的理解能力和数学基础感到忧虑。上一篇LP球,我是占个位置,因为我还不是很清楚LP与稀疏性的关系,留着以后想清楚了再补充吧。
//正文
论文中的方法说明
示意图:
步骤:
1.获得观测值(与MUSIC算法得到Y=AS+N类似)
2.对观测值进行奇异值分解(SVD())
3.对观测值和稀疏信号进行降维
观测值(阵元数*快拍数)->(阵元数*信元数)
稀疏信号(网格数(可能存在的位置数,如果你认为从1:180都有可能存在,那他的网格数就是181,对应文中的N_θ)*快拍数)->(网格数*信元数)
由于信元数远远小于快拍数(在多快拍的情况下),所以大大减小了这两个矩阵的维度
【这里想快拍数到底什么意思想了很久,现在认为快拍就是对信号进行若干次采样,采样次数多了,得到的结果比较准确,就好比检测产品合格率,你不能拿到一个好的产品就认为这一批产品全是好的,也不能拿到一个坏的就认为这批产品全坏了,因此要选取大量的样本,才能保证结果的准确性。这里暂时理解到这里,回头跟我导讨论一下再来补充】
4.计算Ssv每一行的L2范数
这块卡了很久,没看明白‘spatial index i’到底什么意思?空间索引?很方,当时觉得怎么出现了空间这个问题,后来看了几遍才看明白,i表示的是
N_θ,也就是可能出现点的位置。也就是说如果我们是从0:180°以1°为步长画网格,那我们的i就是(1:1:181)的意思,就是位置,他要是用location表达,可能我就不会看这么费劲,还是要好好学英语啊
5.用约束条件进行最优化
其中约束条件用框图中的直接写是解不出来的,我试了。。。(好尴尬)
要用附录中的说明方法解
其中Ysv就是我们降维后的矩阵(阵元数*信源数),已知
"1"代表(网格数*1)的全1向量
p,q,r,Ssv是变量
pqr都是网格数的实变量
Ssv是(网格数*信源数)的复变量
lamda的选值很复杂,严格来讲应该一次一次试,通过迭代得到一个结果最稳定最适合这个模型的值,当然我是随便写的,2
贴代码:
%2017.6.6 %author:Lola %功能:实现稀疏矩阵l1-SVDDOA估计 %reference:A Sparse Signal Reconstruction Perspective for Source %Localization With Sensor Arrays [Malioutov] clc clear all close all M = 8; %阵元数 K = 2; %信源数 L= 100; %快拍数 d_lamda =0.5; %阵元间距半波长 fc=[2*10^3 ,5*10^3]; %信号频率 fs=[20*10^3]; %采样频率 theta1 = [0 15]; %信号来向 snr=20; %信噪比 for k=1:K s=sqrt(10.^(snr/10))*[randn(K,L)+1j*2*pi*fc'*[1:L]/fs]; %信号(信源数*快拍数) for kk=1:M A(kk,k)=exp(-1j*2*pi*(kk-1)*d_lamda*sin(theta1(k)*pi/180)); %阵列流型(阵元数*信源数) end end X=A*s; Y=awgn(X,snr); %加入高白噪声,得到观测信号 Searching_doa=-90:90; for m=1:length(Searching_doa) AA(:,m) = exp(-1j*(0:M-1)*d_lamda*2*pi*sin(Searching_doa(m)*pi/180)); %构造一度为网格间距的完备基矩阵 end Y=X; %得到观测数据矩阵 Dk1=eye(K); Dk2 = zeros(K,L-K); Dk = [Dk1,Dk2].'; [U,D,V] = svd(Y); %进行奇异值分解 Ysv = Y*V'*Dk; Sumvector=ones(length(Searching_doa),1); cvx_begin quiet variables p q variables r(length(Searching_doa)) variable SSV1(length(Searching_doa),K) complex; expression xsv(length(Searching_doa),1) expressions Zk(M,K) complex minimize(p+2*q); subject to Zk = Ysv-AA*SSV1; Zvec = vec(Zk); %把矩阵转化为向量 norm(Zvec) <=p; %第一个不等式约束 Sumvector'*r<=q; %第二个不等式约束 for i=1:length(Searching_doa) %第三个不等式约束 xsv(i,:)=norm(SSV1(i,:)); end for i=1:length(Searching_doa) xsv(i)<=r(i); end cvx_end power=10*log10(abs(SSV1(:,1))/max(abs(SSV1(:,1)))); plot(Searching_doa,power,'r'); xlabel('DOA/degree'); ylabel('PowerdB');
继续好好学数学
论文不给附录我肯定是不知道那个问题要用二阶锥解得,就算告诉我用二阶锥方法解我也不知道怎么解。这个仿真从研究到写到调试用了差不多五个小时。。。理解能力太差,受到了深深的打击。。。