# Help with Implementation of Transfer Function using Python or MATLAB

Relatively new to the DSP side and wanted some help to implement this approach using Python (with NumPy, SciPy libraries) or via MATLAB.

Background of the problem:

I'm running a linear dynamic loads model using ANSYS (commercial finite element software) and the run time is on the order of 4-6 days. The input is an acceleration time history and we're able to get time history data (acceleration, displacements, forces, moments, etc.) at specified locations.

Since the model is linear, I was able to create a transfer function using a SISO approach ($$x$$: input, $$y$$:output) using Python: $$H(f) = \frac{P_{yx}}{P_{xx}}$$

Where $$P_{yx}$$ is the cross-power spectral density and $$P_{xx}$$ is the power spectral density.

Question: As mentioned earlier, the time it takes to run a case is 4-6 days. We are in an iterative process where our input acceleration time histories are changing but our model remains the same. Therefore, I was curious what the "correct" approach is to leverage the existing transfer function, $$H(f)$$, and calculate a revised $$Y(f)$$ for a new $$X(f)$$.
The intent is that we can use the transfer function to generate new dynamic loads on a shorter duration than the full transient analysis.

If anyone has any guidance on this "process", any help is welcome!

Assuming you are modeling and measuring an LTI (Lime Time Invariant) System, we simply have

$$Y(\omega) = H(\omega)\cdot X(\omega)$$

$$y(t) = h(t)*x(t) = \int_{-\infty}^{+\infty} h(\tau)x(t-\tau)d\tau$$

The output signal can be obtained either by frequency domain multiplication or time domain convolution.

That's easy enough in theory but the devil is often in the details. The equations above are in continuous time but to implement this in a computer you need to sample in both time and frequency and use the DFT (Discrete Fourier Transform) instead.

Turns out if you want to accurately represent a continuous transfer function in the discrete domain, you need to avoid aliasing in BOTH domains, i.e. the sampling rate must be higher than twice the bandwidth and the DFT length must be longer than the length of the impulse response.

Nit-picky sidenote: this can actually NOT be done perfectly. A signal that's finite in time is infinite in frequency and vice versa. It can't be limited in both domains. However, we assume here that we can make the residual aliasing error small enough for our application.

The next point to consider with actual measurement is that your Signal to Noise ratio needs to be sufficiently large at ALL frequencies including DC and Nyquist. That can be problematic if your excitation signal $$X(\omega)$$ doesn't have a relatively smooth full range spectrum.

• Thank you for the help Hilmar; you are correct that we are modeling and measuring an LTI. From your experience, is there any recommended (or existing) Python or MATLAB functions/scripting that performs the frequency domain multiplication or time domain convolution? Sep 11, 2023 at 14:45
• @ganondorf29: there are quite a few and the best choice depends on your details. I would perhaps start with fftfilt() Sep 12, 2023 at 2:05