Total Harmonic Distortion Calculation

Hello!

I'm working on a script to calculate total harmonic distortion (THD) both for voltage and current.

I've created a model in SIMULINK that contains a section that generates a clean harmonic signal.

The SIMULINK THD block calculates an answer that I can prove by hand and with my scripts.

Now when I run another data set for a "real" AC output from a simulated h-bridge inverter, I can't get my scripts to match the SIMULINK THD block calculations. I can hand prove my outcomes, but I can't arrive at the SIMULINK outcomes.

My script uses the following method: I find the peak of the data set, assume that peak is the funamental. Grab it's indicie and then add 18 bins to search for the next harmonic. I used the ideal harmonic waveform to count the number of bins between harmonics. My fundamental frequency is 60hz. So I walk out a number of harmonics, capture their amplitudes and calculate the THD.

Fs= 120;  %Input Nyquist sampling frequency for signal
Ts=1/Fs; %sampling period
fund=60;  %Fundamental frequency
cycles=1; %Number of cycles evaluated
L=length(var); %grid imported from simulation run

%scaled fft of signal, looking for magnitude ONLY
Xspect=abs(Ts*fft(var))/3000;  %3000 to scale values
%take first half of signal (for positive values of the signal
X=Xspect(1:401);
%define frequency range for first half of signal;
f=Fs*(0:400);

%plot for visual
plot(f,X);
%tabulate for utility
f=f.'; %changed to column

tablex=table(f,X);
figure(1)
stem(f,X)
set(gca,'yscal','log')

freq=[0:round(L/2)-1]*fund/cycles;
freq=freq(1:401);
figure(2)
plot(freq,X)
set(gca,'yscal','log')

window=0;
n=18;       %bins between harmonics - will change with frequency
i=1;
num_harm=10;        % number of harmonics to be included in sum
while le(i,num_harm)     %set number of harmonics, default 6
    if i==1
        [peak(i),fundMaxIndex] =(max(X));   %capture fundamental
    else
         peak(i)=max(X(i*n+1:i*n+n,1))
         peak(i)=(peak(i)/peak(1))^2;
    end
    i=i+1;
end

sumPeak = sum(abs(peak(2:num_harm)));
THD = 100*sqrt(sumPeak)

When using this actual source, these results below are calculated by the build in THD block.

My outcomes compared to these values are the following: 1.55%(8.34%) , 0.5938%(4.37%) , 0.9696% (5.61%)

Any recommendations on my script or the approach i'm using? I'm trying to validate this calculation for my project.