验证12GHztestdxdyclear all;close all;clcN=60;x方向阵元个数M=60;y方向阵元个数f=12;信号频率 单位GHzc=29979210^2;光速 单位mslambda=cf;波长lambdag = lambda 088 波导中波长dx=lambdag5;x方向单元间距dy=lambdag5;y方向单元间距phi=linspace-9090181; 方位角范

以下是修改后的程序：

clear all;
close all;
clc

N = 60; % x方向阵元个数
M = 60; % y方向阵元个数
f = 12; % 信号频率 单位GHz
c = 2.99792 * 10^8; % 光速 单位m/s
lambda = c / (f * 10^9); % 波长，需要将GHz转换为Hz
lambdag = lambda * 0.88; % 波导中波长
dx = lambdag / 5; % x方向单元间距
dy = lambdag / 5; % y方向单元间距
phi = linspace(-90, 90, 181); % 方位角范围
theta = linspace(-90, 90, 181); % 俯仰角范围
theta0 = 1; % 目标俯仰角
phi0 = 180; % 预先设定的方向，目标方位角
Amax = 1; % 用于幅值调整
Amin = 0.5; % 用于幅值调整
k = 1; % 用于循环变量
p = zeros(3, M * N); % 生成3行，m*n列的矩阵
for iii = 1 : N
    for jjj = 1 : M
        p(:, k) = [-(N + 1) / 2 * dx + iii * dx, -(M + 1) / 2 * dy + jjj * dy, 0]'; % 单行时应把对应坐标调为0 即在dx/dy上乘一个0
        k = k + 1;
    end
end   % 阵元的位置信息

% 在三维图中绘制所有阵元
figure(1);
plot3(p(1,:), p(2,:), p(3,:), 'ko');
hold on;
xlabel('\it x');
ylabel('\it y');
zlabel('\it z');

% 二维全息幅值计算
A0 = p; % 获取所有辐射点位置坐标信息可以去掉
P0 = 2 * pi / lambda * [sin(theta0 * pi / 180) * cos(phi0 * pi / 180), sin(theta0 * pi / 180) * sin(phi0 * pi / 180), cos(theta0 * pi / 180)].'; % 目标波除位置信息外信息
Q0 = P0.' * p; % 目标波全部信息
q = 1; % 用于循环变量
for qq = 1 : (M * N)
    if p(2, q) > 0
        R0(:, q) = Q0(:, q) - 2 * pi / lambdag * sqrt(p(1, q)^2 + p(2, q)^2) - pi; % 将一半部分补180°相位防止凹陷
    else
        R0(:, q) = Q0(:, q) - 2 * pi / lambdag * sqrt(p(1, q)^2 + p(2, q)^2);
    end
    q = q + 1;
end
R = R0; % 检查是否补180°相位  代入PCAAD时无需考虑反向 因为相位设置是延轴向设置
M0 = Amax + Amin * cos(R0); % 幅值进行缩放
m = 1; % 用于循环变量
n = length(M0);
for mm = 1 : length(M0)
    if M0(:, m) > 1
        M0(:, m) = 1;
    else
        M0(:, m) = 0;
    end

    m = m + 1;
end
m0 = M0;
% 以上对幅值进行归一化
MM0 = M0.'; % 用于pcaad幅值文件
MM1 = reshape(MM0, N, M); % 用于绘制黑白图

figure(5); % 绘制黑白格
[a, b] = size(MM1);
MM11 = zeros(a + 1, b + 1);
MM11(2:a + 1, 1:b) = MM1;
pcolor(flipud(1 - MM11))
colormap(gray(2))
axis square

% 以下为绘制方向图函数程序
v = zeros(m, n);
for ii = 1 : length(theta)
    for jj = 1 : length(phi)
        k = 2 * pi / lambda * [sin(theta(ii) * pi / 180) * cos(phi(jj) * pi / 180) + lambda / lambdag, sin(theta(ii) * pi / 180) * sin(phi(jj) * pi / 180) + lambda / lambdag, 0].';
        v = exp(1i * k.' * p); % x方向阵因子
        b(ii, jj) = abs(sum(v .* m0)) / sqrt(M * N);
    end
end
B1 = max(abs(b(:))); % 获取在目标方向的最大值用于验证
B = b(ii, jj); % 以下二维方向图中所有方向信息对数化--dB形式

figure(2); % 绘制方位角方向
plot(phi, 20 * log10(abs(B(:, theta0))), 'k-')
grid on, hold on
xlabel('方位角(deg)'), ylabel('阵列方向图(dB)')
axis([-90 90 -50 10]); % 可为x轴和y轴设置一个极限范围

figure(3); % 绘制俯仰角方向图
plot(theta, 20 * log10(abs(B(phi0, :))))
grid on, hold on
xlabel('俯仰角(deg)'), ylabel('阵列方向图(dB)')
axis([-90 90 -50 10]); % 可为x轴和y轴设置一个极限范围

figure(4); % 绘制3D方向图
mesh(B)
title('圆阵动态30度指向方向图');
xlabel('俯仰角/（\circ）');
ylabel('方位角/（\circ）');
zlabel('阵列方向图(dB)');
axis([-90 90 -90 90 -50 10]); % 可为x轴和y轴设置一个极限范围
grid on

请注意，修改后的程序中，将频率f的单位由GHz转换为Hz，即乘以10^9。并修正了计算幅值归一化的代码