Own Welch procedure vs pwelch function

I have a 36001 long time series (time step 0.1s) representing the wind speed during 1 hour at a position:

Dummybodies.midwind(:,1)

This wind speed is represented by a Kaimal spectrum.

I have coded my own Welch procedure and then found out there is a Maltab built-in function. Thus, I thought it'd be interesting to compare them. However, it does not fit and I'd like to know what I've done wrong.

Here is a summary of the Welch procedure

In my case N=36001, and I've created 23 segments of L=2048 samples. The first 161 points of the time serie are discarded.

%% Own analysis
Fs=10; % sample freq
L=2048; % number of points in each segment
segmentn=23; % number of segments
overlap=0.25; % overlap between segments
thrownpts=161; % points discarded from the initial time serie

S1test=Dummybodies.midwind(:,1); % signal/time serie to analyse

segments=zeros(L,segmentn);
P1_1=zeros(segmentn,1025);

S1test=S1test(thrownpts:end);
% Segments
for k=1:segmentn
    segments(:,k)=S1test(1+(k-1)*(1-overlap)*L:k*L-(k-1)*L*overlap); 
    s=segments(:,k)';
    w=hamming(L)'; % Hamming window
    w=sqrt(L)*w/sqrt(sum(w.^2)); % Normalization
    s=w.*s;
    Y=fft(s); 

    f = Fs*(0:(n/2))/n;
    P2 = abs(Y/L);
    P1 = P2(:,1:n/2+1);
    P1(:,2:end-1) = 2*P1(:,2:end-1);
    P1_1(k,:) = P1; % line k of P1_1 is the single sided amplitude of segment k
end

PSDwperso=mean(P1_1); % average on all the segments

Here is how I've done the analysis using Matlab built-in function pwelch

%% Matlab analysis
S1test=Dummybodies.midwind(:,1);
Fs=10;
L=2048;
overlap=0.25;
thrownpts=161;

S1test=S1test(thrownpts+1:end);

[PSDwMatlab,f2] = pwelch(S1test,L,overlap*L,L,Fs);

Answer 1

[Mathieu NOE]

hello

for me , your function compute the autospectrum and not a PSD

to do a PSD :

PSD = Y.^2/df;

with

df = Fs/nfft;