www.pudn.com > array.zip > circ_array.m, change:2003-07-07,size:5959b


function [pattern,amn] = circ_array(N,dolxr,dolyr,theta0,phi0,winid,win,nbits); 
%%%%%%%%%%%%%%%%%%%% ************************ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This function computes the 3-D directive gain patterns for a circular plannar array 
% This function uses the fft2 to compute its output. It assumes that there are the same number  
% of elements along the major x- and y-axes 
%%%%%%%%%%%%%%%%%%%% ************  INPUTS ************ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% N ==> number of elements along x-aixs or y-axis 
% dolxr ==> element spacing in x-direction; dolyr ==> element spacing in y-direction. Both are in lambda units 
% theta0 ==> elevation steering angle in degrees, phi0 ==> azimuth steering angle in degrees 
% This function uses the function (rec_to_circ) which computes the circular array from a square  
% array (of size NXN) using the notation developed by ALLEN,J.L.,"The Theory of Array Antennas  
% (with Emphasis on Radar Application)" MIT-LL Technical Report No. 323,July, 25 1965.  
% winid ==> window identifier; winid negative ==> no window ; winid positive ==> use window given by win 
% win ==> input window function (2-D window) MUST be of size (Nxr X Nyr) 
% nbits is the number of nbits used in phase quantization; nbits negative ==> NO quantization 
%%%%%%%%%%%%%%%%%%%% *********** OUTPUTS ************* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% amn ==> array of ones and zeros; ones indicates true element location on the grid 
% zeros mean no elements at that location; pattern ==> directive gain pattern 
%%%%%%%%%%%%%%%%%%%% ************************ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
eps = 0.0001; 
nx = 0:N-1; 
ny = 0:N-1; 
i = sqrt(-1); 
 
% check that window size is the same as the array size 
[nw,mw] = size(win); 
if winid >0 
   if mw ~= N 
      fprintf('STOP == Window size must be the same as the array') 
       return 
   end 
   if nw ~= N 
      fprintf('STOP == Window size must be the same as the array') 
       return 
   end 
end 
 
%if dol is > 0.5 then; choose dol = 0.5 and compute new N 
if(dolxr <=0.5) 
   ratiox = 1 ; 
   dolx = dolxr ; 
   Nx = N ; 
else 
   ratiox = ceil(dolxr/.5) ; 
   Nx = (N-1) * ratiox + 1 ; 
   dolx = 0.5 ; 
end 
if(dolyr <=0.5) 
   ratioy = 1 ; 
   doly = dolyr ; 
   Ny = N ; 
else 
   ratioy = ceil(dolyr/.5);  
   Ny = (N-1)*ratioy + 1 ; 
   doly = 0.5 ; 
end 
 
% choose proper size fft, for minimum value choose 256X256 
Nrx = 10 * Nx;  
Nry = 10 * Ny; 
nfftx = 2^(ceil(log(Nrx)/log(2))); 
nffty = 2^(ceil(log(Nry)/log(2))); 
if nfftx < 256 
    nfftx = 256; 
end 
if nffty < 256 
    nffty = 256; 
end 
 
% generate array of elements with or without window 
if winid < 0  
   array = ones(N,N); 
else 
    array = win; 
end 
 
% convert steering angles (theta0, phi0) to radians 
theta0 = theta0 * pi / 180; 
phi0 = phi0 * pi / 180; 
% convert steering angles (theta0, phi0) to U-V sine-space 
u0 = sin(theta0) * cos(phi0); 
v0 = sin(theta0) * sin(phi0); 
 
% Use formula thetal = (2*pi*n*dol) * sin(theta0) divided into 2^m levels 
% and rounded to the nearest qunatization level 
if nbits < 0 
    phasem = exp(i*2*pi*dolx*u0 .* nx * ratiox); 
    phasen = exp(i*2*pi*doly*v0 .* ny * ratioy); 
else 
    levels = 2^nbits; 
    qlevels = 2.0*pi / levels; % compute quantization levels 
    sinthetaq = round(dolx .* nx * u0 * levels * ratiox) .* qlevels; % vector of possible angles 
    sinphiq = round(doly .* ny * v0 * levels *ratioy) .* qlevels; % vector of possible angles 
    phasem = exp(i*sinthetaq); 
    phasen = exp(i*sinphiq);      
