Interpreting Phase Bode Plot from Empirical Transfer Function

I have an experiment where a vehicle body is supported and a freely hanging wheel with an added known imbalance mass on the tire. The tire is driven up to 40 mph and allowed to spin down from there or held at a constant speed depending on the test. This system is treated as a mass-spring damper where the input to the system is the vertical component of the known imbalance mass's centrifugal force calculated as

Taking the fourier-transform of the measured raw wheel acceleration in the vertical direction, , and the calculated vertical force, , then creating an empirical transfer function, , yields a magnitude plot which consistently matches a "theoretical" bode plot made using nominal vehicle suspension values and the tf() function from the following transfer function.

The phase plot of this empirical transfer function is obtained using the following code :

Fs = 5000;

L = length(t);
xfft = Fs*(0:(L/2))/L;

% Calculate FFT for measured acceleration
ff_a_z = fft(a_z);
a_z_ff = ff_a_z/L;
a_z_ff_p1 = a_z_ff(1:L/2+1);
a_z_ff_p1(2:end-1) = 2*a_z_ff_p1(2:end-1);

% Calculate FFT for calculated force
ff_F_t = fft(F_t);
F_t_ff = ff_F_t/L;
F_t_ff_p1 = F_t_ff(1:L/2+1);
F_t_ff_p1(2:end-1) = 2*F_t_ff_p1(2:end-1);

% Empirical Transfer Function
G = a_z_ff_p1./F_t_ff_p1';

%Get phase of TF in degrees
figure
G_ang = unwrap(angle(G))*180/pi;
semilogx(xfft*2*pi, G_ang)
hold on
%semilogx(wout, phase)
legend(["Empirical TF", "Theoretical TF"])

The empirical phase of consistently takes this shape with 2 peaks but with the peaks varying up to 600 degrees. This does not take the 2nd overdamped phase form that I was expecting to see. I am aware that the empirical transfer function will not be accurate at frequencies that the system has not been excited in, but in the range of interest, 10 to 52 rad/s, there is significant excitation. I noticed that the individual phases of the vertical acceleration and force are every large, but they grow larger together and somewhat "cancel" eachother out.

That being said, because I am dividing a complex number by a complex number in , then the following is also true:

Trying this method shows a much more "consistent" phase of but only after "rewrapping" the phase so it lies within a reasonable range

% Phase of individual F(s) and A_z(s) spectra
F_t_ang = unwrap(angle(F_t_ff_p1))*180/pi;
a_z_ang = unwrap(angle(a_z_ff_p1))*180/pi;

% Pure Difference
G_ang2 = a_z_ang - F_t_ang';

