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')
```

Answer 1

First of all - I'm not familiar with Python and didn't check your code.

1a) For FFTs it is the usual that the frequency resolution equals the inverse length of your data record.

1b) As far as I understand your code (I may be wrong here), you take the FFT of the step response. The Bode diagram usually shows the FFT of the impulse response. Take care - the result may not be what you expect.

1c) In practice (due to transients, noise,...) the usual approach to measure the FRF is

• excite the system with a broadband signal, estimate the spectral density of the input and the output signal, divide the output spectral density by the input spectral density. This approach is well covered in the literature (and called "spectral method" or something similar

• exite the system with a sinuosoid, wait until the transient has died out, and read the magnitude and phase from your data. Repeat for every frequency where you want to know the FRF.

• There are other approaches to directly measure the FRF, like the Local Polynomial Method (or Local Rational Method), TRIMM, and more - they are covered in (scientific) publications, together with a lot of beautiful maths. The approaches I mentioned are called non-parametric methods (there are also parametric methods, but that may lead to far away for your use case)

2. The phase jump at Nyquist frequency is to be expected. The FFT spectrum of a rational signal will be symmetric to the DC-frequency - the magnite is axisymmetric, the phase point symmetric. For some FFTs the "negative" and "positive" frequencies are wrapped (this may depend on your FFT implementation). The result is a "jump" the phase.

3. The input signal should excite the entire frequency range of interest - check the spectral density function of the input signal (i.e. the step).