复现太赫兹IRS论文图2：散射场幅度平方曲线不符求排查

问题概述 #

我正在复现论文《Intelligent Reflecting Surfaces at Terahertz Bands: Channel Modeling and Analysis》中的图2（散射场幅度平方随观测角度变化曲线），当前用MATLAB生成的曲线趋势及y轴范围与论文原图存在明显差异，需要排查代码错误。

论文图2展示三条曲线：$L_x = L_y = \lambda/2$（Eq.9）、$L_x = L_y = \lambda/2$（Eq.10）、$L_x = L_y = 5\lambda$（Eq.9），反映散射场幅度平方随$\theta_r$的变化规律。

所用MATLAB代码 #

% Intelligent reflecting surfaces at terahertz bands: channel modelling and
% analysis
% figure 2

close all
clear all
clc

% setting latex inerpreter
set(groot,'defaulttextinterpreter','latex'); 
set(groot,'DefaultTextFontname', 'CMU Serif');
set(groot,'DefaultAxesFontName', 'CMU Serif');
set(groot,'DefaultTextFontSize',14)
set(groot,'DefaultAxesFontSize',14)

% phi: azimuth angle, theta: polar angle
% t: transmitter, r: receiver
% d: distance from origin
theta_r = 0:.1:90;   % receiver polar angel
theta_t = 30;        % transmitter polar angle
phi_r = 60;          % receiver azimuth angle
% amplitude of electric (E-field) of the incident plane
amplitudeE_i_squared = 1;            
f = 300e9;          % carrier frequency
D_r = 4;            % distances measured from the (0,0)th IRS element to Rx
c = 3e8;            % velocity of light
lambda = c/f;       % wavelength

%% 1. lx = ly = lambda/2, eq. 9
L_x = lambda/2;     % size of reflecting element in x
L_y = lambda/2;     % size of reflecting element in y

X = ((pi*L_x)/lambda) * sind(theta_r) * cosd(phi_r);
Y = ((pi*L_y)/lambda) * ( (sind(theta_r) * sind(phi_r)) - sind(theta_t));

F = cosd(theta_t)^2 * ( (cosd(theta_r).^2 * cosd(phi_r)^2) + ...
    sind(phi_r)^2);

E_s_NormSquared_eq9 = 10*log10( ((L_x*L_y)/lambda)^2 * ...
    (amplitudeE_i_squared/D_r^2) * F .* sinc(X).^2 .* sinc(Y).^2);

plot(theta_r,E_s_NormSquared_eq9,'b');
hold on

%% 2. lx = ly = lambda/2, eq. 10
E_s_NormSquared_eq10 = 10*log10( ((L_x*L_y)/lambda)^2 * ...
    (amplitudeE_i_squared/D_r^2) * F);

plot(theta_r,E_s_NormSquared_eq10,'k-.');

%% 3. lx = ly = 5*lambda, eq. 9
L_x = 5*lambda;
L_y = 5*lambda;

% lambda = 2*lambda;

X = ((pi*L_x)/lambda) * sind(theta_r) * cosd(phi_r);
Y = ((pi*L_y)/lambda) * ( (sind(theta_r) * sind(phi_r)) - sind(theta_t));

F = cosd(theta_t)^2 * ( (cosd(theta_r).^2 * cosd(phi_r)^2) + ...
    sind(phi_r)^2);

E_s_NormSquared_eq9 = 10*log10( ((L_x*L_y)/lambda)^2 * ...
    (amplitudeE_i_squared/D_r^2) * F .* ...
    sinc(X).^2 .* sinc(Y).^2);

% E_s_NormSquared_eq10 = 10*log10( ((L_x*L_y)/lambda)^2 * ...
% (amplitudeE_i_squared/D_r^2) * F);

plot(theta_r,E_s_NormSquared_eq9,'r');
% plot(theta_r,E_s_NormSquared_eq10,'r');

xlabel('$\theta_r [degrees]$');
ylabel('$||E_s||^2 [dB]$');
hl = legend('$L_x = L_y = \lambda/2$ (Eq. 9)',...
    '$L_x = L_y = \lambda/2$ (Eq. 10)',...
    '$L_x = L_y = 5\lambda$');
set(hl, 'Interpreter','latex')
xlim([min(theta_r) max(theta_r)]);

% saveas(gca, ['f2_', datestr(now,'ddmmyyyy'), '.png']);

核心错误点及修正方案 #

• MATLABsinc函数定义不匹配
  论文中的sinc函数标准定义为$\text{sinc}(x) = \frac{\sin(x)}{x}$，但MATLAB自带的sinc函数是$\text{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$，直接使用会导致计算结果偏差。
  修正：将所有sinc(X)替换为sin(X)./X，同时添加X(X==0)=eps;避免除零错误。

• 缺少归一化处理
  论文图中的曲线大概率做了归一化（将每条曲线的最大值设为0dB），而代码直接计算绝对dB值，导致y轴范围与趋势不符。
  修正：对每条曲线执行归一化，例如：

  E_s_NormSquared_eq9 = E_s_NormSquared_eq9 - max(E_s_NormSquared_eq9);
  

• 距离参数$D_r$的物理意义偏差
  代码中$D_r$定义为从(0,0)号IRS元素到接收机的距离，但论文中可能指从IRS阵列中心到接收机的距离，当阵列尺寸较大（如$5\lambda$）时，该差异会影响幅度计算。
  修正：确认论文中$D_r$的定义，若为阵列中心距离，需调整代码中距离项的计算逻辑。

内容的提问来源于stack exchange，提问作者Huá dé ní 華得尼