UPA beforming Matlab 仿真

未考虑时延

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
clear;
clc;
% clf;
% close all;

% coordinate theta: 0-180; phi = -180--180
dB = 30;
M = 100;
theta_range = linspace(-pi/2,pi/2,M); % elevation angle 
phi_range = linspace(-pi/2,pi/2,M);

DUT_array = 'UPA';

param.ph.c = 299792458; % In m/s
param.ph.freq = 10e9; 
param.ph.lambda = param.ph.c/param.ph.freq;


%%% target and interfering signal directions
phi=[-40];
phi_des = deg2rad(phi);
theta=[30];
theta_des = deg2rad(theta);

%%% DUT array configuration

P = 100;%100x100
DUT_3D.N = P;
DUT_3D.d = 0.015; 
dd = (0:1:(P-1))*DUT_3D.d;
DUT_3D.x = zeros(1,DUT_3D.N);
DUT_3D.y = dd;
DUT_3D.z = dd;

w = zeros(1,P^2);

% beforming权重计算
for i=1:length(theta_des)
    for j=1:length(phi_des)
        p_y=exp(-1j*2*pi/param.ph.lambda*cos(theta_des(i))*sin(phi_des(j))*DUT_3D.y);
        p_z=exp(-1j*2*pi/param.ph.lambda*sin(theta_des(i))*DUT_3D.z);
        a=kron(p_y, p_z);%1x10000
        w = w + a;
    end
end


