How can I use tfestimate with multi-experiment data?

I am trying to use tfestimate() with inputs and outputs garnered from multiple experiments. However, I'm not sure what the most appropriate way (i.e. best, mathematically rigorous) is to combine the experiments. I searched for an answer but couldn't find anything, so I've tried some approaches myself.

My primary question is this: Is there an established way to combine experimental spectral analysis data using the built-in MATLAB spectral analysis functions (tfestimate, periodogram, etc.), and/or can I tweak either of my approaches to both work well and be mathematically rigorous.

I've attached a script that compares the methods I've tried -- the important excerpts are shown below, but running the script will just produce three plots. First, I define a second order system and run some chirp signals on it.

N = 10001;
Time = linspace(0,10,N);
sys = tf([20 400], [1 20 400]);
x = zeros(N,3); y = zeros(N,3);

x(:,1) = chirp(Time,0,max(Time),100);
y(:,1) = lsim(sys,x(:,1),Time);

x(:,2) = chirp(Time,0,max(Time),10);
y(:,2) = lsim(sys,x(:,2),Time);

x(:,3) = chirp(Time,10,max(Time),40);
y(:,3) = lsim(sys,x(:,3),Time);

Then my first approach simply concatenates the inputs and outputs in sequence:

xvec = [x(:,1); x(:,2); x(:,3)]; 
yvec = [y(:,1); y(:,2); y(:,3)];

[Txy, F] = tfestimate(xvec,yvec,[],[],[],Fs);
[Cxy, Fc] = mscohere(xvec,yvec,[],[],[],Fs);

This returns great results; however, I think it's not entirely numerically rigorous, since MATLAB has no way of knowing that the signal is concatenated, so the transitions could lead to false readings. I assume that as long as the data is significantly larger than the number of concatenated signals, such effects are negligible, but I'm not well-versed enough in spectral analysis to really support that claim mathematically.

For my second approach, I get the transfer function sequence for each text separately, and then average the results.

[Txy1, F1] = tfestimate(x(:,1),y(:,1),[],[],[],Fs);
[Txy2, F2] = tfestimate(x(:,2),y(:,2),[],[],[],Fs);
[Txy3, F3] = tfestimate(x(:,3),y(:,3),[],[],[],Fs);
Txy = (Txy1 + Txy2 + Txy3)/3;
F = (F1 + F2 + F3)/3;

This approach provides acceptable results, but is nowhere near as nice as my first approach, which is conceivably performed using the same data (only the averaging occurs earlier in the analysis.