等距网格细化的单快拍稀疏矩阵的DOA估计

方法：

这里采用的是等距划分。

想了很久也试了很多，最后是二分法给我的启发

二分法的方法参照：

function [i,x,fx]=erfenfa_new_use1(f,bot,top,err)  
%输入：f为要求解的方程函数表达式，bot为求解区间的下界，top为求解区间的上界，err为所要求的误差范围  
%输出：i为二分法求解的次数，x为最终的根，fx为方程的根x对应的函数值（精度要求内接近于零或等于零）  
if f(bot)*f(top)>0  
    disp('\n注意：f(bot)*f(top)>0,无法继续运算，请重新调整区间端点bot和top.\n');  
    return  
end  
%wucha=abs(top_new-bot_new)/2;  % put that into for loop || while loop  
n=ceil((log(top-bot)- log(err))/log(2))-1; %n为二分法运算总的次数；ceil是上取整 ,相对应的floor是下取整  
for i=0:1:n    % 步长默认为1，中间的1加不加都可以  
    x=(bot+top)/2;  
    fx=f(x);  
    if fx==0  
        bot=x; top=x;  
    elseif f(bot)*fx<0  
        top=x;  
    else  
        bot=x;  
    end  
%     if abs(top-bot)<=err  % 结合题意，加不加1/2都行，只是有着多二分运算一次的影响  
%         %此外，实际这里的if判断err的语句part，针对for循环是无必要的  
%         break  % 这里用break，不用return  
%     end  
end  
fprintf('\nThe result:\n二分法运算次数i=%d;方程的根x=%.4f;f(x)=%.4f\n',i,x,fx);%保留4位小数 

所以在想，可不可以把二分法扩展成多分法，反正也是等距的网格细化。

然后就试了一下，调了若干次吧，最后好在是成功了，不过这种方法有缺陷

1.只能搜索既定目标

2.要给出既定目标的大致范围

我没有试可不可以搜索两个以上的（相距不要太近）的目标，我认为两个目标如果角度差比较大的话还是可以实现的，这里就是学会这种细化网格的方法。

贴代码：

%2017.6.4
%author:Lola
%功能：实现单快拍稀疏矩阵DOA估计，采用网格细化减少运算量，单信号
%reference:A Sparse Signal Reconstruction Perspective for Source
%Localization With Sensor Arrays [Malioutov]
clc
clear all
close all

M = 8;                                                   %阵元数
K = 1;                                                   %信源数
L= 1;                                                    %快拍数
d_lamda =0.5;                                            %阵元间距半波长
w = [pi/4 ]';                                            %信号频率
theta1 = [7.25];                                         %信号来向
snr=20;                                                  %信噪比
grid_number=4;                                           %网格划分：每次在给定区间内画2等分
N=grid_number+1;
top=16;
boot=0;
err=0.01;

for k=1:K
    s=sqrt(10.^(snr/10))*[1+1j*w];                       %信号（信源数*快拍数）
    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;
X=awgn(X,snr);                                          %加入高白噪声

n=ceil((log(top-boot)- log(err))/log(grid_number))-1;   %n为网格细化总次数
step1=[];                                               %每次细化的步长
top1=[];                                                %网格最大点
boot1=[];                                               %网格最小点
for i=1:n                                               % 细化过程
    if theta1>=top || theta1<=boot;                     %判断给出区间是否包含要搜索的点
    disp('\n细化范围有误，无法继续运算，请重新调整区间端点.\n');
    return
end
    if i==1&&i~=n;
        top1(i)=top;
        top=top1(:,i);
        boot1(i)=boot;
        boot=boot1(:,i);
        step1(i)=(top1-boot1)/grid_number;                %第1次细化的步长
        step=step1(:,i);
    else if i~=n;
           step1(i)=step1(:,i-1)/grid_number;
           step=step1(:,i);                               %第i次细化的步长
           boot1(i)=Y(:,i-1);                             
           boot=boot1(:,i)-0.5*step1(:,i-1);            
           top1(i)=Y(:,i-1);
           top=top1(:,i)+0.5*step1(:,i-1);
        else
            YY=Y(:,i-1);
        end
    end   
    theta=boot:step:top;
    AA=[];                                                                %构造过完备基
    for kkk= 1:length(theta)
        g=exp(-1j*2*pi*[0:M-1]'*d_lamda*sin(theta(kkk)/180*pi));
        AA=[AA,g];
    end
    cvx_begin
    variable x(N) complex;
    minimize(square_pos(norm(X-AA*x,2))+2*norm(x,1));
    cvx_end
    Y=[];
    m=1;
    x1=conj(x);
    power=x1.*x;                                          %求信号能量
    xx=[(20*log10(power))].';                             %转化成功率
    [v,I]=findpeaks(xx);                                  %峰值搜索
    [vv,II]=sort(v,'descend');
    m=II;
    Y(i)=(I(m)-1)*step+boot;                              %第g次搜索的结果起点
end
fprintf('The estimation of DOA is YY=%.2f',YY);

结果：