Interpreting Phase Bode Plot from Empirical Transfer Function

8 views (last 30 days)
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
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:
  • 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).

Accepted Answer

William Rose
William Rose on 3 Nov 2022
@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
William Rose on 5 Nov 2022
@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).

Sign in to comment.

More Answers (0)




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!