Phase shift of transfer function using MATLAB

Question

I'm trying to compute the phase shift of a transfer function (accelerance) in MATLAB, but I'm having issues.

The transfer function was computed by dividing the magnitude of the system output by the input (please see here).

I've been using pwelch() to compute all of the functions in the frequency domain, but from the best of my understanding, the Welch PSD estimate outputs a real-valued vector, so there's no way to get a phase shift.

When I try to use the fft() function, the result (magnitude) doesn't look right, so I'm not sure if I would trust the phase shift. Below is the code I used to compute the transfer function:

% Welch Acceleration
[pa,f] = pwelch(ay,1000,10,0:10:10000,Fs);
% Welch Force
[pf,f] = pwelch(fy,1000,10,0:10:10000,Fs);
% Transfer function
H = pa./pf;

If someone could point me in the right direction it would be appreciated.

Answer 1

You are correct that using only pwelch would only yield real results, thus no phase information. If your signals do not contain any noise and when taking the fft they don't show Gibbs-like phenomenon (due to misalignment of the start and end values of the signals) one should be able to obtain both good phase and magnitude information (provided that the input signal excites all frequencies).

These two assumptions regarding the signals are kind of solved by using pwelch. Namely, by cutting the signal into multiple segments the noise should hopefully average to zero when adding their fft results and by multiply each segment (still in the time domain) by a window function the start and end points are better aligned (though might introduce some spectral leakage). However, pwelch also multiplies each component with it's complex conjugate in order to obtain the spectral density, which is indeed real valued.

Instead of the spectral density from pwelch, one could also calculate the cross spectral density. This is also used by the tfestimate function from matlab. By dividing the cross spectral density of the input and output signals by the spectral density of the input signal, both using the Welch method/segments and windows, then (somewhat) reliable phase can also be obtained even in the presence of noise. Though, you might want to play around with the number of segments. Namely, the more segments the more likely it hopefully is that the noise averages to zero, but more shorter segments also means less spectral resolution.