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:
- 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.
- 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?
- 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')
```