Bode plot from from step response data

Question

For quality control, I want to measure the frequency response of hardware I built. What I can easily generate on the hardware is a step response. Thus, I want to derive the magnitude and phase from the step response. This works well in general, however I am looking for ways to improve and understand basic relationships.

Remark: Quality control is probably done best by comparing the measured step response directly with a "golden" step response. However, I want to provide also numbers and visuals of magnitude and phase as the original system specification is given in these terms.

Find below the python code of what I have, comparing my approach with a direct Bode calculation from the transfer function. Questions:

1. Is there a formula or rule of thumb of how the duration/length in time of the recorded step response and the sample rate used affect the quality/uncertainty of magnitude and phase? One obvious thing is, that Bode data for low frequencies with longer period than the recorded duration will have lower quality.
2. I always get a phase jump at the nyquist frequency, regardless of the sample rate and signal duration used. I am struggling to identify the root cause of this. Does it have an explanation from DSP theory or is it my implementation?
3. What other things should I pay attention to to get high quality Bode data?

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import signal # to generate example data

def stepResp2Bode(stepResponse, fS, minFreqRes=None):
    if minFreqRes is not None:
        if fS/stepResponse.size > minFreqRes:
            print('Minimum frequency resolution not achieved!')
            return # not enough samples for required frequency resolution
    impResp = np.diff(stepResponse, axis=0, prepend=0) # get impulse response by differentiating step response
    fftResult = np.fft.rfft(impResp, axis=0) # do fft
    bodeF = np.fft.rfftfreq(impResp.size, d=1/fS) # create frequency vector for fft-Results
    bodeAmplDB = 20 * np.log10(np.abs(fftResult)) # extract amplitude in dB
    bodePhase = np.angle(fftResult) # extract phase in radians
    results = pd.DataFrame({'frequency':bodeF,
                        'amplitude in dB':bodeAmplDB,
                        'phase in radians':bodePhase,
                        })
    return results

if __name__ == "__main__":
    # generate example data - in the application this would be real world measured data
    fS = 50 # sampling frequency
    wc = (2*np.pi*1) # corner frequency of example low pass
    t = np.arange(0, 5, 1/fS) # vector of time values
    sys = signal.lti([1.0], [1/wc, 1.0]) # generate lowpass filter
    stepResponse = signal.step(sys, T=t)[1] # calculate step response of low pass filter at given time values
    # Calculation
    bodeStepResults = stepResp2Bode(stepResponse, fS, minFreqRes=1) # from step response to magnitude and phase
    w, mag, phase = signal.bode(sys) # bode plot directly from system description
    bodeScipyResults = pd.DataFrame({'frequency':(w/(2*np.pi)),
                    'amplitude in dB':mag,
                    'phase in radians':(phase/180*np.pi),
                    })
    # plot results
    fig = plt.figure(1)
    fig.clear()
    ax1 = fig.add_subplot(211)
    ax1.semilogx(bodeStepResults['frequency'], bodeStepResults['amplitude in dB'], '-', label='from step response')
    ax1.semilogx(bodeScipyResults['frequency'], bodeScipyResults['amplitude in dB'], '-', label='from signal.bode')
    ax1.set_xlabel('Frequency  in Hz')
    ax1.set_ylabel('Amplitude  in dB')
    ax1.grid(which='both', axis='both')
    ax1.legend(loc='upper center')
    ax2 = fig.add_subplot(212)
    ax2.semilogx(bodeStepResults['frequency'], bodeStepResults['phase in radians'], '-', label='from step response')
    ax2.semilogx(bodeScipyResults['frequency'], bodeScipyResults['phase in radians'], '-', label='from signal.bode')
    ax2.set_xlabel('Frequency  in Hz')
    ax2.set_ylabel('Phase in radians')
    ax2.grid(which='both', axis='both')
    ax2.legend(loc='upper center')
```