方法:
这里采用的是等距划分。
想了很久也试了很多,最后是二分法给我的启发
二分法的方法参照:
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);
结果: