RIS辅助太赫兹通信信道特征建模与MATLAB仿真分析

目录

1.太赫兹信道基础损耗原理

2.RIS级联信道原理

3. 系统性能指标原理

4. MATLAB实现程序

可重构智能表面（RIS）由大量无源反射单元构成，通过独立调控各单元相移，对入射太赫兹波进行相位补偿与波束整形，将发射端信号经反射链路定向汇聚至接收端。太赫兹通信核心瓶颈为空间传播损耗与分子选择性吸收，RIS通过构建虚拟视距链路，绕过遮挡、补偿高损耗，显著提升信道增益与链路可靠性。

1.太赫兹信道基础损耗原理

RIS是由大量可编程超材料单元组成的平面阵列，每个单元可独立调控入射电磁波的相位和幅度。在太赫兹频段(0.1-10THz)，信号传播面临严重的路径损耗、分子吸收和阻挡效应。RIS通过智能反射建立虚拟视距(LoS)链路，补偿太赫兹信道的严重衰落。系统模型包含三段链路：发射端(Tx)→RIS链路、RIS反射、RIS→接收端(Rx)链路。RIS包含𝑁个反射单元，每个单元施加相移𝜃𝑛 ，形成对角相移矩阵。

太赫兹信号传播损耗包含空间扩散损耗与分子吸收损耗，水蒸气等分子在特定频段形成共振吸收峰，导致信道频率选择性衰减。

2.RIS级联信道原理

RIS信道为发射端-RIS-接收端双跳级联链路，信道响应为两段信道与相移矩阵的级联。

3. 系统性能指标原理

4. MATLAB实现程序

clear; clc; close all; %% ========== 参数设置 ========== c = 3e8; % 光速 f_center = 0.3e12; % 中心频率300GHz lambda = c / f_center; % 波长 d_s = lambda / 2; % RIS单元间距 B = 10e9; % 带宽10GHz Pt_dBm = 20; % 发射功率20dBm Pt = 10^((Pt_dBm-30)/10); sigma2_dBm = -80; % 噪声功率 sigma2 = 10^((sigma2_dBm-30)/10); K_factor = 10; % 莱斯因子 Nt = 4; % 发射天线数 num_monte = 500; % 蒙特卡洛次数 %% ========== 1. 太赫兹分子吸收系数(简化模型) ========== f_range = linspace(0.1e12, 1e12, 2000); % 水蒸气吸收简化模型(多个吸收峰) f_res = [0.557e12, 0.752e12, 0.988e12]; % 共振频率 S = [0.3, 0.15, 0.4]; % 吸收强度 alpha_half = [15e9, 20e9, 25e9]; % 半宽度 kappa_abs = zeros(size(f_range)); for i = 1:length(f_res) kappa_abs = kappa_abs + S(i) * (alpha_half(i)./((f_range - f_res(i)).^2 + alpha_half(i)^2)); end kappa_abs = kappa_abs * 1e-7; % 标度调整 figure('Name','分子吸收系数','Position',[100 100 800 500]); plot(f_range/1e12, kappa_abs, 'b-', 'LineWidth', 2); xlabel('频率 (THz)'); ylabel('吸收系数 \kappa_{abs} (m^{-1})'); title('太赫兹频段分子吸收系数'); grid on; set(gca,'FontSize',12); %% ========== 2. 路径损耗 vs 距离(不同频率) ========== d_range = 1:0.5:100; freqs_test = [0.1e12, 0.3e12, 0.5e12, 0.8e12, 1e12]; colors = lines(length(freqs_test)); figure('Name','路径损耗vs距离','Position',[100 100 900 600]); for fi = 1:length(freqs_test) f_tmp = freqs_test(fi); % 计算该频率的吸收系数 k_tmp = 0; for i = 1:length(f_res) k_tmp = k_tmp + S(i)*(alpha_half(i)/((f_tmp-f_res(i))^2+alpha_half(i)^2)); end k_tmp = k_tmp * 1e-7; PL_dB = zeros(size(d_range)); for di = 1:length(d_range) d = d_range(di); L_spread = (4*pi*f_tmp*d/c)^2; L_abs = exp(k_tmp * d); PL_dB(di) = 10*log10(L_spread * L_abs); end plot(d_range, PL_dB, '-', 'LineWidth', 2, 'Color', colors(fi,:)); hold on; end xlabel('距离 (m)'); ylabel('路径损耗 (dB)'); title('太赫兹路径损耗随距离变化(含分子吸收)'); legend('0.1 THz','0.3 THz','0.5 THz','0.8 THz','1.0 THz','Location','southeast'); grid on; set(gca,'FontSize',12); %% ========== 3. 路径损耗 vs 频率(不同距离) ========== distances_test = [5, 10, 20, 50]; figure('Name','路径损耗vs频率','Position',[100 100 900 600]); for di = 1:length(distances_test) d = distances_test(di); PL_freq = zeros(size(f_range)); for fi = 1:length(f_range) f_tmp = f_range(fi); L_spread = (4*pi*f_tmp*d/c)^2; k_tmp = 0; for i = 1:length(f_res) k_tmp = k_tmp + S(i)*(alpha_half(i)/((f_tmp-f_res(i))^2+alpha_half(i)^2)); end k_tmp = k_tmp * 1e-7; L_abs = exp(k_tmp * d); PL_freq(fi) = 10*log10(L_spread * L_abs); end plot(f_range/1e12, PL_freq, '-', 'LineWidth', 2); hold on; end xlabel('频率 (THz)'); ylabel('路径损耗 (dB)'); title('太赫兹路径损耗随频率变化'); legend('d=5m','d=10m','d=20m','d=50m','Location','southeast'); grid on; set(gca,'FontSize',12); %% ========== 4. RIS辅助信道建模函数 ========== % 定义辅助函数 generate_array_response = @(N, angle, d_spacing, wavelength) ... exp(1j * 2*pi * d_spacing / wavelength * (0:N-1).' * sin(angle)) / sqrt(N); %% ========== 5. SNR vs RIS单元数 ========== N_list = [4, 8, 16, 32, 64, 100, 144, 196, 256]; d_TR = 30; d_RR = 15; % Tx-RIS距离, RIS-Rx距离 % 计算路径损耗 f0 = f_center; k0 = 0; for i = 1:length(f_res) k0 = k0 + S(i)*(alpha_half(i)/((f0-f_res(i))^2+alpha_half(i)^2)); end k0 = k0 * 1e-7; PL_TR = (4*pi*f0*d_TR/c)^2 * exp(k0*d_TR); PL_RR = (4*pi*f0*d_RR/c)^2 * exp(k0*d_RR); SNR_random = zeros(length(N_list), num_monte); SNR_optimal = zeros(length(N_list), num_monte); SNR_noRIS = zeros(1, num_monte); % 直接链路(被遮挡,额外30dB损耗) d_direct = 40; PL_direct = (4*pi*f0*d_direct/c)^2 * exp(k0*d_direct) * 1e3; for mi = 1:num_monte % 无RIS直接链路 h_direct = sqrt(1/(2*PL_direct)) * (randn(1,Nt) + 1j*randn(1,Nt)); w_mrt = h_direct' / norm(h_direct); SNR_noRIS(mi) = Pt * abs(h_direct * w_mrt)^2 / sigma2; for ni = 1:length(N_list) N = N_list(ni); % Tx->RIS信道 (N x Nt) phi_t = pi/6; phi_r = pi/4; psi_r = pi/3; a_RIS = generate_array_response(N, phi_r, d_s, lambda); a_Tx = generate_array_response(Nt, phi_t, d_s, lambda); H_LoS = a_RIS * a_Tx'; H_NLoS = (randn(N,Nt) + 1j*randn(N,Nt)) / sqrt(2); H_TR = sqrt(N/PL_TR) * (sqrt(K_factor/(K_factor+1))*H_LoS + sqrt(1/(K_factor+1))*H_NLoS); % RIS->Rx信道 (N x 1) phi_rx = pi/5; a_Rx = generate_array_response(N, phi_rx, d_s, lambda); h_NLoS_rr = (randn(N,1) + 1j*randn(N,1)) / sqrt(2); h_RR = sqrt(N/PL_RR) * (sqrt(K_factor/(K_factor+1))*a_Rx + sqrt(1/(K_factor+1))*h_NLoS_rr); % 波束成形向量 w = ones(Nt,1)/sqrt(Nt); % 随机相移 theta_rand = 2*pi*rand(N,1); Phi_rand = diag(exp(1j*theta_rand)); h_cas_rand = h_RR.' * Phi_rand * H_TR * w; SNR_random(ni,mi) = Pt * abs(h_cas_rand)^2 / sigma2; % 最优相移 v = h_RR .* (H_TR * w); theta_opt = -angle(v); Phi_opt = diag(exp(1j*theta_opt)); h_cas_opt = h_RR.' * Phi_opt * H_TR * w; SNR_optimal(ni,mi) = Pt * abs(h_cas_opt)^2 / sigma2; end end SNR_rand_avg = 10*log10(mean(SNR_random, 2)); SNR_opt_avg = 10*log10(mean(SNR_optimal, 2)); SNR_noRIS_avg = 10*log10(mean(SNR_noRIS)); figure('Name','SNR vs RIS单元数','Position',[100 100 900 600]); plot(N_list, SNR_opt_avg, 'ro-', 'LineWidth', 2, 'MarkerSize', 8); hold on; plot(N_list, SNR_rand_avg, 'bs--', 'LineWidth', 2, 'MarkerSize', 8); plot(N_list, SNR_noRIS_avg * ones(size(N_list)), 'k-.', 'LineWidth', 2); xlabel('RIS单元数 N'); ylabel('平均SNR (dB)'); title('SNR随RIS单元数变化'); legend('RIS最优相移','RIS随机相移','无RIS(遮挡)','Location','northwest'); grid on; set(gca,'FontSize',12); %% ========== 6. 可达速率 vs 发射功率 ========== Pt_dBm_range = -10:2:30; N_fixed = 64; Rate_opt = zeros(size(Pt_dBm_range)); Rate_rand = zeros(size(Pt_dBm_range)); Rate_noRIS = zeros(size(Pt_dBm_range)); for pi_idx = 1:length(Pt_dBm_range) Pt_tmp = 10^((Pt_dBm_range(pi_idx)-30)/10); snr_opt_tmp = zeros(1, num_monte); snr_rand_tmp = zeros(1, num_monte); snr_no_tmp = zeros(1, num_monte); for mi = 1:num_monte N = N_fixed; phi_t = pi/6; phi_r = pi/4; a_RIS = generate_array_response(N, phi_r, d_s, lambda); a_Tx = generate_array_response(Nt, phi_t, d_s, lambda); H_LoS = a_RIS * a_Tx'; H_NLoS = (randn(N,Nt) + 1j*randn(N,Nt)) / sqrt(2); H_TR = sqrt(N/PL_TR) * (sqrt(K_factor/(K_factor+1))*H_LoS + sqrt(1/(K_factor+1))*H_NLoS); phi_rx = pi/5; a_Rx = generate_array_response(N, phi_rx, d_s, lambda); h_NLoS_rr = (randn(N,1) + 1j*randn(N,1)) / sqrt(2); h_RR = sqrt(N/PL_RR) * (sqrt(K_factor/(K_factor+1))*a_Rx + sqrt(1/(K_factor+1))*h_NLoS_rr); w = ones(Nt,1)/sqrt(Nt); % 最优 v = h_RR .* (H_TR * w); theta_opt = -angle(v); Phi_opt = diag(exp(1j*theta_opt)); h_opt = h_RR.' * Phi_opt * H_TR * w; snr_opt_tmp(mi) = Pt_tmp * abs(h_opt)^2 / sigma2; % 随机 Phi_rand = diag(exp(1j*2*pi*rand(N,1))); h_rand = h_RR.' * Phi_rand * H_TR * w; snr_rand_tmp(mi) = Pt_tmp * abs(h_rand)^2 / sigma2; % 无RIS h_direct = sqrt(1/(2*PL_direct)) * (randn(1,Nt) + 1j*randn(1,Nt)); snr_no_tmp(mi) = Pt_tmp * abs(h_direct * w)^2 / sigma2; end Rate_opt(pi_idx) = B * mean(log2(1 + snr_opt_tmp)) / 1e9; Rate_rand(pi_idx) = B * mean(log2(1 + snr_rand_tmp)) / 1e9; Rate_noRIS(pi_idx) = B * mean(log2(1 + snr_no_tmp)) / 1e9; end figure('Name','可达速率vs发射功率','Position',[100 100 900 600]); plot(Pt_dBm_range, Rate_opt, 'ro-', 'LineWidth', 2, 'MarkerSize', 8); hold on; plot(Pt_dBm_range, Rate_rand, 'bs--', 'LineWidth', 2, 'MarkerSize', 8); plot(Pt_dBm_range, Rate_noRIS, 'g^-.', 'LineWidth', 2, 'MarkerSize', 8); xlabel('发射功率 (dBm)'); ylabel('可达速率 (Gbps)'); title('可达速率随发射功率变化 (N=64, f=0.3THz)'); legend('RIS最优相移','RIS随机相移','无RIS','Location','northwest'); grid on; set(gca,'FontSize',12); %% ========== 7. 信道增益CDF ========== N_cdf = 64; gain_opt = zeros(1, 2000); gain_rand = zeros(1, 2000); gain_noRIS = zeros(1, 2000); for mi = 1:2000 N = N_cdf; a_RIS = generate_array_response(N, pi/4, d_s, lambda); a_Tx = generate_array_response(Nt, pi/6, d_s, lambda); H_TR = sqrt(N/PL_TR)*(sqrt(K_factor/(K_factor+1))*a_RIS*a_Tx' + ... sqrt(1/(K_factor+1))*(randn(N,Nt)+1j*randn(N,Nt))/sqrt(2)); a_Rx = generate_array_response(N, pi/5, d_s, lambda); h_RR = sqrt(N/PL_RR)*(sqrt(K_factor/(K_factor+1))*a_Rx + ... sqrt(1/(K_factor+1))*(randn(N,1)+1j*randn(N,1))/sqrt(2)); w = ones(Nt,1)/sqrt(Nt); v = h_RR .* (H_TR * w); Phi_opt = diag(exp(-1j*angle(v))); gain_opt(mi) = abs(h_RR.' * Phi_opt * H_TR * w)^2; Phi_rand = diag(exp(1j*2*pi*rand(N,1))); gain_rand(mi) = abs(h_RR.' * Phi_rand * H_TR * w)^2; h_d = sqrt(1/(2*PL_direct))*(randn(1,Nt)+1j*randn(1,Nt)); gain_noRIS(mi) = abs(h_d * w)^2; end figure('Name','信道增益CDF','Position',[100 100 900 600]); [f1,x1] = ecdf(10*log10(gain_opt)); [f2,x2] = ecdf(10*log10(gain_rand)); [f3,x3] = ecdf(10*log10(gain_noRIS)); plot(x1, f1, 'r-', 'LineWidth', 2); hold on; plot(x2, f2, 'b--', 'LineWidth', 2); plot(x3, f3, 'g-.', 'LineWidth', 2); xlabel('信道增益 (dB)'); ylabel('CDF'); title('信道增益累积分布函数 (N=64)'); legend('RIS最优相移','RIS随机相移','无RIS','Location','southeast'); grid on; set(gca,'FontSize',12); %% ========== 8. RIS波束图案 ========== N_beam = 64; Ny = 8; Nz = 8; % 8x8 RIS target_az = pi/6; target_el = pi/4; % 计算最优相移用于指向目标方向 phase_config = zeros(Ny, Nz); for ny = 1:Ny for nz = 1:Nz phase_config(ny,nz) = -2*pi*d_s/lambda * ((ny-1)*sin(target_az)*cos(target_el) + (nz-1)*sin(target_el)); end end az_scan = linspace(-pi/2, pi/2, 360); el_scan = linspace(-pi/2, pi/2, 360); beam_az = zeros(size(az_scan)); beam_2d = zeros(length(el_scan), length(az_scan)); for ai = 1:length(az_scan) AF = 0; for ny = 1:Ny for nz = 1:Nz phase_steer = 2*pi*d_s/lambda * ((ny-1)*sin(az_scan(ai))*cos(target_el) + (nz-1)*sin(target_el)); AF = AF + exp(1j*(phase_config(ny,nz) + phase_steer)); end end beam_az(ai) = abs(AF)^2 / N_beam^2; end for ai = 1:length(az_scan) for ei = 1:length(el_scan) AF = 0; for ny = 1:Ny for nz = 1:Nz phase_steer = 2*pi*d_s/lambda*((ny-1)*sin(az_scan(ai))*cos(el_scan(ei)) + (nz-1)*sin(el_scan(ei))); AF = AF + exp(1j*(phase_config(ny,nz) + phase_steer)); end end beam_2d(ei,ai) = abs(AF)^2 / N_beam^2; end end figure('Name','RIS方位面波束图','Position',[100 100 900 600]); plot(az_scan*180/pi, 10*log10(beam_az+1e-10), 'b-', 'LineWidth', 2); xlabel('方位角 (度)'); ylabel('归一化增益 (dB)'); title('RIS反射波束方位面方向图 (8\times8, 目标30°)'); ylim([-40, 0]); grid on; set(gca,'FontSize',12); xline(target_az*180/pi, 'r--', 'LineWidth', 1.5); legend('波束图','目标方向'); figure('Name','RIS二维波束图','Position',[100 100 800 700]); imagesc(az_scan*180/pi, el_scan*180/pi, 10*log10(beam_2d+1e-10)); colorbar; caxis([-30 0]); xlabel('方位角 (度)'); ylabel('俯仰角 (度)'); title('RIS二维波束图 (8\times8阵列)'); set(gca,'FontSize',12); colormap('jet'); %% ========== 9. 可达速率 vs RIS-Rx距离 ========== d_RR_range = 5:2:80; N_dist = 64; Rate_dist_opt = zeros(size(d_RR_range)); Rate_dist_rand = zeros(size(d_RR_range)); for di = 1:length(d_RR_range) d_rr = d_RR_range(di); PL_rr_tmp = (4*pi*f0*d_rr/c)^2 * exp(k0*d_rr); snr_o = zeros(1, num_monte); snr_r = zeros(1, num_monte); for mi = 1:num_monte N = N_dist; a_RIS = generate_array_response(N, pi/4, d_s, lambda); a_Tx = generate_array_response(Nt, pi/6, d_s, lambda); H_TR = sqrt(N/PL_TR)*(sqrt(K_factor/(K_factor+1))*a_RIS*a_Tx' + ... sqrt(1/(K_factor+1))*(randn(N,Nt)+1j*randn(N,Nt))/sqrt(2)); a_Rx = generate_array_response(N, pi/5, d_s, lambda); h_RR = sqrt(N/PL_rr_tmp)*(sqrt(K_factor/(K_factor+1))*a_Rx + ... sqrt(1/(K_factor+1))*(randn(N,1)+1j*randn(N,1))/sqrt(2)); w = ones(Nt,1)/sqrt(Nt); v = h_RR .* (H_TR * w); Phi_opt = diag(exp(-1j*angle(v))); snr_o(mi) = Pt * abs(h_RR.' * Phi_opt * H_TR * w)^2 / sigma2; Phi_rand = diag(exp(1j*2*pi*rand(N,1))); snr_r(mi) = Pt * abs(h_RR.' * Phi_rand * H_TR * w)^2 / sigma2; end Rate_dist_opt(di) = B * mean(log2(1 + snr_o)) / 1e9; Rate_dist_rand(di) = B * mean(log2(1 + snr_r)) / 1e9; end figure('Name','速率vs RIS-Rx距离','Position',[100 100 900 600]); plot(d_RR_range, Rate_dist_opt, 'ro-', 'LineWidth', 2, 'MarkerSize', 8); hold on; plot(d_RR_range, Rate_dist_rand, 'bs--', 'LineWidth', 2, 'MarkerSize', 8); xlabel('RIS-Rx 距离 (m)'); ylabel('可达速率 (Gbps)'); title('可达速率随RIS-Rx距离变化 (N=64, f=0.3THz)'); legend('最优相移','随机相移','Location','northeast'); grid on; set(gca,'FontSize',12); %% ========== 10. 相位量化比特数影响 ========== N_quant = 64; bits_list = [1, 2, 3, 4, 6, 8]; Pt_scan = -10:2:30; Rate_quant = zeros(length(bits_list)+1, length(Pt_scan)); for pi_idx = 1:length(Pt_scan) Pt_tmp = 10^((Pt_scan(pi_idx)-30)/10); for mi = 1:num_monte N = N_quant; a_RIS = generate_array_response(N, pi/4, d_s, lambda); a_Tx = generate_array_response(Nt, pi/6, d_s, lambda); H_TR = sqrt(N/PL_TR)*(sqrt(K_factor/(K_factor+1))*a_RIS*a_Tx' + ... sqrt(1/(K_factor+1))*(randn(N,Nt)+1j*randn(N,Nt))/sqrt(2)); a_Rx = generate_array_response(N, pi/5, d_s, lambda); h_RR = sqrt(N/PL_RR)*(sqrt(K_factor/(K_factor+1))*a_Rx + ... sqrt(1/(K_factor+1))*(randn(N,1)+1j*randn(N,1))/sqrt(2)); w = ones(Nt,1)/sqrt(Nt); v = h_RR .* (H_TR * w); theta_cont = -angle(v); % 连续相移(理想) Phi_c = diag(exp(1j*theta_cont)); h_c = h_RR.' * Phi_c * H_TR * w; Rate_quant(end, pi_idx) = Rate_quant(end, pi_idx) + ... B * log2(1 + Pt_tmp*abs(h_c)^2/sigma2) / 1e9 / num_monte; % 量化相移 for bi = 1:length(bits_list) nb = bits_list(bi); levels = 2^nb; theta_q = round(theta_cont / (2*pi/levels)) * (2*pi/levels); Phi_q = diag(exp(1j*theta_q)); h_q = h_RR.' * Phi_q * H_TR * w; Rate_quant(bi, pi_idx) = Rate_quant(bi, pi_idx) + ... B * log2(1 + Pt_tmp*abs(h_q)^2/sigma2) / 1e9 / num_monte; end end end figure('Name','相位量化影响','Position',[100 100 900 600]); markers = {'v-','s-','d-','^-','p-','h-','o-'}; for bi = 1:length(bits_list) plot(Pt_scan, Rate_quant(bi,:), markers{bi}, 'LineWidth', 1.8, 'MarkerSize', 6); hold on; end plot(Pt_scan, Rate_quant(end,:), 'k-', 'LineWidth', 2.5); xlabel('发射功率 (dBm)'); ylabel('可达速率 (Gbps)'); title('不同相位量化比特数下的可达速率 (N=64)'); legend_str = cell(1, length(bits_list)+1); for bi = 1:length(bits_list) legend_str{bi} = sprintf('%d-bit', bits_list(bi)); end legend_str{end} = '连续(理想)'; legend(legend_str, 'Location', 'northwest'); grid on; set(gca,'FontSize',12); %% ========== 11. 莱斯因子影响 ========== K_list = [0, 1, 5, 10, 20, 50]; N_K = 64; SNR_K = zeros(length(K_list), num_monte); for ki = 1:length(K_list) K_tmp = K_list(ki); for mi = 1:num_monte N = N_K; a_RIS = generate_array_response(N, pi/4, d_s, lambda); a_Tx = generate_array_response(Nt, pi/6, d_s, lambda); H_TR = sqrt(N/PL_TR)*(sqrt(K_tmp/(K_tmp+1))*a_RIS*a_Tx' + ... sqrt(1/(K_tmp+1))*(randn(N,Nt)+1j*randn(N,Nt))/sqrt(2)); a_Rx = generate_array_response(N, pi/5, d_s, lambda); h_RR = sqrt(N/PL_RR)*(sqrt(K_tmp/(K_tmp+1))*a_Rx + ... sqrt(1/(K_tmp+1))*(randn(N,1)+1j*randn(N,1))/sqrt(2)); w = ones(Nt,1)/sqrt(Nt); v = h_RR .* (H_TR * w); Phi_opt = diag(exp(-1j*angle(v))); SNR_K(ki,mi) = Pt * abs(h_RR.' * Phi_opt * H_TR * w)^2 / sigma2; end end figure('Name','莱斯因子影响','Position',[100 100 900 600]); SNR_K_avg = 10*log10(mean(SNR_K, 2)); bar(categorical(string(K_list)), SNR_K_avg, 0.6, 'FaceColor', [0.3 0.6 0.9]); xlabel('莱斯因子 K'); ylabel('平均SNR (dB)'); title('不同莱斯因子下的平均SNR (N=64, 最优相移)'); grid on; set(gca,'FontSize',12); %% ========== 12. 信道频率响应(宽带) ========== N_wb = 64; Nf = 512; f_sub = linspace(f_center - B/2, f_center + B/2, Nf); H_freq_opt = zeros(1, Nf); H_freq_rand = zeros(1, Nf); % 固定一次信道实现 a_RIS_c = generate_array_response(N_wb, pi/4, d_s, lambda); a_Tx_c = generate_array_response(Nt, pi/6, d_s, lambda); H_TR_base = sqrt(K_factor/(K_factor+1))*a_RIS_c*a_Tx_c' + ... sqrt(1/(K_factor+1))*(randn(N_wb,Nt)+1j*randn(N_wb,Nt))/sqrt(2); a_Rx_c = generate_array_response(N_wb, pi/5, d_s, lambda); h_RR_base = sqrt(K_factor/(K_factor+1))*a_Rx_c + ... sqrt(1/(K_factor+1))*(randn(N_wb,1)+1j*randn(N_wb,1))/sqrt(2); w = ones(Nt,1)/sqrt(Nt); % 在中心频率设计最优相移 v_c = h_RR_base .* (H_TR_base * w); theta_opt_c = -angle(v_c); Phi_opt_c = diag(exp(1j*theta_opt_c)); Phi_rand_c = diag(exp(1j*2*pi*rand(N_wb,1))); for fi = 1:Nf f_tmp = f_sub(fi); lambda_tmp = c / f_tmp; k_tmp = 0; for i = 1:length(f_res) k_tmp = k_tmp + S(i)*(alpha_half(i)/((f_tmp-f_res(i))^2+alpha_half(i)^2)); end k_tmp = k_tmp * 1e-7; PL_tr = (4*pi*f_tmp*d_TR/c)^2 * exp(k_tmp*d_TR); PL_rr = (4*pi*f_tmp*d_RR/c)^2 * exp(k_tmp*d_RR); % 频率相关的阵列响应 a_r = generate_array_response(N_wb, pi/4, d_s, lambda_tmp); a_t = generate_array_response(Nt, pi/6, d_s, lambda_tmp); H_TR_f = sqrt(N_wb/PL_tr) * (sqrt(K_factor/(K_factor+1))*a_r*a_t' + ... sqrt(1/(K_factor+1))*H_TR_base); a_rx = generate_array_response(N_wb, pi/5, d_s, lambda_tmp); h_RR_f = sqrt(N_wb/PL_rr) * (sqrt(K_factor/(K_factor+1))*a_rx + ... sqrt(1/(K_factor+1))*h_RR_base); H_freq_opt(fi) = abs(h_RR_f.' * Phi_opt_c * H_TR_f * w)^2; H_freq_rand(fi) = abs(h_RR_f.' * Phi_rand_c * H_TR_f * w)^2; end figure('Name','宽带频率响应','Position',[100 100 900 600]); plot((f_sub-f_center)/1e9, 10*log10(H_freq_opt/max(H_freq_opt)), 'r-', 'LineWidth', 1.5); hold on; plot((f_sub-f_center)/1e9, 10*log10(H_freq_rand/max(H_freq_rand)), 'b-', 'LineWidth', 1.5); xlabel('相对中心频率偏移 (GHz)'); ylabel('归一化信道增益 (dB)'); title('RIS辅助太赫兹宽带信道频率响应 (N=64)'); legend('最优相移(中心频率设计)','随机相移'); grid on; set(gca,'FontSize',12); ylim([-30 2]); %% ========== 13. 信道冲激响应(功率延迟分布) ========== h_time_opt = ifft(sqrt(H_freq_opt)); h_time_rand = ifft(sqrt(H_freq_rand)); tau = (0:Nf-1) / B * 1e9; % 时延(ns) figure('Name','功率延迟分布','Position',[100 100 900 600]); stem(tau(1:50), 10*log10(abs(h_time_opt(1:50)).^2/max(abs(h_time_opt).^2)+1e-10), 'r', 'LineWidth', 1.5, 'MarkerSize', 4); hold on; stem(tau(1:50), 10*log10(abs(h_time_rand(1:50)).^2/max(abs(h_time_rand).^2)+1e-10), 'b', 'LineWidth', 1.5, 'MarkerSize', 4); xlabel('时延 (ns)'); ylabel('归一化功率 (dB)'); title('RIS辅助太赫兹信道功率延迟分布'); legend('最优相移','随机相移'); grid on; set(gca,'FontSize',12); ylim([-40 2]); %% ========== 14. 多用户场景 - 不同用户角度的波束增益 ========== N_mu = 100; % 10x10 RIS user_angles = [-40, -20, 0, 20, 40] * pi/180; az_fine = linspace(-pi/2, pi/2, 500); figure('Name','多用户波束','Position',[100 100 900 600]); colors_mu = lines(length(user_angles)); for ui = 1:length(user_angles) % 相移对准第ui个用户 phase_mu = -2*pi*d_s/lambda * (0:N_mu-1).' * sin(user_angles(ui)); beam_mu = zeros(size(az_fine)); for ai = 1:length(az_fine) sv = exp(1j*2*pi*d_s/lambda*(0:N_mu-1).'*sin(az_fine(ai))); beam_mu(ai) = abs(sum(exp(1j*phase_mu) .* sv))^2 / N_mu^2; end plot(az_fine*180/pi, 10*log10(beam_mu+1e-10), '-', 'LineWidth', 2, 'Color', colors_mu(ui,:)); hold on; end xlabel('角度 (度)'); ylabel('归一化增益 (dB)'); title('RIS波束指向不同用户方向 (N=100)'); legend('-40°','-20°','0°','20°','40°','Location','south'); grid on; ylim([-40 0]); set(gca,'FontSize',12); fprintf('\n===== 仿真完成 =====\n'); fprintf('共生成 14 幅图\n'); fprintf('1. 分子吸收系数\n'); fprintf('2. 路径损耗vs距离\n'); fprintf('3. 路径损耗vs频率\n'); fprintf('4. SNR vs RIS单元数\n'); fprintf('5. 可达速率vs发射功率\n'); fprintf('6. 信道增益CDF\n'); fprintf('7. RIS方位面波束图\n'); fprintf('8. RIS二维波束图\n'); fprintf('9. 速率vs RIS-Rx距离\n'); fprintf('10. 相位量化比特数影响\n'); fprintf('11. 莱斯因子影响\n'); fprintf('12. 宽带频率响应\n'); fprintf('13. 功率延迟分布\n'); fprintf('14. 多用户波束\n');

测试结果如下：