Wrong results performing THD calculation with FFT on square wave function in Python

Question

I buitlt a function which should perform a FFT first on a generated function, later on real data. Then I want the THD (Total Harmonic Distortion) from the signal. To proof the System I want to check it with functions like square wave and other waveforms. I have issues to get correct results after the FFT. I expect the amplitude on the output array I calculate are not correct or my THD calculation is wrong. Maybe its possible that the fail stars in the first steps too. In the future this code should handle real data from a measurement card. But first I wanna try it with a dummy code. I expected results around 48.3%. But my Results are about 43%.

I tried to round my timebase ("time" array) because the time array has not precise numbers throuth floating numbers. But that didnt help. I also printed the fft data in a .csv and the maximum points in array are correct so the resulting amplitude should be not.

samplingFrequency   = 100000
seconds             = 10
samplingInterval       = 1 / samplingFrequency
beginTime           = int(0)
endTime             = int(10); 
signal1Frequency     = 50

Multiplier           = 1000000
arranged_timefactor  = int((endTime/(1/signal1Frequency))*2)

time_unround       = np.linspace(beginTime, endTime, samplingFrequency*seconds, endpoint=False)     #Timebase

time  = np.around(time_unround, decimals=5)     #round timebase

forStart   = int(0)                              #define startpoint and endpoint with a number of samples
forEnd     = int(len(time))                      #from the signalperiod 

for x in range (forStart, forEnd):                                                  
    if signal.square(2 * np.pi * signal1Frequency * (x/Multiplier)) > 0.5:
        Startpoint = x
        break

print("Startpoint",Startpoint)

for y in range (forStart, forEnd):
    if signal.square(2 * np.pi * signal1Frequency * (y/Multiplier)) < -0.5:
        Endpoint = y
        break

print("Endpoint",Endpoint)
test = signal.square(2 * np.pi * signal1Frequency * time)
print("Endpointarr:", test[Endpoint])

time_arranged = time[Startpoint:Endpoint*arranged_timefactor]

amplitude1_arranged = signal.square(2 * np.pi * signal1Frequency * time_arranged)

fourierTransform = fftpack.fft(amplitude1_arranged)/len(amplitude1_arranged)           # Normalize amplitude

fourierTransform = fourierTransform[range(int(len(amplitude1_arranged)/2))] # Exclude sampling frequency

 

tpCount     = len(amplitude1_arranged)

values      = np.arange(int(tpCount/2))

timePeriod  = tpCount/samplingFrequency

frequencies = values/timePeriod

X_mag = np.array(np.abs(fourierTransform))
print (X_mag)

Vrms =  (X_mag[10*signal1Frequency])
Vrms1 = (X_mag[20*signal1Frequency])
Vrms2 = (X_mag[30*signal1Frequency])
Vrms3 = (X_mag[40*signal1Frequency])
Vrms4 = (X_mag[50*signal1Frequency])
Vrms5 = (X_mag[60*signal1Frequency])
Vrms6 = (X_mag[70*signal1Frequency])
Vrms7 = (X_mag[80*signal1Frequency])
Vrms8 = (X_mag[90*signal1Frequency])
Vrms9 = (X_mag[100*signal1Frequency])
Vrms10 = (X_mag[110*signal1Frequency])

THD = (sqrt(Vrms1**2+Vrms2**2+Vrms3**2+Vrms4**2+Vrms5**2+Vrms6**2+Vrms7**2+Vrms8**2+Vrms9**2+Vrms10**2))/Vrms
THDR = (THD/(sqrt(1+(THD**2))))
print ('THD ', THD*100,'%')
print ('THDR ', THDR*100,'%')