end 
 
% add the phase shift terms 
array = array .* (transpose(phasem) * phasen) ; 
 
% determine if interpolation is needed (i.e N > Nr) 
if (Nx > N )| (Ny > N) 
   for xloop = 1 : N 
      temprow = array(xloop, :) ; 
      w( (xloop-1)*ratiox+1, 1:ratioy:Ny) =  temprow ; 
   end 
   array = w; 
else 
    w(1:Nx, :) = array(1:N,:); 
end 
 
% Convert rectangular array into circular using function rec_to_circ 
[m,n] = size(w) ; 
NC = max(m,n);  % Use Allens algorithm 
if Nx == Ny 
    temp_array = w; 
else 
    midpoint = (NC-1)/2 +1 ; 
    midwm = (m-1)/2 ; 
    midwn = (n-1)/2 ; 
    temp_array = zeros(NC,NC); 
    temp_array(midpoint-midwm:midpoint+midwm, midpoint-midwn:midpoint+midwn) = w ; 
end 
amn = rec_to_circ(NC);  % must be rectangular array (Nx=Ny) 
amn = temp_array .* amn ; 
 
% Compute array pattern 
arrayfft = abs(fftshift(fft2(amn,nfftx,nffty))).^2 ; 
 
%compute [su,sv] matrix 
U = [-nfftx/2:(nfftx/2)-1] ./(dolx*nfftx); 
indexx = find(abs(U) <= 1); 
U = U(indexx); 
V = [-nffty/2:(nffty/2)-1] ./(doly*nffty); 
indexy = find(abs(V) <= 1); 
V = V(indexy); 
[SU,SV] = meshgrid(V,U); 
indx = find((SU.^2 + SV.^2) >1); 
arrayfft(indx) = eps/10; 
 
 
%Normalize to generate gain patern 
rbar=sum(sum(arrayfft(indexx,indexy))) / dolx/doly/4./nfftx/nffty; 
arrayfft = arrayfft(indexx,indexy) ./rbar; 
[SU,SV] = meshgrid(V,U); 
indx = find((SU.^2 + SV.^2) >1); 
arrayfft(indx) = eps/10; 
pattern = 10*log10(arrayfft +eps); 
 
figure(1) 
mesh(V,U,pattern); 
xlabel('V') 
ylabel('U'); 
zlabel('Gain pattern - dB') 
 
figure(2) 
contour(V,U,pattern) 
axis image 
grid 
xlabel('V') 
ylabel('U'); 
axis([-1 1 -1 1]) 
 
figure(3) 
 
x0 = (NC+1)/2 ; 
y0 = (NC+1)/2 ; 
radiusx = dolx*((NC-1)/2 + 0.05/dolx) ; 
radiusy = doly*((NC-1)/2 + 0.05/dolx) ; 
theta = 5  ; 
[xxx, yyy]=find(abs(amn)>0); 
xxx = xxx-x0 ; 
yyy = yyy-y0 ; 
plot(yyy*doly, xxx*dolx,'rx') 
axis equal 
hold on 
axis([-radiusy-0.5 radiusy+0.5 -radiusx-0.5  radiusx+0.5]); 
grid 
title('antenna spacing pattern'); 
xlabel('y - \lambda units') 
ylabel('x - \lambda units') 
 
 
[x, y]= makeellip( 0, 0, radiusx, radiusy, theta) ; 
 
plot(y, x) ; 
axis([-radiusy-0.5 radiusy+0.5 -radiusx-0.5  radiusx+0.5]); 
 
 
[xxx0, yyy0]=find(abs(amn)<=0); 
xxx0 = xxx0-x0 ; 
yyy0 = yyy0-y0 ; 
plot(yyy0*doly, xxx0*dolx,'co') 
axis([-radiusy-0.5 radiusy+0.5 -radiusx-0.5  radiusx+0.5]); 
axis equal 
 
hold off ; 
 
return