matlab - 为什么随着散射体数量的变化，我会得到不一致的图形（不是按要求线性递减）？

Question

我正在尝试为相位光时域反射计（phi-OTDR）模拟反向散射信号的功率和衰减。我使用了研究论文中给出的公式：

总功率 p = p1+p2

对于任何 N 值（散射体数量），我应该以线性递减的方式获得长度与总功率的关系图，如下所示。这里 N = 100。

但是当我将 N 更改为任何其他值时，我不会得到相同的结果。相反，我得到的是这样的：这里，N=200：

我的代码如下所示：

NA = 0.13; %numerical aperture
ncore = 1.47 ; %fiber core refractive index
ndiff = 3e-3 ; %difference in index of refraction
lamda = 1.55 ;%Lambda in micrometer
c = 3e5; %speed of light in km/s
f = c/lamda*1e9; %pulse frequency
m = 4.55 ; %for single mode fiber
S = (NA/ncore)^2 / m ; %Capture coefficient 
asdB= (0.76+0.51*ndiff)/lamda^4; % Rayleigh attenuationin dB/km
as=asdB/4.34; % Rayleigh attenuation in km-1
adB= 0.22 ; % attenuationin dB/km
a=adB/4.34; % attenuation in km-1
ng=1.4681; % group refractive index
vg=c/ng; %group velocity in km/s
W = 1e-7; %Pulse width 100 ns
Z = 1; %distance of the fiber in km
N = 100; %number of scatterers
z = 0:Z/N:Z; %scatterers across the fiber length
fa=1.1e8; %Acousto Optic Modulator frequency in Hz
tau = (2*ng*z)/c;%delay time of scattering centers
t1=2*ng*Z/c; %roundtrip time
t=0:t1/100:t1; % at 100m resolution
L=vg*t/2; %length of fiber with 100m resolution- so steps of 0.01
Pin=0.8; % optical power of input lightwave

%% Sum of the optical power of the backscattered lightwave from each independent scatterer without distrbtion
P1=0;
for i=1:N+1
    P1 = P1 + Pin*S*as*exp((-2*a*c*tau(i))/ng)*(heaviside(t-tau(i)+W/2)-heaviside(t-tau(i)-W/2)); 
end
P2=0;
for i=1:N+1
    for j=i+1:N+1
        phi=2*pi*(f+fa)*(tau(i)-tau(j));
        P2 = P2 + Pin*cos(phi)*S*as*exp(-a*c*(tau(i)+tau(j))/ng)*(heaviside(t-tau(i)+W/2)-heaviside(t-tau(i)-W/2))...
            .*(heaviside(t-tau(j)+W/2)-heaviside(t-tau(j)-W/2));
    end
end
P2 = 2*P2;
Ptotal = P1+P2;
figure(1);
plot(L,Ptotal);
grid on;
xlabel('Distance in km ');
ylabel('Output Power in Watts '); 

标签： matlabgraphsignal-processingmatlab-figure