getting group delay by differentiating phase of reflectance

Hi all

I am trying to calculate and draw the group delay of contradirectional bragg grating.

I already calculated the reflectance and got its phase as shown in the figure.

when I differentiate the phase i got unlogical spikes like that figure

while the correct plot is shown in Yariv book as the following figure (trace denoted by taw)

to differentiate i tried two functions diff and gradient and both gave me the same spikes (which is wrong). These are the lines i used in matlab:

Reflection=S11;
R=abs(Reflection).^2;
R_phase2=unwrap(angle(Reflection));
phase_lambda=lambda*1e9;
dphasedlambda1 =-phase_lambda.^2./(2*pi*c).*( gradient(R_phase2) ./ gradient(phase_lambda));  
dphasedlambda2 =-phase_lambda(1:end-1).^2./(2*pi*c).*( diff(R_phase2) ./ diff(phase_lambda));  

Any help with that problem

[KALYAN ACHARJYA]

Could you please attach the complete code so that I can reproduce the results and identify any discrepancies?

[omnia]

sure... here is the code

sure. here is the code

c = 3e8;             % Speed of light in vacuum (m/s)
lambda_B=1560e-9;
L_Bragg=9.6000e-05;    
e=4e-6;
delta_n=0.01;
lambda = (1535:0.05 :1585)*1e-9;  
neff=2.463*ones(1,length(lambda));
V =1.006920722304421e+04;  % sinusoidal profile)
kappa_Bragg=1.006920722304421e+04

kappa_Bragg = 1.0069e+04

% Derived Parameters
beta = 2 * pi * neff ./ lambda;     
Delta_beta = 2 * pi * neff.*(1./lambda-1./lambda_B);

S=sqrt(V.^2-Delta_beta.^2/4);
cosh_SL=cosh(S*L_Bragg);
sinh_SL=sinh(S*L_Bragg);
Q = S .* cosh_SL + 1j*0.5*Delta_beta.* sinh_SL;

% waveguide constants
alpha_db=1.5e2;         
k=alpha_db*lambda/(40*pi*log(e));
gama=exp(-1j*beta*e+2*pi*k*e./lambda);

S11=-1j*kappa_Bragg.*sinh_SL./Q;

Reflection2=S11;
R=abs(Reflection2).^2;
R_phase2=unwrap(angle(Reflection2));
phase_lambda=lambda*1e9;
dphasedlambda2 =-phase_lambda(1:end-1).^2./(2*pi*c).*( diff(R_phase2) ./ diff(phase_lambda));  

% characteristics for one Bragg grating
figure
yyaxis left
plot(lambda*1e9,abs(S11).^2,'r')
ylabel('Reflection');
yyaxis right
plot(lambda*1e9,R_phase2,'b')
% ylabel('delay');
legend('Reflection','delay')
xlabel('wavelength (nm)');
set(gca,"FontSize",14)

●

[omnia]

Any help about that?

[David Goodmanson]

Hi omnia,

Here is the path of Reflection2 = S11 in the complex plane, per plot(Reflection2). As you can see, the curve goes right through the origin four times, each time leading to an instant change of the phase angle by pi. The unwrap function purportedly removes jumps in angle when abs(jump) >pi. In this case two of the jumps are just less than pi, but two of them are just barely greater than pi, and for some reason unwrap does not unwrap them. Even if it did, it would mask the real problem and not help.

I believe the problem here is that the curve passes right through the origin. If the curve passed nearby but not through, then you would get narrow peaks as in the Yariv plot. I think either you are not using his constants or your code is not quite right, maybe missing a loss term of some kind. Fixing that might also fix the magnitude plot, which does not have the flat top of the Yariv plot.

[omnia]

thank you for your answer.

the peak value can be flat if we increase the length so that the reflection will appraoch 1, hence getting the same curve as Yariv.

Although your answer did not solve the problem, it helped me to realize how to solve it.

the solution is simply piecewise differentiation. Finally it gives the same result.

It was very helpful, thank you

Answer 1

[charan]

Hi omnia,

The spikes in the plot of differential of phase is because there are points of non-differentiability in the function for phase. At these points, the derivative of the function shoots up to a large value and as a result the plot shows these values as spikes. The rest of the derivative values are very small when compared to them and hence the remaining portion of the plot appears to be a straight line. The second plot shows the required portion of the graph.

To resolve this issue you can use "ischange" function to extract the points of non-differentiability and calculate the differential at all the other points. The following code might help you:

%Original code
c = 3e8;             % Speed of light in vacuum (m/s)
lambda_B=1560e-9;
L_Bragg=9.6000e-05;    
e=4e-6;
delta_n=0.01;
lambda = (1535:0.05 :1585)*1e-9;  
neff=2.463*ones(1,length(lambda));
V =1.006920722304421e+04;  % sinusoidal profile)
kappa_Bragg=1.006920722304421e+04;
kappa_Bragg = 1.0069e+04;
% Derived Parameters
beta = 2 * pi * neff ./ lambda;     
Delta_beta = 2 * pi * neff.*(1./lambda-1./lambda_B);
S=sqrt(V.^2-Delta_beta.^2/4);
cosh_SL=cosh(S*L_Bragg);
sinh_SL=sinh(S*L_Bragg);
Q = S .* cosh_SL + 1j*0.5*Delta_beta.* sinh_SL;
% waveguide constants
alpha_db=1.5e2;         
k=alpha_db*lambda/(40*pi*log(e));
gama=exp(-1j*beta*e+2*pi*k*e./lambda);
S11=-1j*kappa_Bragg.*sinh_SL./Q;
Reflection2=S11;
R=abs(Reflection2).^2;
R_phase2=unwrap(angle(Reflection2));
phase_lambda=lambda*1e9;
%Code to plot the differential of phase
figure
plot(diff(R_phase2))

●

plot(diff(R_phase2(300:700)))

idx=find(ischange(R_phase2,'linear'))

idx = 1×6

90 288 401 604 718 924

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

phase_diff =diff(R_phase2(1:idx(1)-1));
plot(lambda(1:idx(1)-2)*1e9,phase_diff)
hold on
for i=1:length(idx)-1
    phase_diff =diff(R_phase2(idx(i)+2:idx(i+1)-1));
    plot(lambda(idx(i)+2:idx(i+1)-2)*1e9,phase_diff)
end
    phase_diff = diff(R_phase2(idx(end)+1:end));
        plot(lambda(idx(end)+2:end)*1e9,phase_diff)
hold off

Give the "method" argument of "ischange" function as "linear" to obtain the points of non-differentiability. You can refer to this documentation: https://www.mathworks.com/help/matlab/ref/ischange.html

[omnia]

that is completely true.

Thank you for the code too.