% With "Rewrapping"
G_diff = mod(((a_z_ang - F_t_ang') + 180), 360) - 180;

And a zoomed in view:

Questions:

• Why does my phase consistently show a ~-150 degree phase angle after "rewrapping". I'm expecting a positive angle between 180 and 0?
• Is it appropriate for me to "rewrap" it as I have it here?
• Will any non-linearities have a large consistent affect on the phase plot of this system?
• Is this indicative that the system order is more complex then a mass-spring-damper?
• Can I use knowledge of my system order (i.e. 2 zeros @ 0 Hz and 2 poles) to try and manipulate this phase diagram to make more sense?

The code I'm using was "validated" against a simulated model of the mass-spring-damper using an integrator and the empirical and theoretical magnitude and phase bode plots matched almost perfectly (as expected from a perfectly linear system).

Answer 1

[William Rose]

@Perry, great quesiton.

Could you post a file with the experimental data a_z and F_t, which you used in the claculaitons above?

As you know, unwrap adds or subtracts integr multiples of 360. The phase angle is approximately -150 in phase plot #3. The is no way that unwrap could get there from a phase of 0 to +180. Therefore I think the experimental phase is not what you expect, or the way you are calculating it is not giving a statitically reliable result.

There is a considerable literature on esitating tranfer functions. Your way of esitmating it is the obvious way:

G(f)=Y(f)/X(f), where f=frequency and X and Y are the DFTs of y(t) and x(t). (Discrete Fourier transform (DFT) is the name for the Fourier transform of a sampled signal. The fast Fourier transform is an algorithm for computing the DFT efficiently.) Your way makes sense, but it turns out that it is quote susceptible to noise in the signal. A better esitmate, which is the best esitmate in a least squares error sense, is obtained by choppng x and y into N segments (possibly overlapping), and computing X*(f)Y(f) and X*(f)X(f) for each segment, where X*(f) is the complex conjugate of X(f). Then you average togther the N estimates of X*Y(f) to get an estimate Sxy(f), and you average together the N esitmates of X*X to get an esitmate Sxx(f). THen you compute the transfer funciton estimate . There is a nice analogy between this approach and the estimation of slope, b, in time domain data, where the model is y(t)=a+bx(t).

Matlab has functions to do the above, or you can do it yourself. (I have.) I wonder if that approach, which is the recommended approach, ould give a phase estimate closer to the theoretical value.

[William Rose]

@Perry,

[edit: correct typos in text; code unchanged]

To estimate the transfer function by the method I described above, use the built-in function tfestimate().

[G,f]=tfestimate(x,y,window,noverlap,nfft,fs)

If the length of x and y is a multiple of 16, I would try

N=length(x);
window=N/16; nfft=window; noverlap=window/2;
[G1,f1]=tfestimate(x,y,window,noverlap,nfft,fs);

This will chop x,y into 16+15=31 half-overlapped segments. I would also try

window=N/8; nfft=window; noverlap=window/2;
[G2,f2]=tfestimate(x,y,window,noverlap,nfft,fs);

This will chop x,y into 8+7=15 half-overlapped segments. I would also try

window=N/4; nfft=window; noverlap=window/2;
[G3,f3]=tfestimate(x,y,window,noverlap,nfft,fs);

This will chop x,y into 4+3=7 half-overlapped segments.

The advantage of more segments in the estimate is that the transfer function estimate will have a smaller error bar.* The disadvantage of more segments is that the frequency resolution is lower, i.e. in the spectrum is bigger.

*The estimation of error bars on the spectrum is complicated, and deserves a separate Q and A. The attached Technical Report for the Air Force Flight Dynamics Laboratory has the error bar formulas which have been used by pretty much everyone, including Matlab, since then.

[Perry]

Thanks for your reply William,

The data you asked for is attached. My understanding of unwrap is that it adds or subtracts 2pi, but the angle() function places the angle within a range of [-pi pi] which is reasonable with my "rewrapped" answer.

The tfestimate() function seems particularly useful for what I'm trying to do becuase it allows me to specify resolution and range of frequencies for tf estimation. I got an error from your code, but I slightly adapted it to make a bit more sense with my setup. Interestingly, this gives me a similar answer as what I tried previously.

%% Transfer Function Estimate
start_freq = 51.6764;
stop_freq = 19.1710;

Fs = 5000;
N = length(a_z)
window = N/8; 
nfft = [stop_freq:2.5:start_freq]./(2*pi); 
%nfft = []; 
noverlap = window/2;
[G1, f1] = tfestimate(F_t,a_z, window, noverlap, nfft, Fs)

figure(7)
semilogx(f1*2*pi, 20*log10(abs(G1)))
xlabel("Frequency [rad/s]")
ylabel("Magnitude [dB]")
title("Magnitude")

figure(8)
semilogx(f1*2*pi, angle(G1)*180/pi)
xlabel("Frequency [rad/s]")
ylabel("Phase [degrees]")
title("Phase")

The various window sizes provide pretty similar answers but N/8 made the most sense.

I'm still a bit confused how I'm not getting the positive phase I was expecting from transfer function.

[William Rose]

@Perry, those plots of magnitude and phase look good!

A possible cause of an unexpected phase is that the transducers and hardware used to collect the force and acceleration data have time delay, or phase, responses, which are not equal to each other. This alters the phase of the estimated transfer function. I have had to do careful measurement and correction to remove this effect from transfer funciton estimates derived from experimental data, because the pressure transducers had a different phase and gain response than the flow transducers.

[Perry]

Hey William,

You deserve to know that I made a significant blunder when calculating F(t), I was using a sin when I should have been using a cos. This caused a 90 degree shift in the phase. Here is the corrected graph which lays in the expected range of 0 to 180.

Regardless, your suggestion of using tfestimate() was very helpful. I know you mentioned that the error of this can be evaluated is well. Is there a MATLAB command for doing this? Thanks again!

[William Rose]

@Perry, Good job figuring out the phase!

I thought tfestimate() had an option to return error bars on the transfer funciton estimate, but it does not. pwelch() returns the optional argument pxxc, which is a 2-column vector containing the upper and lower confidence interval for the power spectrum. Interestingly, the confidence interval pxxc, for the estimate pxx, returned by pwelch(), does not follow the formula I expected it to follow - the formula which is in Bendat & Piersol (1971) p.329, and in Otnes & Enochson (1972), p.217. The pwelch() help cites textbooks by Hayes (1996) and by Stoica & Moses (2005). I don't have those books.

I have done some original research to compare the transfer function error distribution from tfestimate() to the formula of N.R. Goodman (1965). The error distribution from tfestimate() does not agree exactly with the prediction of Goodman, but they are not too far off. Interestingly, Welch did not publish his method until 1967, so the method for estimtating the numerator and denominator of the gain estimate, which was assumed by Goodman, was not the Welch method we use now. Furthermore, the Matlab help for tfestimate cites Vold, Crowley, & Rocklin, "New ways of estimating frequency response functions", (1984). I don't have that paper. How did Vold et al. estimate the transfer function? Does their paper give a formula for a confidence interval for the TF estimate?

I really should get my hands on the "modern" (1996, 2005) textbooks and Vold et al (1984).