function [edoa, noPeaksFound] = doapeak(doaSpec, method, in1, in2)%DOAPEAK Finds the peaks of a DOA spectrum and the number of peaks.%%--------%Synopsis:%  [edoa, noPeaksFound] = doapeak(doaSpec, method, ...)%%Description:%  Estimates doas and number of doas from any DOA spectrum. All methods%  searches for peaks in the DOA spectrum.%%Output:%  edoa  (Vector of DoaT): Estimated directions [radians]. The number of%    DOA:s in this vector is equal to the value of the input parameter%    noSrc unless no spectrum peaks were found, in which case the length%    of this vector is 1.%  noPeaksFound (IntScalarT): Number of peaks found. This number is always%    less than or equal to the value of the input parameter "noSrc" for%    the methods 'amir' and 'amirbin' (see below).%%Input:%  doaSpec (DoaSpecT): Input DOA-spectrum. It is a quadratic measure of the%    presence of sources in different directions. It can be either a power%    spectrum (e.g. cbf) or a pseudo spectrum (e.g. music).%  method (StringT): See below.%    = 'amir': Amir Heydarkhan's method.%    = 'amirbin': Modified Amir Heydarkhan's method (Binary Search).%    = 'MA1': Moving average filtered search.%    = 'fam': Fredrik Athely's method.%%--------%Modified Amir Heydarkhan's method (Binary Search) ('amirbin')%Synopsis%  [edoa, noPeaksFound] = doapeak(doaSpec, 'amirbin', noSrc, startTreshold)%%Description:%  Same as the method 'amir' but with improved search scheme (binary).%  This may leed to different result then with the method 'amir'. This is due%  to that the methods have different convergence speed but with the same%  number of iterations. Compared to 'amir' there are for 'amirbin' some%  extra input parameters.%%Input:%  noSrc [D](IntScalarT): Number of sources or peaks we want to find.%  startTreshold [D](RealScalarT): Normalized start threshold. A number%    between 0 and 1.0 specifying the start threshold. This is then made%    higher and lower by the method. If the threshold falls below %    "startTreshold", the method is terminated.%%Algoritm:%  The algoritm only looks at those parts of the input DOA-spectrum%  that have values larger than a threshold.  Peaks in those parts%  are detected as having a derivative (from the Matlab function "diff")%  with a sign change + to -. If too few peaks are found, the threshold%  is lowered. If too many peaks are found, the threshold is increased.%  The threshold is changed a maximum of 100 times.%%--------%Amir Heydarkhan's method  ('amir')%Synopsis%  [edoa, noPeaksFound] = doapeak(doaSpec, 'amir', noSrc)%%Description:%  Finds the peaks of a DOA spectrum and the number of peaks (if less or%  equal to the value of the input parameter noSrc).%%  Iterates 100 times with a search algoritm that uses steps of%  size 1.2 when decreasing the search value and 1.1 when increasing the value%  i.e. value/1.2 respectively value*1.1.%%  A tip is to take the logaritm of the DOA spectrum when using Capon's %  beamformer before calling "doapeak".%%  It is suggested to use the method 'amirbin' instead of 'amir', since%  the former executes faster and safer and any improvement is done on this %  (former) method.%%Input:%  noSrc (IntScalarT): Number of sources or peaks we want to find.%%Algoritm:%  The algoritm only looks at those parts of the input DOA-spectrum%  that have values larger than a threshold.  Peaks in those parts%  are detected as having a derivative (from the Matlab function "diff")%  with a sign change + to -. If too few peaks are found, the threshold%  is lowered to 83 percent of the previous value of the threshold and%  a new search for peaks is performed. If too many peaks are found, the%  threshold is increased to 110 percent of the previous value  and a%  new search is performed. The threshold is changed a maximum of 100 times.%  If there are still too many peaks, a warning message is issued. If%  there are still too few peaks, the DOA of the last found peak is copied%  to the places of the not found peaks in the output parameter "edoa".%  In this way the size of "edoa" is the expected (given by "noSrc").%  In this case of to few found DOA:s, the number of found DOA: is given%  by the output parameter "noPeaksFound". The start value of the threshold%  is 1 percent of the maximum value of the DOA-spectrum.%%Known Bugs:%  The following execution error occured once in the methods 'amir':%    ???  Index exceeds matrix dimensions.%    Error in ==> /home/magic/svabj/lob/dbt2/dbt/doapeak.m%    On line 246  ==>           while (sdQ(k)==0)%    ...%    Error in ==> /home/magic/svabj/lob/dgaantal/run46.m%    On line 23  ==> run8cpi('run46',0,2,noTrial,snr,41)%%--------%Moving average filtered search ('MA1')%Synopsis%  [edoa, noPeaksFound] = doapeak(doaSpec, 'MA1', sampleT)%%Description:%  Searches for peaks in the DOA spectrum. (filtered)%%Input:%  sampleT (IntScalarT): Sample period in radians. i.e. angle between samples.%%Algoritm:%  The spectrum is first filtered with the requested filter then a search%  is performed to find the peaks of the curve by loooking at sign change%  of the derivative of the DOA spectrum. The filtering is used to reduce%  peaks that due to noise and likewise.%%--------%Fredrik Athley's Method ('fam').%Synopsis:%  [edoa, noPeaksFound] = doapeak(doaSpec, 'fam')%  [edoa, noPeaksFound] = doapeak(doaSpec, 'fam', thFix, thAdapt)%%Description:%  Estimates directions from a directional spectrum.%%Input:%  thFix [D](RealScalarT): The fixed threshold [dB].%  thAdapt  [D](RealScalarT): The relative adaptive threshold [dB].%%Algoritm:%  Only the values of the logaritmic DOA spectra above a theshold are%  considered. The theshold is determined as the largest of a fixed%  level and an adaptive level. The latter is set to be "thAdapt" (default=20)%  dB below the highest peak of the spectrum. The default value for the fixed %  level is 2 dB but can be changed by the input parameter "thFix".%  The estimated directions are given by the local maxima above the theshold. %  The code in this method is copied from the function "bearingextractor"%  with minor modifications.%%--------%Notations:%  Data type names are shown in parentheses and they start with a capital%  letter and end with a capital T. Data type definitions can be found in [1]%  or by "help dbtdata".%  [D] = This parameter can be omitted and then a default value is used.%  When the [D]-input parameter is not the last used in the call, it must be%  given the value [], i.e. an empty matrix.%  ... = There can be more parameters. They are explained under respective%  metod or choice.%%Software Quality:%  (About what is done to ascertain software quality. What tests are done.)%  This function was used in the work done by A. Heydarkhan [4].%%Known Bugs:%%References:%  [1]: Bj鰎klund S.: "DBT, A MATLAB Toolbox for Radar Signal Processing.%    Reference Guide", FOA-D--9x-00xxx-408--SE, To be published.%  [4]: Heydarkhan A.: "Model Based Direction of Arrival Estimation Methods%    Applied to Experimental Antenna Data", Master's Thesis, Chalmers%    University of Technology 1997.%  [5]: Athley F.:"Integration av modellbaserade metoder i radar, ett%    till鋗pningsexempel", Rapport FY/P-98:015, Ericsson Microwave Systems 1998.%%See Also:%  spanosrc, sdoaspc, doaspc1%   *  DBT, A Matlab Toolbox for Radar Signal Processing  *% (c) FOA 1994-2000. See the file dbtright.m for copyright notice.%%  Start        : 970613 Amir Heydarkhan (amihey).%  Latest change: $Date: 2000/10/16 15:20:42 $ $Author: svabj $.%  $Revision: 1.27 $% *****************************************************************************% ----------------------------------------------------------------------- %% Handle input parameters% ----------------------------------------------------------------------- %arginNo=2;if (nargin < arginNo)  dbterror('To few input parameters.')endarginNo = arginNo +1;% ****************** Add missing input parameters ******************% ****************** Default values ******************% ****************** Error check input parameters ******************chkdtype(doaSpec, 'DoaSpecT')chkdtype(method, 'StringT')% ----------------------------------------------------------------------- %% Method: amir% ----------------------------------------------------------------------- %if strcmp(method,'amir')  noSrc = in1;  % ****************** Pick out fields from input parameters. ******************  Q = doaSpec.specSmpl;  if noSrc == 1,    [Y, edoaIx] = max(Q);  else    edoaIx = 0;    l=0;    jumplimit = 100;    jump = 0;    startTreshold = .5;    treshold = startTreshold;    Qmax = max(Q);    %interval = 2    %Q(1:interval) = diag(zeros(interval));     %smooth the ends    %Q(size(Q,1)-interval+1:size(Q,1)) %= diag(zeros(interval));    while ((size(edoaIx,1) ~= noSrc) & (jump < jumplimit))      edoaIx = 0;      l=0;      for k = 1 : size(Q,1),        %Smoothing        if Q(k) < Qmax*treshold,          Q2(k)=0;        else          Q2(k)=Q(k);        end      end      dQ = diff(Q2);      dQ(size(dQ,1)+1) = 0;      sdQ = sign(dQ);                  %Searching      for k = 2:length(sdQ)-1,        if sdQ(k-1) == 1,          while (sdQ(k)==0)            k=k+1;          end % while.          if sdQ(k) == -1,            l=l+1;            edoaIx(l) = k;          end        end      end      %fprintf('%d st har hittat\n',size(edoaIx,1));      if length(edoaIx) > noSrc,        %disp('謐ar treshold')        treshold = treshold*1.1;        if treshold > 1,          jump = jumplimit;        end      elseif length(edoaIx) < noSrc,        %disp('Minskar treshold')          if treshold < startTreshold,            jump = jumplimit;          else,            treshold = treshold/1.2;          end      end      jump = jump + 1;    end%while  end%ifif (length(edoaIx)>noSrc)  %error('DBT-Error: More peaks found than desired.'); dbtwarning('More peaks found than desired.');end%if% ----------------------------------------------------------------------- %% Method: amirbin. Same as amir but with binary search scheme.% ----------------------------------------------------------------------- %elseif strcmp(method,'amirbin') % ****************** Add missing input parameters ******************  if (nargin < arginNo)    in1 = [];  end  arginNo = arginNo +1;  if (nargin < arginNo)    in2 = [];  end  arginNo = arginNo +1;   % ****************** Default values ******************  noSrc = in1;  startTreshold = in2;  if isempty(noSrc)    noSrc = 4;  end%if  if isempty(startTreshold)    startTreshold = 0.5;  end%if  % ****************** Error check input parameters ******************  chkdtype(noSrc, 'IntScalarT')    % ****************** Do the work ******************  Q = doaSpec.specSmpl;  if noSrc == 1,    [Y, edoaIx] = max(Q);  else    edoaIx = 0;    l=0;    jumplimit = 100;    jump = 0;    %startTreshold = .5;    treshold = startTreshold;    Qmax = max(Q);    while ((length(edoaIx) ~= noSrc) & (jump < jumplimit))      edoaIx = 0;      l=0;      for k = 1 : size(Q,1),        %Smoothing        % the parts of the spectrum that are below the current theshold        % are set to zero.        if Q(k) < Qmax*treshold,          Q2(k)=0;        else          Q2(k)=Q(k);        end      end      % Derivative.      dQ = diff(Q2);      dQ(size(dQ,1)+1) = 0;      sdQ = sign(dQ);      %Searching      for k = 2:length(sdQ)-1,        if sdQ(k-1) == 1,          while ((sdQ(k)==0) & (k < length(sdQ)))            k=k+1;          end % while.          if sdQ(k) == -1,            l=l+1;            edoaIx(l) = k;          end        end      end      if length(edoaIx) > noSrc,        %disp('Increase treshold')        treshold = treshold/2 + 0.5;        if treshold > 1,          jump = jumplimit;        end%if      elseif length(edoaIx) < noSrc,        %disp('Decrease treshold')          if treshold < startTreshold,            jump = jumplimit;          else,            treshold = 1.5*treshold-0.5;          end%if      end%if      jump = jump + 1;    end %while  end %if% ----------------------------------------------------------------------- %% Method: Moving Average Filter 1.% ----------------------------------------------------------------------- %elseif strcmp(method,'MA1')  Q = doaSpec.specSmpl;    % The sampled DOA Spectrum.  sampleT = in1;  % Sample period time in degrees.  MALen = 0.5;    % Moving Average filter length in degrees.  MACoeff = ones(1,MALen/sampleT)*(sampleT/MALen);  Qmax=max(Q);  for k = 1 : length(Q),    if Q(k) < Qmax*0.01,      Q(k)=0;    end  end  Qf = filter(MACoeff,1,Q);  dQf = diff(Qf);  dQf(length(dQf)+1) = 0;  sdQf = sign(dQf);  l=0;  for k = 2:length(sdQf)-1,    if sdQf(k-1) == 1,      while (sdQf(k)==0)        k=k+1;      end % while.      if sdQf(k) == -1,        l=l+1;        edoaIx(l) = k;      end    end  end% ----------------------------------------------------------------------- %% Method: fam, Fredrik Athley's Method%% Start: 9xxxxx Fredrik Athley (freath).% ----------------------------------------------------------------------- %elseif strcmp(method,'fam')  %function [azest,notgt] = bearingextractor(P,theta,thFix)  % ****************** Add missing input parameters ******************  if (nargin < arginNo)    in1 = [];  end  arginNo = arginNo +1;  if (nargin < arginNo)    in2 = [];  end  arginNo = arginNo +1;  % ****************** Default values ******************  thFix = in1;  thAdapt = in2;  if isempty(thFix)   thFix = 2;		% Fixed threshold [dB].  end%if  if isempty(thAdapt)   thAdapt = 20;	% The adaptive threshold is "thAdapt" dB below the                        % highest peak in the spectrum.  end%if  % ****************** Error check input parameters ******************  chkdtype(thFix, 'RealScalarT')  chkdtype(thAdapt, 'RealScalarT')  % ****************** Do the work ******************  edoaIx = 0;  P = 10*log10(doaSpec.specSmpl);    % The sampled DOA spectrum in dB.  theta = doaSpec.doaPos;  specLen = length(P);		% Length of DOA spectrum.	  notgt = 0;  pmax = max(P);  th = max(thFix,pmax-thAdapt);	% Applied threshold in dB.  for n = 2:specLen-1    exwin =[P(n-1);P(n+1)];    if (P(n)>max(exwin)) & (P(n)>th)      notgt = notgt + 1;      edoaIx(notgt) = n;    end%if  end%for n% ----------------------------------------------------------------------- %% Method:% ----------------------------------------------------------------------- %end % if method.% ----------------------------------------------------------------------- %% Deliverance of calculated values.% ----------------------------------------------------------------------- %% edoaIx (IndexT): Index for the found peaks in the spectrum.% If no peaks are found edoaIx == 0 here.if (edoaIx ~= 0)  ang = doaSpec.doaPos;  etheta = ang(edoaIx);  edoa=[etheta; zeros(1,length(etheta))];  noPeaksFound = length(edoaIx);else  edoa=0;  noPeaksFound = 0;%  plot(Q);end%if%End Of File