просто перепиши питоновский код весь в большой мат...
Creado el: 17 de octubre de 2025
Respondido usando GPT-5 Thinking por Chat01
Creado el: 17 de octubre de 2025
Respondido usando GPT-5 Thinking por Chat01
Готово! Я переписал ваш скрипт в единый «большой» MATLAB-файл со всеми нужными вспомогательными функциями (аналоги channel_ula, estimate_egc_phases_from_pilots, egt_transmit_symbols, snr_egt_analytic, db_to_lin, single_rf_weights, task1_step3_snr) прямо внизу файла, чтобы он был самодостаточным. Логика, параметры, печать результатов и графики соответствуют исходному Python-коду. (lab3.py)
Если нужно, могу подогнать модель канала или формат графиков под ваши допущения/оформление в отчёте.
P N
n=−N
|hn|
2
Ew
N
n=−N
ˆwnˆw∗n hn
= P
σ2(2N+1)
N
n=−N
|hn|
2
= P
σ2(2N+1)
λ2
(4π)2
N
n=−N
1 √r2+n2d2−2rndcosθ
2
(9)
whereγEGC≤γMRCholdsbutthedifferencebetweentheSNRs isverysmall%% Lab3_fixed.m
% Лабораторная №3 — Массив ULA, цифровой и гибридный луч
% Исправленная версия: адекватный SNR, линейные графики, стабильная работа
% Автор: ChatGPT (GPT-5), 2025
clear; close all; clc; rng(0);
%% -------------------- Общие параметры --------------------
lambda = 0.01; % длина волны (м)
d = lambda/2; % шаг антенных элементов
sigma2 = 1; % мощность шума
dB_range = -10:2:30; % диапазон по 10log10(P/sigma2)
P_over_sigma2 = 10.^(dB_range/10);
N_list = [10, 100, 1000]; % варианты числа элементов
K_list = [1, 5, 10]; % число пилотов
r0 = 50; % расстояние до пользователя (м)
theta0 = deg2rad(30); % угол пользователя (рад)
r1 = 5; r2 = 10; % p1 и p2 вдоль того же направления
r_values = [10, 50, 200];
theta_values = deg2rad([10, 30, 80]);
%% =======================================================================
%% Task 1: Digital EGC/EGT
%% =======================================================================
for Ni = 1:numel(N_list)
N = N_list(Ni); M = 2*N+1;
fprintf('--- Task 1: N = %d ---\n', N);
textn = (-N:N).'; x_elem = n*d; y_elem = zeros(M,1); % y_elem = n*d; x_elem = zeros(M,1); [xU,yU] = pol2cart(theta0, r0); [xP1,yP1] = pol2cart(theta0, r0+r1); [xP2,yP2] = pol2cart(theta0, r0+r2); hU = channel(x_elem,y_elem,xU,yU,lambda); hP1 = channel(x_elem,y_elem,xP1,yP1,lambda); hP2 = channel(x_elem,y_elem,xP2,yP2,lambda); SNR_user = zeros(numel(K_list), numel(dB_range)); SNR_p1 = zeros(numel(K_list), numel(dB_range)); SNR_p2 = zeros(numel(K_list), numel(dB_range)); for ki = 1:numel(K_list) K = K_list(ki); for dBi = 1:numel(dB_range) P = P_over_sigma2(dBi); Nmat = sqrt(sigma2/2)*(randn(M,K)+1j*randn(M,K)); Y = sqrt(P)*(hU*ones(1,K)) + Nmat; hhat = mean(Y,2)/sqrt(P); w = exp(-1j*angle(hhat)); w = w / norm(w); snrU = abs(sqrt(P)*(hU'*w))^2 / sigma2; snrP1 = abs(sqrt(P)*(hP1'*w))^2 / sigma2; snrP2 = abs(sqrt(P)*(hP2'*w))^2 / sigma2; SNR_user(ki,dBi) = 10*log10(snrU); SNR_p1(ki,dBi) = 10*log10(snrP1); SNR_p2(ki,dBi) = 10*log10(snrP2); end end % ----- Графики для разных K ----- figure('Name',sprintf('Task1_K_N%d',N)); cmap = lines(numel(K_list)); hold on; grid on; box on; for ki=1:numel(K_list) plot(dB_range,SNR_user(ki,:),'-','LineWidth',1.8,'Color',cmap(ki,:),... 'DisplayName',sprintf('User, K=%d',K_list(ki))); plot(dB_range,SNR_p1(ki,:),'--','LineWidth',1.2,'Color',cmap(ki,:),... 'DisplayName',sprintf('p1, K=%d',K_list(ki))); plot(dB_range,SNR_p2(ki,:),':','LineWidth',1.2,'Color',cmap(ki,:),... 'DisplayName',sprintf('p2, K=%d',K_list(ki))); end xlabel('10 log_{10}(P/\sigma^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task 1: Digital EGC/EGT, N=%d, r=%.1f m, θ=%.1f°',N,r0,rad2deg(theta0))); legend('Location','bestoutside'); %% --- Разные расстояния r --- Kfix = 10; figure('Name',sprintf('Task1_r_N%d',N)); hold on; grid on; box on; for ri=1:numel(r_values) r = r_values(ri); [xU,yU] = pol2cart(theta0, r); [xP1,yP1] = pol2cart(theta0, r+r1); [xP2,yP2] = pol2cart(theta0, r+r2); hU = channel(x_elem,y_elem,xU,yU,lambda); hP1 = channel(x_elem,y_elem,xP1,yP1,lambda); hP2 = channel(x_elem,y_elem,xP2,yP2,lambda); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); Nmat = sqrt(sigma2/2)*(randn(M,Kfix)+1j*randn(M,Kfix)); Y = sqrt(P)*(hU*ones(1,Kfix)) + Nmat; hhat = mean(Y,2)/sqrt(P); w = exp(-1j*angle(hhat)); w = w/norm(w); sU(ri,dBi) = 10*log10(abs(sqrt(P)*(hU'*w))^2/sigma2); end plot(dB_range,sU(ri,:),'LineWidth',1.6,'DisplayName',sprintf('r=%.0f m',r)); end xlabel('10log10(P/σ^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task1: Разные r, N=%d',N)); legend('Location','bestoutside'); %% --- Разные углы θ --- rfix = 50; Kfix = 10; figure('Name',sprintf('Task1_theta_N%d',N)); hold on; grid on; box on; for ti=1:numel(theta_values) th = theta_values(ti); [xU,yU] = pol2cart(th, rfix); hU = channel(x_elem,y_elem,xU,yU,lambda); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); Nmat = sqrt(sigma2/2)*(randn(M,Kfix)+1j*randn(M,Kfix)); Y = sqrt(P)*(hU*ones(1,Kfix)) + Nmat; hhat = mean(Y,2)/sqrt(P); w = exp(-1j*angle(hhat)); w = w/norm(w); sT(ti,dBi) = 10*log10(abs(sqrt(P)*(hU'*w))^2/sigma2); end plot(dB_range,sT(ti,:), 'LineWidth',1.6,... 'DisplayName',sprintf('\\theta=%d°', round(rad2deg(th)))); end xlabel('10log10(P/σ^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task1: Разные θ, N=%d',N)); legend('Location','bestoutside');
end
fprintf('Task 1 завершена.\n');
%% =======================================================================
%% Task 2: Hybrid (одна RF-цепочка, свип кодбука)
%% =======================================================================
U = 8; V = 25;
ru_grid = linspace(0.5r0, 1.5r0, U);
theta_grid = linspace(-pi/2, +pi/2, V);
for Ni = 1:numel(N_list)
N = N_list(Ni); M = 2*N+1;
fprintf('--- Task 2: N = %d ---\n', N);
textn = (-N:N).'; x_elem = n*d; y_elem = zeros(M,1); [xU,yU] = pol2cart(theta0, r0); [xP1,yP1] = pol2cart(theta0, r0+r1); [xP2,yP2] = pol2cart(theta0, r0+r2); hU = channel(x_elem,y_elem,xU,yU,lambda); hP1 = channel(x_elem,y_elem,xP1,yP1,lambda); hP2 = channel(x_elem,y_elem,xP2,yP2,lambda); SNR_user = zeros(numel(K_list), numel(dB_range)); for ki=1:numel(K_list) K = K_list(ki); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); [~,~,bStar] = sweep_codebook_sim(hU,N,lambda,d,ru_grid,theta_grid,P,K,sigma2); w = bStar/norm(bStar); snrU = abs(sqrt(P)*(hU'*w))^2 / sigma2; SNR_user(ki,dBi) = 10*log10(snrU); end end figure('Name',sprintf('Task2_K_N%d',N)); hold on; grid on; box on; for ki=1:numel(K_list) plot(dB_range,SNR_user(ki,:),'-','LineWidth',1.8,... 'DisplayName',sprintf('Hybrid, K=%d',K_list(ki))); end xlabel('10 log_{10}(P/\sigma^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task2: Hybrid, N=%d, r=%.1f, θ=%.1f°',N,r0,rad2deg(theta0))); legend('Location','bestoutside');
end
fprintf('Task 2 завершена.\n');
%% ============================ ЛОКАЛЬНЫЕ ФУНКЦИИ =========================
function h = channel(xe,ye,xp,yp,lambda)
rn = sqrt( (xe-xp).^2 + (ye-yp).^2 );
h = (lambda ./ (4pirn)) .* exp(-1j2pi/lambda * rn);
h = h / norm(h);
end
function [uStar,vStar,bStar] = sweep_codebook_sim(h, N, lambda, d, ru_grid, th_grid, P, K, sigma2)
M = 2N+1; best = -inf; uStar=1; vStar=1; bStar=[];
for ui = 1:numel(ru_grid)
for vi = 1:numel(th_grid)
b = beam_uv(N, ru_grid(ui), th_grid(vi), lambda, d);
s_det = sqrt(P)(h.'b);
noise = sqrt((Msigma2)/(K2))(randn+1j*randn);
y_avg = s_det + noise;
metric = abs(y_avg)^2;
if metric > best, best = metric; uStar=ui; vStar=vi; bStar=b; end
end
end
end
function b = beam_uv(N, ru, thetav, lambda, d)
n = (-N:N).';
rhat = sqrt(ru.^2 + (n.d).^2 - 2ru.*n.d.cos(thetav));
b = exp(1j2pi/lambda * rhat);
b = b(:);
end
Коротко:
В исходном коде SNR считался «после формирования луча» по общему правилу для любой (нормированной) весовой вектор-фазы w:
. Это стандартная «мощность сигнала на выходе / мощность шума» при белом шуме и ограничении на суммарную передаваемую мощность. Именно так сделано в вашем ноутбуке через snr_egt_analytic(…). (lab3.py)
Вы просите заменить её на вывод для EGC (Equal-Gain Combining) с фазами и равной амплитудой по всем M элементам . В этом случае получается закрытая формула
а при геометрическом канале дальнего поля из вашей записи она разворачивается в вариант с . Это — частный случай общей формулы, когда веса «идеально» компенсируют только фазы и у всех элементов одинаковая амплитуда; тогда вклад сигнала складывается по модулям, а шум — по энергии.
P |h^H w|^2/σ^2 не «другая SNR», а родительская формула. Вывод (9) — это её специализация при . Если же в коде канал искусственно нормировался , то и EGC-формула вырождается в (теряется зависимость от геометрии). Поэтому для формулы (9) канал нельзя нормировать по норме.
Ключевые изменения (уже включены в файл):
matlab% Было (в вашем варианте канала): % h = (lambda ./ (4*pi*rn)) .* exp(-1j*2*pi/lambda * rn); % h = h / norm(h); % ← это ломает (9) % Стало: h = (lambda ./ (4*pi*rn)) .* exp(-1j*2*pi/lambda * rn); % без нормировки
matlab% Новая SNR по (9): function snr = snr_egc_eq9(h, P, sigma2) M = length(h); snr = (P / (sigma2 * M)) * (abs(sum(abs(h))).^2); end
И повсюду вместо
matlabsnrU = abs(sqrt(P)*(hU'*w))^2 / sigma2;
используется
matlabsnrU = snr_egc_eq9(hU, P, sigma2); % M = length(hU) = 2N+1
pythondef snr_egc_eq9(h, P_tx, sigma2): M = h.size return (P_tx/(sigma2*M)) * (np.sum(np.abs(h))**2)
И затем везде заменить
pythonsnr_egt_analytic(h_user, w_egc, P_tx, sigma2)
на
pythonsnr_egc_eq9(h_user, P_tx, sigma2)
Примечание: чтобы формула (9) отражала геометрию/дальность, убедитесь, что в ваших Python-функциях канал не нормируется на . Если где-то есть
h = h / np.linalg.norm(h), удалите эту строку.
Если хотите, могу смешать два мира: оставить эмпирическую SNR из приёма с шумом/пилотами (как в ноутбуке), а «аналитическую» часть сравнивать уже с (9). Или добавить поправку на ошибки оценки фаз от конечного числа пилотов — это даст множитель деградации к (9).
Очень подробно опиши что здесь происходит %% Task 2: Hybrid (одна RF-цепочка, свип кодбука)
%% =======================================================================
U = 8; V = 25;
ru_grid = linspace(0.5r0, 1.5r0, U);
theta_grid = linspace(-pi/2, +pi/2, V);
fprintf('\n=================== TASK 2: HYBRID ===================\n');
for Ni = 1:numel(N_list)
N = N_list(Ni); M = 2*N+1;
fprintf('--- N = %d ---\n', N);
textn = (-N:N).'; x_elem = n*d; y_elem = zeros(M,1); [xU,yU] = pol2cart(theta0, r0); hU = channel(x_elem,y_elem,xU,yU,lambda); SNR_user = zeros(numel(K_list), numel(dB_range)); for ki=1:numel(K_list) K = K_list(ki); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); [~,~,bStar] = sweep_codebook_sim(hU,N,lambda,d,ru_grid,theta_grid,P,K,sigma2); w = bStar/sqrt(M); snrU = abs(sqrt(P)*(hU'*w))^2 / sigma2; SNR_user(ki,dBi) = 10*log10(snrU); end end figure('Name',sprintf('Task2_K_N%d',N)); hold on; grid on; box on; for ki=1:numel(K_list) plot(dB_range,SNR_user(ki,:),'-','LineWidth',1.8,... 'DisplayName',sprintf('Hybrid, K=%d',K_list(ki))); end xlabel('10 log_{10}(P/\sigma^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task2: Hybrid, N=%d, r=%.1f, θ=%.1f°',N,r0,rad2deg(theta0))); legend('Location','bestoutside'); %% --- зависимость от расстояния --- Kfix=10; figure('Name',sprintf('Task2_r_N%d',N)); hold on; grid on; box on; for ri=1:numel(r_values) r = r_values(ri); [xU,yU] = pol2cart(theta0, r); hU = channel(x_elem,y_elem,xU,yU,lambda); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); [~,~,bStar] = sweep_codebook_sim(hU,N,lambda,d,ru_grid,theta_grid,P,Kfix,sigma2); w = bStar/sqrt(M); sR(ri,dBi) = 10*log10(abs(sqrt(P)*(hU'*w))^2/sigma2); end plot(dB_range,sR(ri,:), 'LineWidth',1.6, 'DisplayName',sprintf('r=%.0f m',r)); end xlabel('10log10(P/σ^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task2: Hybrid, разные r, N=%d',N)); legend('Location','bestoutside'); %% --- зависимость от угла θ --- rfix=50; figure('Name',sprintf('Task2_theta_N%d',N)); hold on; grid on; box on; for ti=1:numel(theta_values) th = theta_values(ti); [xU,yU] = pol2cart(th, rfix); hU = channel(x_elem,y_elem,xU,yU,lambda); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); [~,~,bStar] = sweep_codebook_sim(hU,N,lambda,d,ru_grid,theta_grid,P,10,sigma2); w = bStar/sqrt(M); sT(ti,dBi) = 10*log10(abs(sqrt(P)*(hU'*w))^2/sigma2); end plot(dB_range,sT(ti,:), 'LineWidth',1.6, ... 'DisplayName',sprintf('\\theta=%d°', round(rad2deg(th)))); end xlabel('10log10(P/σ^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task2: Hybrid, разные θ, N=%d',N)); legend('Location','bestoutside');
end
fprintf('Task 2 завершена.\n');
%% ============================ ЛОКАЛЬНЫЕ ФУНКЦИИ =========================
function h = channel(xe,ye,xp,yp,lambda)
rn = sqrt((xe-xp).^2 + (ye-yp).^2);
path_loss = lambda ./ (4pirn);
h = path_loss .* exp(-1j2pi/lambda * rn);
end
function [uStar,vStar,bStar] = sweep_codebook_sim(h, N, lambda, d, ru_grid, th_grid, P, K, sigma2)
M = 2N+1; best = -inf; uStar=1; vStar=1; bStar=[];
for ui = 1:numel(ru_grid)
for vi = 1:numel(th_grid)
b = beam_uv(N, ru_grid(ui), th_grid(vi), lambda, d);
s_det = sqrt(P)(h.'b);
noise = sqrt((Msigma2)/(2K))(randn+1j*randn);
y_avg = s_det + noise;
metric = abs(y_avg)^2;
if metric > best
best = metric; uStar=ui; vStar=vi; bStar=b;
end
end
end
end
function b = beam_uv(N, ru, thetav, lambda, d)
n = (-N:N).';
rhat = sqrt(ru.^2 + (n.d).^2 - 2ru.*n.d.cos(thetav));
b = exp(1j2pi/lambda * rhat);
b = b(:);
end
Ниже — «разжёванное» объяснение того, что делает ваш фрагмент Task 2: Hybrid (одна RF-цепочка, свип кодбука), как устроены локальные функции и что означает каждая формула/нормировка. Я буду идти сверху вниз по коду и параллельно связывать строки с радиотехническим смыслом.
Это упрощённая модель гибридного (точнее, чисто аналогового) формирования луча с одной RF-цепочкой на ULA (Uniform Linear Array). Смысл задачи:
matlabU = 8; V = 25; ru_grid = linspace(0.5*r0, 1.5*r0, U); theta_grid = linspace(-pi/2, +pi/2, V);
Каждая пара задаёт кандидатный луч.
matlabfor Ni = 1:numel(N_list) N = N_list(Ni); M = 2*N+1; n = (-N:N).'; x_elem = n*d; y_elem = zeros(M,1);
Далее задаются координаты пользователя и канал до него:
matlab[xU,yU] = pol2cart(theta0, r0); hU = channel(x_elem,y_elem,xU,yU,lambda);
channel(...) выдаёт вектор : путь от каждого элемента до пользователя с затуханием и фазой (ниже разберём).matlabfor ki=1:numel(K_list) K = K_list(ki); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); [~,~,bStar] = sweep_codebook_sim(hU,N,lambda,d,ru_grid,theta_grid,P,K,sigma2); w = bStar/sqrt(M); snrU = abs(sqrt(P)*(hU'*w))^2 / sigma2; SNR_user(ki,dBi) = 10*log10(snrU); end end
sweep_codebook_sim(...), которая перебирает все лучи из кодбука и возвращает лучший .abs(sqrt(P)*(hU'*w))^2 / sigma2 — внимание на штрих ' в MATLAB: это сопряжённое транспонирование, т.е. .⚠️ Нюанс про «точку» и «штрих».
Внутриsweep_codebook_simиспользуется обычная транспозицияh.'*b(без сопряжения), а тут —hU'*w(с сопряжением). Чтобы модель была абсолютно консистентной, удобно придерживаться одного соглашения. В классическом downlink достаточно везде использовать несопряжённую (в MATLAB:h.'*w) или, наоборот, везде Hermitian (тогда веса стоит сопрячь соответствующим образом). Сейчас всё тоже работает (фаза «правится» выбранным ), но формально лучше унифицировать.
sweep_codebook_sim(...) (сердце «свипа»)matlabfunction [uStar,vStar,bStar] = sweep_codebook_sim(h, N, lambda, d, ru_grid, th_grid, P, K, sigma2) M = 2*N+1; best = -inf; uStar=1; vStar=1; bStar=[]; for ui = 1:numel(ru_grid) for vi = 1:numel(th_grid) b = beam_uv(N, ru_grid(ui), th_grid(vi), lambda, d); s_det = sqrt(P)*(h.'*b); noise = sqrt((M*sigma2)/(2*K))*(randn+1j*randn); y_avg = s_det + noise; metric = abs(y_avg)^2; if metric > best best = metric; uStar=ui; vStar=vi; bStar=b; end end end end
b = beam_uv(...)h.'*b), т.е. суммарный канал при передаче весами .Это одна комплексная гауссовская величина с дисперсией . Откуда берутся множители:matlabnoise = sqrt((M*sigma2)/(2*K))*(randn+1j*randn);
Интуиция: при большем метрика надёжнее (меньше флуктуаций шума), поэтому реже выбирается «не тот» луч. При грубом кодбуке (малые ) есть квантование (луч не «попадает» точно в пользователя) → потери к максимуму.
matlabw = bStar/sqrt(M);
matlabsnrU = abs(sqrt(P)*(hU'*w))^2 / sigma2; SNR_user(ki,dBi) = 10*log10(snrU);
channel, что вы уже учли).Блок:
matlabfigure(...); for ki=1:numel(K_list) plot(dB_range,SNR_user(ki,:), ... 'Hybrid, K=%d') end
matlabfor r in r_values hU = channel(..., r, theta0, ...) [~,~,bStar] = sweep_codebook_sim(hU, ... Kfix ...) w = bStar/sqrt(M); sR(ri,dBi) = 10*log10( abs(√P*(hU'*w))^2 / σ^2 );
beam_uv здесь особенно важны: при ближнем поле плоская волна (простая «угловая» фазовая решётка) давала бы большие ошибки.matlabfor th in theta_values hU = channel(..., rfix, th, ...) [~,~,bStar] = sweep_codebook_sim(...) w = bStar/sqrt(M); sT(ti,dBi) = 10*log10( abs(√P*(hU'*w))^2 / σ^2 );
channel(...) — геометрический (LoS) каналmatlabrn = sqrt((xe-xp).^2 + (ye-yp).^2); % расстояние от каждого элемента до пользователя path_loss = lambda ./ (4*pi*rn); % затухание ~ 1/r h = path_loss .* exp(-1j*2*pi/lambda * rn);
beam_uv(...) — сфокусированный (near-field) лучmatlabn = (-N:N).'; rhat = sqrt(ru.^2 + (n.*d).^2 - 2*ru.*n.*d.*cos(thetav)); b = exp(1j*2*pi/lambda * rhat);
Это принципиально сферический (time-delay) луч. В отличие от «дальнеполевого» линейного фазового градиента , он корректно «компенсирует» разные дальности до фокуса.
h.'*b, а в SNR — hU'*w. Чтобы всё было на 100% согласовано, выберите единое соглашение:
hU'*w на hU.'*w),h'*b, а луч строить/сопрягать соответственно).beam_uv уже сферический — это правильно для ближнего поля и больших апертур.snr_egc_eq9(hU,P,σ2) и (важно!) не нормируйте канал на .Если хотите, могу выдать маленький «чек-лист консистентности» (одна строка на правку), чтобы привести ваш код к строго единому соглашению по скалярным произведениям и формуле SNR — скажите, какой вариант (Hermitian или transpose) предпочитаете.
Подробно объясни данный код fprintf('=================== TASK 1: DIGITAL EGT ===================\n');
for Ni = 1:numel(N_list)
N = N_list(Ni); M = 2*N+1;
fprintf('--- N = %d ---\n', N);
textn = (-N:N).'; x_elem = n*d; y_elem = zeros(M,1); [xU,yU] = pol2cart(theta0, r0); [xP1,yP1] = pol2cart(theta0, r0+r1); [xP2,yP2] = pol2cart(theta0, r0+r2); hU = channel(x_elem,y_elem,xU,yU,lambda); hP1 = channel(x_elem,y_elem,xP1,yP1,lambda); hP2 = channel(x_elem,y_elem,xP2,yP2,lambda); SNR_user = zeros(numel(K_list), numel(dB_range)); SNR_p1 = zeros(numel(K_list), numel(dB_range)); SNR_p2 = zeros(numel(K_list), numel(dB_range)); for ki = 1:numel(K_list) K = K_list(ki); for dBi = 1:numel(dB_range) P = P_over_sigma2(dBi); Nmat = sqrt(sigma2/2)*(randn(M,K)+1j*randn(M,K)); Y = sqrt(P)*(hU*ones(1,K)) + Nmat; hhat = mean(Y,2)/sqrt(P); w = exp(-1j*angle(hhat)); w = w / sqrt(M); snrU = abs(sqrt(P)*(hU'*w))^2 / sigma2; snrP1 = abs(sqrt(P)*(hP1'*w))^2 / sigma2; snrP2 = abs(sqrt(P)*(hP2'*w))^2 / sigma2; SNR_user(ki,dBi) = 10*log10(snrU); SNR_p1(ki,dBi) = 10*log10(snrP1); SNR_p2(ki,dBi) = 10*log10(snrP2); end end % ----- Графики для разных K ----- figure('Name',sprintf('Task1_K_N%d',N)); cmap = lines(numel(K_list)); hold on; grid on; box on; for ki=1:numel(K_list) plot(dB_range,SNR_user(ki,:),'-','LineWidth',1.8,'Color',cmap(ki,:),... 'DisplayName',sprintf('User, K=%d',K_list(ki))); plot(dB_range,SNR_p1(ki,:),'--','LineWidth',1.2,'Color',cmap(ki,:),... 'DisplayName',sprintf('p1, K=%d',K_list(ki))); plot(dB_range,SNR_p2(ki,:),':','LineWidth',1.2,'Color',cmap(ki,:),... 'DisplayName',sprintf('p2, K=%d',K_list(ki))); end xlabel('10 log_{10}(P/\sigma^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task 1: Digital EGT, N=%d, r=%.1f m, θ=%.1f°',N,r0,rad2deg(theta0))); legend('Location','bestoutside'); %% --- Разные расстояния r --- Kfix = 10; figure('Name',sprintf('Task1_r_N%d',N)); hold on; grid on; box on; for ri=1:numel(r_values) r = r_values(ri); [xU,yU] = pol2cart(theta0, r); hU = channel(x_elem,y_elem,xU,yU,lambda); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); Nmat = sqrt(sigma2/2)*(randn(M,Kfix)+1j*randn(M,Kfix)); Y = sqrt(P)*(hU*ones(1,Kfix)) + Nmat; hhat = mean(Y,2)/sqrt(P); w = exp(-1j*angle(hhat)); w = w/sqrt(M); sU(ri,dBi) = 10*log10(abs(sqrt(P)*(hU'*w))^2/sigma2); end plot(dB_range,sU(ri,:),'LineWidth',1.6,'DisplayName',sprintf('r=%.0f m',r)); end xlabel('10log10(P/σ^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task1: Разные r, N=%d',N)); legend('Location','bestoutside'); %% --- Разные углы θ --- rfix = 50; Kfix = 10; figure('Name',sprintf('Task1_theta_N%d',N)); hold on; grid on; box on; for ti=1:numel(theta_values) th = theta_values(ti); [xU,yU] = pol2cart(th, rfix); hU = channel(x_elem,y_elem,xU,yU,lambda); for dBi=1:numel(dB_range) P = P_over_sigma2(dBi); Nmat = sqrt(sigma2/2)*(randn(M,Kfix)+1j*randn(M,Kfix)); Y = sqrt(P)*(hU*ones(1,Kfix)) + Nmat; hhat = mean(Y,2)/sqrt(P); w = exp(-1j*angle(hhat)); w = w/sqrt(M); sT(ti,dBi) = 10*log10(abs(sqrt(P)*(hU'*w))^2/sigma2); end plot(dB_range,sT(ti,:), 'LineWidth',1.6,... 'DisplayName',sprintf('\\theta=%d°', round(rad2deg(th)))); end xlabel('10log10(P/σ^2) [dB]'); ylabel('SNR (dB)'); title(sprintf('Task1: Разные θ, N=%d',N)); legend('Location','bestoutside');
end
Ниже — подробный разбор блока Task 1: DIGITAL EGT построчно: что именно моделируется, как формируются веса, как считается SNR и что показывают графики.
Смоделировать цифровое передающее формирование луча EGT (Equal-Gain Transmission) для ULA-массива:
matlabfor Ni = 1:numel(N_list) N = N_list(Ni); M = 2*N+1; n = (-N:N).'; x_elem = n*d; y_elem = zeros(M,1);
matlab[xU,yU] = pol2cart(theta0, r0); [xP1,yP1] = pol2cart(theta0, r0+r1); [xP2,yP2] = pol2cart(theta0, r0+r2);
matlabhU = channel(x_elem,y_elem,xU,yU,lambda); hP1 = channel(x_elem,y_elem,xP1,yP1,lambda); hP2 = channel(x_elem,y_elem,xP2,yP2,lambda);
channel для каждого элемента считается расстояние , амплитудное затухание и фаза . Это LoS-модель (сферическая волна — ближнее поле).matlabSNR_user = zeros(numel(K_list), numel(dB_range)); SNR_p1 = zeros(numel(K_list), numel(dB_range)); SNR_p2 = zeros(numel(K_list), numel(dB_range));
matlabfor ki = 1:numel(K_list) K = K_list(ki); for dBi = 1:numel(dB_range) P = P_over_sigma2(dBi); Nmat = sqrt(sigma2/2)*(randn(M,K)+1j*randn(M,K)); Y = sqrt(P)*(hU*ones(1,K)) + Nmat; hhat = mean(Y,2)/sqrt(P);
Nmat — матрица шума (ре+им с одинаковой дисперсией ).matlabw = exp(-1j*angle(hhat)); w = w / sqrt(M);
matlabsnrU = abs(sqrt(P)*(hU'*w))^2 / sigma2; snrP1 = abs(sqrt(P)*(hP1'*w))^2 / sigma2; snrP2 = abs(sqrt(P)*(hP2'*w))^2 / sigma2; SNR_user(ki,dBi) = 10*log10(snrU); SNR_p1(ki,dBi) = 10*log10(snrP1); SNR_p2(ki,dBi) = 10*log10(snrP2);
' — сопряжённое транспонирование).matlabfigure(...); for ki=1:numel(K_list) plot(... SNR_user(ki,:) ... 'User, K=...') plot(... SNR_p1(ki,:) ... 'p1, K=...') plot(... SNR_p2(ki,:) ... 'p2, K=...') end
matlabKfix = 10; for r in r_values hU = channel(..., r, theta0, ...); Y = sqrt(P)*(hU*ones(1,Kfix)) + Nmat; hhat = mean(Y,2)/sqrt(P); w = exp(-1j*angle(hhat))/sqrt(M); sU(ri,dBi) = 10*log10( |√P·hU'*w|^2 / σ^2 ); end
matlabrfix=50; Kfix=10; for th in theta_values hU = channel(..., rfix, th, ...); ... строим w из пилотов при этом угле ... sT(ti,dBi) = 10*log10( |√P·hU'*w|^2 / σ^2 ); end
TDD-рекурсивность: в коде предполагается, что восходящий и нисходящий каналы совпадают по фазам/амплитудам (идеальная TDD-рекупроцитет и калибровка), поэтому фазы, оценённые по пилотам, годятся для downlink EGT.
Почему нормировка обязательна: так мы сравниваем режимы при неизменной суммарной передаваемой мощности . Без нормировки луч «рос» бы с некорректно.
Почему : это классическая «родительская» формула SNR после линейного формирования луча. Ваш «EGC-вариант» получается из неё как частный случай при и без нормировки канала на (что вы уже учли в своих последних версиях).
Роль : увеличивает точность фазовой оценки (дисперсия ошибки ). На графиках это даёт более высокие кривые SNR.
Соседние точки p1/p2: показывают уровень межнаправленных утечек (sidelobes/near-field mismatch) при использовании весов, заточенных под .
Идеологически это соответствует «цифровому EGT» в вашем ноутбуке: оценка фаз по пилотам, равноамплитудные веса и расчёт SNR по для пользователя и соседей. (lab3.py)
Если хотите, могу дать короткий чек-лист для полной «консистентности» (одна строка на правку): единый выбор vs , вариант формулы SNR (общая vs EGC-специализация), и где именно нельзя нормировать канал, чтобы (9) работала «как в теории».