# Bode plot from from step response data

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,
})
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,
})
# plot results
fig = plt.figure(1)
fig.clear()
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.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.grid(which='both', axis='both')
ax2.legend(loc='upper center')
$$$$
`
• For 1, if it's IIR, then the longer the better, if it's FIR then... the step response should be clearly stable over a finite number of samples. For 2, it's possible that you end up with a zero due to IFFT, otherwise it's entirely up to the transfer function whether there is a zero at Nyquist or not. For 3, I would imagine plenty of padding would do the trick, though it would not magically extend the bandwidth. Aug 20 '20 at 13:33

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)

1. 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.

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