beampatt.m

书籍Detection, Estimation, and Modulation Theory的Matlab源代码；

    s=get(e_w,'string');
    q=eval(s);
    set(e_w,'UserData',q,'String',num2str(q,'%2.2f'));
    BP_DAT(6)=q;

elseif strcmp(action,'redraw'),
    N=BP_DAT(1);
    d=BP_DAT(2);
    f=BP_DAT(3);
    us=BP_DAT(4);
    wtype=BP_DAT(5);
    q=BP_DAT(6);
    h1=BP_DAT(7);
    h2=BP_DAT(8);
    h3=BP_DAT(9);
    h4=BP_DAT(10);
    h5=BP_DAT(11);


    D  = [-(N-1)/2:1:(N-1)/2].';  % sensor locations

    vs = exp(j*2*pi*d*f*D*us);
    u = [-3:0.01:3];
    v = exp(j*2*pi*d*f*D*u);

    th = [-1:0.01:1]*pi;
    vth = exp(j*2*pi*d*f*D*cos(th));
    if wtype == 1               % uniform
      Wq = ones(N,1)/N;
    elseif wtype ==2            % uniform difference  
      if rem(N,2)==0                      % even
        Wq = [-ones(N/2,1);ones(N/2,1)]/N;
      else
        Wq = [-ones((N-1)/2,1);0;ones((N-1)/2,1)]/N;
      end
    elseif wtype ==3            % triangular (sampling)
        Wq = ones(N,1)-abs(D)*2/N;
        Wq = Wq/sum(Wq);
     elseif wtype ==4            % triangular (root matching)
        Wq = ones(N,1)-abs(D)*2/(N+1);
        Wq = Wq/sum(Wq);
     elseif wtype ==5            % raised cosine 
        p=q;
        if p<0
           p=0;
        elseif p>1
        		p=1;   
        end           
        Wq = (1-p)*(cos(pi*D/N))+p;
        Wq = Wq/sum(Wq);
     elseif wtype ==6            % cosine ^m 
        m=q;
        if m<0
           m=0;
		  end           
        Wq = (cos(pi*D/N)).^m;
        Wq = Wq/sum(Wq);
     elseif wtype ==7            % raised cosine^2 
        p=q;
        if p<0
           p=0;
        elseif p>1
        		p=1;   
        end           
        Wq = (1-p)*(cos(pi*D/N)).^2+p;
        Wq = Wq/sum(Wq);
     elseif wtype ==8            % Hamming 
        Wq = 0.46*(cos(2*pi*D/N))+0.54;
        Wq = Wq/sum(Wq);
     elseif wtype ==9            % Blackman-Harris 
        Wq = 0.42+0.5*(cos(2*pi*D/N))+0.08*(cos(4*pi*D/N));
        Wq = Wq/sum(Wq);
     elseif wtype ==10            % DPSS 
        u0 = q;
        if u0 <=0
           q = 0;
           Wq = ones(N,1)/N;
        else
           r = 2*pi*u0*sinc([0:1:N-1]*u0);
           A = toeplitz(r,r);
           [vv,l]=eig(A);
           [l,ind]=sort(diag(l));
           vv=vv(:,ind);
           Wq=vv(:,N);
        end
        Wq = Wq/sum(Wq);
     elseif wtype ==11            % Kaiser 
         beta=q;
        if beta<0
           beta=0;
        end           
        Wq = besseli(0,beta*sqrt(1-(2*D/(N-1)).^2));
        Wq = Wq/sum(Wq);
    end 
    BP1 = real((vs.*Wq)'*v);
    BP2 = imag((vs.*Wq)'*v);
    BP3 = db(abs((vs.*Wq)'*v).^2);
    BP4 = angle((vs.*Wq)'*v)/pi;
    BP5 = db(abs((vs.*Wq)'*vth).^2);

    % set -pi to pi
    eps = 1e-5;
    [y,ind]=find(BP4<= -1+eps);
    n1=length(ind);
    BP4(ind) = ones(1,n1);

    % set angle at BP=0 to avg of other terms
    [y,ind]=find(BP3<-100);
    n0 = length(ind);
    for n=1:n0
     if ind(n) == 1
        BP4(1)=BP4(2);
     elseif ind(n)==length(u);
        BP4(ind(n)) = BP4(ind(n)-1);
     else
        BP4(ind(n)) = 0.5*(BP4(ind(n)-1)+BP4(ind(n)+1));
     end
    end 
   
    %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  Initialize Simulation and plots
    %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    FS=12;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Print Values
    %FS = 6;
    %TS = 6;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fillc = 0.8*[1 1 1];
       axes(h3)
       fill([1 3 3 1],[-1 -1 1 1],fillc)
       hold on 
       fill(-[1 3 3 1],[-1 -1 1 1],fillc)
     %  plot([-1:0.25:1].'*[1 1],[-1 1],':k')
       plot([-3 3],[0 0],'--k')
       p1=plot(u,BP1);
       axis([-3 3 -1 1])
       plot([us us],[-1 1],'g');
       set(p1,'LineWidth',1.0)
       plot([1 1],[-1 1],'r')
       plot([-1 -1],[-1 1],'r')
       hold off
       grid on
   %    title('Real','FontSize',TS)
       ylabel('Real')
      set(gca,'FontSize',FS)
       set(gca,'Layer','top')
       set(gca,'Xtick',[-3:0.5:3] )
       get(gca,'FontSize');
       axes(h4)
       fill([1 3 3 1],[-1 -1 1 1],fillc)
       hold on 
       fill(-[1 3 3 1],[-1 -1 1 1],fillc)
   %    plot([-1:0.25:1].'*[1 1],[-1 1],':k')
       plot([-3 3],[0 0],'--k')
       p2=plot(u,BP2);
       get(p2,'LineWidth');
       set(p2,'LineWidth',1.0);
       axis([-3 3  -1 1])
       plot([us us],[-1 1],'g')
       plot([1 1],[-1 1],'r')
       plot([-1 -1],[-1 1],'r')
       hold off
       grid on
   %    title('Imaginary','FontSize',TS)
       ylabel('Imaginary')
    set(gca,'FontSize',FS)
       set(gca,'Layer','top')
       set(gca,'Xtick',[-3:0.5:3] )
       xlabel('u=cos(theta)')

       axes(h1)
       fill([1 3 3 1],[-60 -60 0 0],fillc)
       axis([-3 3 -60 0])
       hold on 
       fill(-[1 3 3 1],[-60 -60 0 0],fillc)
       p3=plot(u,BP3);
       hold on
       plot([us us],[-60 0],'g')
       plot([1 1],[-60 0],'r')
       plot([-1 -1],[-60 0],'r')
   %    plot([-1:0.25:1].'*[1 1],[-60 0],':k')
       hold off
       grid on
       ylabel('dB')
       %title('dB','FontSize',TS)
       set(gca,'FontSize',FS)
       set(gca,'Layer','top')
       set(gca,'Xtick',[-3:0.5:3] )

       axes(h2)
       fill([1 3 3 1],[-1.1 -1.1 1.1 1.1],fillc)
       hold on 
       fill(-[1 3 3 1],[-1.1 -1.1 1.1 1.1],fillc)
       p4=plot(u,BP4);
%       p4=plot(u,abs(BP));
       axis([-3 3 -1.1 1.1])
       plot([us us],[-1.1 1.1],'g')
       plot([1 1],[-1.1 1.1],'r')
       plot([-1 -1],[-1.1 1.1],'r')
    %   plot([-1:0.25:1].'*[1 1],[-1.1 1],':k')
       hold off
       grid on
       ylabel('Phase')
      % title('Phase','FontSize',TS)
       set(gca,'FontSize',FS)
       set(gca,'Layer','top')
       set(gca,'Ytick',[-1 -0.5 0 0.5 1])
       set(gca,'YTickLabel',['-p  '; '-p/2'; ' 0  '; ' p/2'; ' p  '],...
	'FontName','Symbol')
       set(gca,'Xtick',[-3:0.5:3] )

       axes(h5)
       p5=polardb(th,BP5,-40);
       hold on
       polardb(acos([us us]),[-40 0],-40,'g')
       polardb(-acos([us us]),[-40 0],-40,'g')
       hold off
       xlabel('theta')

    drawnow;


elseif strcmp(action,'done'),
    close(gcf);
    clear global BP_DAT

end

function y = db(x)
y = 10*log10(x);


function y=sinc(x)
y=ones(size(x));
i=find(x);
y(i)=sin(pi*x(i))./(pi*x(i));


function hpol=polardb(theta,rho,lim,line_style)
% polardb(theta,rho,lim,linestyle)
% Polardb is a modified version of the regular 'polar' function.
% Plotting range is -180 to 180 deg with zero at top.
% Theta increases in clockwise manner.
% Inputs:
%   theta - angles in RADIANS (although axes labeled in degrees)
%   rho   - plot value in dB
%   lim   - lower limit for plot in dB (e.g. -40)
%   line_style - string indicating line style, e.g. '-g' (optional)
% Example polardb(theta,beampatt,-40,'-r')
% last modified 9/23/98 K. Bell
%
%POLAR	Polar coordinate plot.
%	POLAR(THETA, RHO) makes a plot using polar coordinates of
%	the angle THETA, in radians, versus the radius RHO.
%	POLAR(THETA,RHO,S) uses the linestyle specified in string S.
%	See PLOT for a description of legal linestyles.
%
%	See also PLOT, LOGLOG, SEMILOGX, SEMILOGY.

%	Copyright (c) 1984-94 by The MathWorks, Inc.

const=0.7;

if nargin < 3
	error('Requires 3 or 4 input arguments.')
elseif nargin == 3 
        [n,n1]=size(theta);
	if isstr(rho)
		line_style = rho;
		rho = theta;
		[mr,nr] = size(rho);
		if mr == 1
			theta = 1:nr;
		else
			th = (1:mr)';
			theta = th(:,ones(1,nr));
		end
	else
		line_style = 'auto';
	end
elseif nargin == 2
	line_style = 'auto';
	rho = theta;
	[mr,nr] = size(rho);
	if mr == 1
		theta = 1:nr;
	else
		th = (1:mr)';
		theta = th(:,ones(1,nr));
	end
end
if isstr(theta) | isstr(rho)
	error('Input arguments must be numeric.');
end
if any(size(theta) ~= size(rho))
	error('THETA and RHO must be the same size.');
end

nr = size(rho,2);
tck = -floor(lim/10);
lim = -tck*10;
I=find(rho<lim);
ni=size(I,2);
rho(I)=lim*ones(1,ni);
rho = rho/10+tck*ones(1,nr);

% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;

% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');

% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle  = get(cax, 'DefaultTextFontAngle');
fName   = get(cax, 'DefaultTextFontName');
fSize   = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
	'DefaultTextFontName',   get(cax, 'FontName'), ...
	'DefaultTextFontSize',   get(cax, 'FontSize'), ...
	'DefaultTextFontWeight', get(cax, 'FontWeight') )

% only do grids if hold is off
if ~hold_state

% make a radial grid
	hold on;
	hhh=plot([0 max(theta(:))],[0 max(abs(rho(:)))]);
	v = [get(cax,'xlim') get(cax,'ylim')];
	ticks = length(get(cax,'ytick'));
	delete(hhh);
% check radial limits and ticks
	rmin = 0; rmax = v(4); rticks = ticks-1;
	if rticks > 5	% see if we can reduce the number
		if rem(rticks,2) == 0
			rticks = rticks/2;
		elseif rem(rticks,3) == 0
			rticks = rticks/3;
		end
	end

% define a circle
	th = 0:pi/50:2*pi;
	xunit = cos(th);
	yunit = sin(th);
% now really force points on x/y axes to lie on them exactly
    inds = [1:(length(th)-1)/4:length(th)];
    xunits(inds(2:2:4)) = zeros(2,1);
    yunits(inds(1:2:5)) = zeros(3,1);

%	rinc = (rmax-rmin)/rticks;
        rinc = const;
%	for i=(rmin+rinc):rinc:rmax
	for i=[1:1:tck]
%	for i=rmax:rmax
		plot(const*xunit*i,const*yunit*i,'--','color',tc,'linewidth',0.5);
%		text(0,i+rinc/20,['  ' num2str(10*(i-tck))],'verticalalignment','bottom' );
		text(const*(-i+rinc/100),0,['  ' num2str(10*(i-tck))],'verticalalignment','bottom' );
	end

% plot spokes
	th = (1:6)*2*pi/12;
%	th = (1:2)*2*pi/4;
	cst = cos(th); snt = sin(th);
	cs = [-cst; cst];
	sn = [-snt; snt];
	plot(const*rmax*cs,const*rmax*sn,'--','color',tc,'linewidth',0.5);
       
% annotate spokes in degrees
%	rt = 1.1*rmax;
	rt = 1.15*rmax;
	for i = 1:max(size(th))
	        text(const*rt*snt(i),const*rt*cst(i),int2str(i*30),'horizontalalignment','center' );
	        
  		loc = int2str(i*30-180);
%		if i == max(size(th))
%			loc = int2str(0);
% 		end
		text(-const*rt*snt(i),-const*rt*cst(i),loc,'horizontalalignment','center' );
	end

% set viewto 2-D
	view(0,90);
% set axis limits
	axis(rmax*[-1 1 -1.1 1.1]);
end

% Reset defaults.
set(cax, 'DefaultTextFontAngle', fAngle , ...
	'DefaultTextFontName',   fName , ...
	'DefaultTextFontSize',   fSize, ...
	'DefaultTextFontWeight', fWeight );

% transform data to Cartesian coordinates.
yy = const*rho.*cos(theta);
xx = const*rho.*sin(theta);

% plot data on top of grid
if strcmp(line_style,'auto')
	q = plot(xx,yy);
else
	q = plot(xx,yy,line_style);
end
set(q,'LineWidth',1.0);
if nargout > 0
	hpol = q;
end
if ~hold_state
	axis('equal');axis('off');
end

% reset hold state
if ~hold_state, set(cax,'NextPlot',next); end