%% 3D方向图
% 初始化3D方向图
pattern_3D = zeros(length(theta_range), length(phi_range));
% 3D方向图
for i=1:length(theta_range)
    for j=1:length(phi_range)
        p_y=exp(1j*2*pi/param.ph.lambda*cos(theta_range(i))*sin(phi_range(j))*DUT_3D.y);
        p_z=exp(1j*2*pi/param.ph.lambda*sin(theta_range(i))*DUT_3D.z);
        AF=kron(p_y, p_z);
        pattern_3D(i,j)=abs(w*AF.');
    end
end

% 归一化
pattern_3D = pattern_3D / max(pattern_3D(:));

% 转换为dB
pattern_3D_dB = 20*log10(pattern_3D);
%% 截面图
figure(1);
[~, a] = min(abs(theta_range - theta_des(1))); 
plot(rad2deg(phi_range),pattern_3D_dB(a,:),'r','LineWidth', 1.1);
grid on;
xlabel('方位角 [deg]','fontsize',14)
ylabel('归一化功率 [dB]','fontsize',14)
title(sprintf('截面波束方向图（仰角 %d° ）', theta(1)),'FontSize', 16);

figure(2);
[~, b] = min(abs(phi_range - phi_des(1))); 
plot(rad2deg(theta_range),pattern_3D_dB(:,b),'r','LineWidth', 1.1);
grid on;
xlabel('仰角 [deg]','fontsize',14)
ylabel('归一化功率 [dB]','fontsize',14)
title(sprintf('截面波束方向图（方位角 %d° ）',phi(1)),'FontSize', 16);
%% 3D绘图
figure(3);
[Phi, Theta] = meshgrid(rad2deg(phi_range), rad2deg(theta_range));
surf(Phi, Theta, pattern_3D_dB, 'EdgeColor', 'none');
xlabel('方位角 [deg]', 'FontSize', 12);
ylabel('仰角 [deg] ', 'FontSize', 12);
zlabel('归一化频率 [dB]', 'FontSize', 12);
title('UPA波束方向图', 'FontSize', 16);
colorbar;
axis tight;

考虑时延，多径

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
% UPA阵列参数
n_beishu=1;
Mx = 4*n_beishu;         % x方向阵元数 (方位角维度)
Mz = 8*n_beishu;         % z方向阵元数 (俯仰角维度)
fc = 3.5e9;     % 中心频率 (Hz)
c = 3e8;        % 光速
lambda = c/fc;  % 波长
d = lambda/2;   % 阵元间距

% 频域参数
B = 50e6;       % 带宽 (Hz)
N_freq = 1001;  % 频点数
frequencies = linspace(fc - B/2, fc + B/2, N_freq);

% 多径目标参数：[延迟(s), 俯仰角(°), 方位角(°), 功率(dB)]
paths = [
    160e-9,  -10, 30, 0;   % 第一条径
    860e-9,   0,   0, -5   % 第二条径
];
num_paths = size(paths, 1);
num_array = Mx*Mz;%天线个数

paths_power = 10.^(paths(:,4)/10);
paths_gain = sqrt(paths_power) ; 

%% 3D steering vector generate
% 阵列位置定义 (位于XOZ平面，中心在原点)
% x方向位置: [-1.5d, -0.5d, 0.5d, 1.5d]
% z方向位置: [-3.5d, -2.5d, ..., 3.5d]
x_pos = (-(Mx-1)/2 : (Mx-1)/2) * d;
z_pos = (-(Mz-1)/2 : (Mz-1)/2) * d;

[X, Z] = meshgrid(x_pos, z_pos);
% 生成网格位置 [X, Z] = meshgrid(x_pos, z_pos)
% 导向矢量数学模型:
% a(θ,φ) = exp(-j*(2π/λ)*[x*sinφ*cosθ + z*sinθ])
% 其中:
%   θ: 俯仰角 (从xoy平面起算，z+为90°, z-为-90°)
%   φ: 方位角 (y+为0°, x+为90°)
X_vec = reshape(X,[],1);
Z_vec = reshape(Z,[],1);
% figure
% scatter(X_vec,Z_vec)

theta_range = -90:1:90;    % 俯仰角范围 (°)
phi_range = -90:1:90;     % 方位角范围 (°)
N_theta = length(theta_range);
N_phi = length(phi_range);
spatial_steering_matrix =zeros(num_array,N_theta,N_phi);%每一个天线在相同方位角，俯仰角的相位调整
for theta_idx = 1: N_theta
    for phi_idx = 1: N_phi
        % 计算波矢量分量
        kx = sind(phi_range(phi_idx)) * cosd(theta_range(theta_idx));  % x方向波数分量
        ky = cosd(phi_range(phi_idx)) * cosd(theta_range(theta_idx));  % x方向波数分量
        kz = sind(theta_range(theta_idx));             % z方向波数分量
        for array_idx = 1: num_array
              x_pos_i = X_vec(array_idx);
              z_pos_i = Z_vec(array_idx);
            % spatial_steering_matrix(array_idx,theta_idx,phi_idx)=exp(-1j*(2*pi/lambda)*(x_pos_i*kx + z_pos_i*kz));
            spatial_steering_matrix(array_idx,theta_idx,phi_idx)=exp(-1j*(2*pi/lambda)*[x_pos_i,0,z_pos_i]*[kx;ky;kz ]);
        end
    end
end

%% 生成频域响应 (CFR)
H = zeros(num_array, N_freq); 

for p = 1:num_paths
    delay = paths(p,1);
    theta = (paths(p,2));  % 俯仰角 (弧度)
    phi = (paths(p,3));    % 方位角 (弧度)
    gain = paths_gain(p);
    [~,search_theta] = min(abs(theta-theta_range));
    [~,search_phi] = min(abs(phi-phi_range));
    steering_matrix = spatial_steering_matrix(:,search_theta,search_phi);
    % 添加时延和复增益
    phase_delay = exp(-1j*2*pi*frequencies*delay);
    H = H + gain * steering_matrix * phase_delay;
end

%% 转换到时域CIR
% 对每个天线进行IFFT
h = ifft(H,N_freq, 2);  %沿行做傅里叶变换
% 时延轴生成
delta_t = 1/(frequencies(end)-frequencies(1)); % 时延分辨率
tau = (0:N_freq-1)*delta_t;
tau_max = 1000e-9;  % 最大时延 (1000 ns),N_freq/B


% tau_max = 1/(frequencies(2)-frequencies(1));  % 最大时延 (1000 ns)
valid_idx = tau <= tau_max;
tau = tau(valid_idx);  % 截取有效时延范围
h = h(:,valid_idx);  % 截取有效时延响应
N_tau = length(tau);
%% Estimate
Power_estimated = zeros(N_theta, N_phi, N_tau); 
% for theta_idx = 1: N_theta
%     for phi_idx = 1: N_phi
%        Power_estimated(theta_idx,phi_idx,:)= h.' *conj(squeeze(spatial_steering_matrix(:,theta_idx,phi_idx)));
%     end
% end



for t = 1:N_tau
    for theta_idx = 1: N_theta
        for phi_idx = 1: N_phi
            snapshot = h(:,t);  % [num_array x 1]
            steering_vector =spatial_steering_matrix(:,theta_idx,phi_idx);
            Power_estimated(theta_idx,phi_idx,t) = abs(steering_vector' * snapshot)^2;
        end
    end
end
Power_estimated_dB = 10*log10(abs(Power_estimated));
%%
% 生成网格坐标（实际值）
[theta_mesh, phi_mesh, tau_mesh] = ndgrid(theta_range, phi_range, tau);

% 绘制散点图
figure;
scatter3(theta_mesh(:), phi_mesh(:), tau_mesh(:)*1e9, 30, Power_estimated_dB(:), 'filled');

% 设置坐标轴标签（带单位）
xlabel('Theta (度)'); 
ylabel('Phi (度)'); 
zlabel('时延 (ns)');
title('多径信道功率分布 (dB)');
colorbar;
colormap('parula'); % 或 'turbo'/'jet'
grid on;

% 标记峰值路径
[peak_power, peak_idx] = max(Power_estimated_dB(:));
[theta_idx, phi_idx, tau_idx] = ind2sub(size(Power_estimated_dB), peak_idx);
fprintf('主径: Theta=%.2f°, Phi=%.2f°, 时延=%.2f ns, 功率=%.2f dB\n', theta_range(theta_idx), phi_range(phi_idx), tau(tau_idx)*1e9, peak_power);


%%
% [~,tau_1_idx]= min(abs(paths(1,1)-tau));
% AAP_1 = Power_estimated_dB(:,:,tau_1_idx);
% figure;
% imagesc( phi_range ,theta_range, AAP_1);  % X:方位角 Y:俯仰角
% xlabel('Azimuth Angle \phi (°)');
% ylabel('Elevation Angle \theta (°)');
% title('Angle-Angle Power at Path 1 Delay');
% colorbar;
% axis xy;  % 保证角度方向正常（非图像坐标）
% [max_val, max_idx] = max(AAP_1(:));
% [theta_idx_max, phi_idx_max] = ind2sub(size(AAP_1), max_idx);
% theta_peak = theta_range(theta_idx_max);
% phi_peak = phi_range(phi_idx_max);
% disp(['主瓣方向为：Elevation=', num2str(theta_peak), '°, Azimuth=', num2str(phi_peak), '°']);
% 
% 
% [~,tau_2_idx]= min(abs(paths(2,1)-tau));
% AAP_2 = Power_estimated_dB(:,:,tau_2_idx);
% figure;
% imagesc( phi_range ,theta_range, AAP_2);  % X:方位角 Y:俯仰角
% xlabel('Azimuth Angle \phi (°)');
% ylabel('Elevation Angle \theta (°)');
% title('Angle-Angle Power at Path 2 Delay');
% colorbar;
% axis xy;  % 保证角度方向正常（非图像坐标）



%%
[~, tau_1_idx] = min(abs(paths(1,1) - tau));
AAP_1 = Power_estimated_dB(:,:,tau_1_idx);
AAP_1=AAP_1-max(max(AAP_1));
figure;
imagesc(phi_range, theta_range, AAP_1);  % X:方位角，Y:俯仰角
xlabel('Azimuth Angle \phi (°)');
ylabel('Elevation Angle \theta (°)');
title('Angle-Angle Power at Path 1 Delay [Simulation]');
colorbar;
clim([-40 0]);
axis xy;  % 保证图像角度方向一致

% 找主瓣
[max_val, max_idx] = max(AAP_1(:));
[theta_idx_max, phi_idx_max] = ind2sub(size(AAP_1), max_idx);
theta_peak = theta_range(theta_idx_max);
phi_peak = phi_range(phi_idx_max);
disp(['主瓣方向为：Elevation = ', num2str(theta_peak), '°, Azimuth = ', num2str(phi_peak), '°']);

% 添加红圈标记主瓣
hold on;
plot(phi_peak, theta_peak, 'ro', 'MarkerSize', 10, 'LineWidth', 2);

% 添加文字注释主瓣方向
text(phi_peak + 5, theta_peak, ['\leftarrow (', num2str(theta_peak), '°, ', num2str(phi_peak), '°)'], 'Color', 'red', 'FontSize', 12, 'FontWeight', 'bold');

考虑时延，多径，多普勒

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
clear;
clc;
% UPA阵列参数
n_beishu=1;
Mx = 4*n_beishu;         % x方向阵元数 (方位角维度)
Mz = 8*n_beishu;         % z方向阵元数 (俯仰角维度)
fc = 3.5e9;     % 中心频率 (Hz)
c = 3e8;        % 光速
lambda = c/fc;  % 波长
d = lambda/2;   % 阵元间距

% 时域参数
v=100;          %运动速度（m/s）
v_phi=20;       %运动方位角
v_theta=0;      %运动仰角
v_vector=[cosd(v_theta)*cosd(v_phi),cosd(v_theta)*sind(v_phi),sind(v_theta)];%运动方向矢量
f_max=v/lambda; %最大多普勒频移（Hz）
N_time=1002;    %时域点数
T_delta=1/(2*f_max);%时域采样间隔
doppler_delta=2*f_max/N_time;%多普勒域采样间隔
time=(0:N_time-1)*T_delta;

% 频域参数
B = 50e6;       % 带宽 (Hz)
N_freq = 1001;  % 频点数
frequencies = linspace(fc - B/2, fc + B/2, N_freq);%1x1001

% 多径目标参数：[延迟(s), 俯仰角(°), 方位角(°), 功率(dB)，多普勒频移(Hz)]
paths = [
    160e-9,  -10, 30, 0,  0;   % 第一条径
    860e-9,   0,   0, -5, 0    % 第二条径
];
num_paths = size(paths, 1);
for i=1:num_paths
    path_vector=[cosd(paths(i,2))*cosd(paths(i,3)),cosd(paths(i,2))*sind(paths(i,3)),sind(paths(i,2))];
    paths(i,5)=f_max*path_vector*v_vector.';
end
num_array = Mx*Mz;%天线个数
paths_power = 10.^(paths(:,4)/10);
paths_gain = sqrt(paths_power) ; 

%% 3D steering vector generate
% 阵列位置定义 (位于XOZ平面，中心在原点)
% x方向位置: [-1.5d, -0.5d, 0.5d, 1.5d]
% z方向位置: [-3.5d, -2.5d, ..., 3.5d]
x_pos = (-(Mx-1)/2 : (Mx-1)/2) * d;
z_pos = (-(Mz-1)/2 : (Mz-1)/2) * d;

[X, Z] = meshgrid(x_pos, z_pos);
% 生成网格位置 [X, Z] = meshgrid(x_pos, z_pos)
% 导向矢量数学模型:
% a(θ,φ) = exp(-j*(2π/λ)*[x*sinφ*cosθ + z*sinθ])
% 其中:
%   θ: 俯仰角 (从xoy平面起算，z+为90°, z-为-90°)
%   φ: 方位角 (y+为0°, x+为90°)
X_vec = reshape(X,[],1);
Z_vec = reshape(Z,[],1);
% figure
% scatter(X_vec,Z_vec)

theta_range = -90:1:90;    % 俯仰角范围 (°)
phi_range = -90:1:90;     % 方位角范围 (°)
N_theta = length(theta_range);
N_phi = length(phi_range);
spatial_steering_matrix =zeros(num_array,N_theta,N_phi);%每一个天线在相同方位角，俯仰角的相位调整
for theta_idx = 1: N_theta
    for phi_idx = 1: N_phi
        % 计算波矢量分量
        kx = sind(phi_range(phi_idx)) * cosd(theta_range(theta_idx));  % x方向波数分量
        ky = cosd(phi_range(phi_idx)) * cosd(theta_range(theta_idx));  % x方向波数分量
        kz = sind(theta_range(theta_idx));             % z方向波数分量
        for array_idx = 1: num_array
            x_pos_i = X_vec(array_idx);
            z_pos_i = Z_vec(array_idx);
            % spatial_steering_matrix(array_idx,theta_idx,phi_idx)=exp(-1j*(2*pi/lambda)*(x_pos_i*kx + z_pos_i*kz));

            spatial_steering_matrix(array_idx,theta_idx,phi_idx)=exp(-1j*(2*pi/lambda)*[x_pos_i,0,z_pos_i]*[kx;ky;kz ]);

        end
    end
end

%% 生成频域响应 (CFR),h(t,f)
H = zeros(num_array, N_freq,N_time); 

for p = 1:num_paths
    delay = paths(p,1);
    theta = (paths(p,2));  % 俯仰角 (弧度)
    phi = (paths(p,3));    % 方位角 (弧度)
    gain = paths_gain(p);
    [~,search_theta] = min(abs(theta-theta_range));
    [~,search_phi] = min(abs(phi-phi_range));
    steering_matrix = spatial_steering_matrix(:,search_theta,search_phi);%num_arrayx1
    % 添加时延和复增益
    phase_delay = exp(-1j*2*pi*frequencies*delay);%1x1001
    H_1=gain * steering_matrix * phase_delay;%num_arrayx1001
    % 添加时间轴
    phase_dopppler=exp(1j*2*pi*paths(p,5)*time);
    doppler_reshape = reshape(phase_dopppler, 1, 1, []);%1x1x501
    H_2 = bsxfun(@times, H_1, doppler_reshape);%num_arrayx1001x501
    H = H +H_2;
end

%% 转换到时域CIR，h(t,τ)
% 对每个天线进行IFFT
h = ifft(H,N_freq, 2);  %沿行做傅里叶变换
% 时延轴生成
delta_t = 1/(frequencies(end)-frequencies(1)); % 时延分辨率
tau = (0:N_freq-1)*delta_t;
tau_max = 1000e-9;  % 最大时延 (1000 ns),N_freq/B

% tau_max = 1/(frequencies(2)-frequencies(1));  % 最大时延 (1000 ns)
valid_idx = tau <= tau_max;
tau = tau(valid_idx);  % 截取有效时延范围
h = h(:,valid_idx,:);  % 截取有效时延响应
N_tau = length(tau);

%% 转换到s(v,τ)
s=fft(h,N_time,3);%-fm~fm
s=fftshift(s,3);
%多普勒频谱轴生成
doppler=(-N_time/2:N_time/2-1)*doppler_delta;
doppler_max=max(paths(:,5))+10;%截取与多径的多普勒频移相关的一部分
doppler_min=min(paths(:,5))-10;
doppler_idx=find(doppler>=(doppler_min)&doppler<=(doppler_max));
doppler=doppler(doppler_idx);
s=s(:,:,doppler_idx);
N_doppler=length(doppler);



%% Estimate
Power_estimated = zeros(N_theta, N_phi, N_tau,N_doppler); 
% for theta_idx = 1: N_theta
%     for phi_idx = 1: N_phi
%        Power_estimated(theta_idx,phi_idx,:)= h.' *conj(squeeze(spatial_steering_matrix(:,theta_idx,phi_idx)));
%     end
% end


for f=1:N_doppler
    for t = 1:N_tau
        for theta_idx = 1: N_theta
            for phi_idx = 1: N_phi
                snapshot = squeeze(s(:,t,f));  % [num_array x 1]
                steering_vector =spatial_steering_matrix(:,theta_idx,phi_idx);
                Power_estimated(theta_idx,phi_idx,t,f) = abs(steering_vector' * snapshot)^2;%共轭转置
            end
        end
    end
end
Power_estimated_dB = 10*log10(abs(Power_estimated));
% 标记峰值路径
[peak_power, peak_idx] = max(Power_estimated_dB(:));
[theta_idx, phi_idx, tau_idx,doppler_idx] = ind2sub(size(Power_estimated_dB), peak_idx);
fprintf('主径: Theta=%.2f°, Phi=%.2f°, 时延=%.2f ns, 多普勒频移=%.2f Hz, 功率=%.2f dB\n', theta_range(theta_idx), phi_range(phi_idx), tau(tau_idx)*1e9,doppler(doppler_idx), peak_power);