×
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Contact US

Log In

Come Join Us!

Are you an
Engineering professional?
Join Eng-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Frequency response in MatLab question

Frequency response in MatLab question

Frequency response in MatLab question

(OP)
I’m doing my very first steps in MatLab (5 days since I had a first look at it).
I’d like to ask for a bit of help:
A bit of pretext: I’m trying to learn about vehicle ride dynamics using Adams car/ride 4 post shaker rig. I created a model representing amg gt4 car and excite the car with chirp input (constant velocity sine sweep).
I intend to postprocess the results in MatLab. My aim is to be able to get frequency response functions (FRF) – in a way described in attached pdf file (p.8, 7 post rig data analysis)

So far I’m able to import data to MatLab (exporting as table/spreadsheet from A/car) and create plots of power spectrum, psd and transfer function.
In that attached pdf author is using mathlab TFE function to get FRF. I’m using TFESTIMATE (but I tried TFE as well). Problem is that his results (Y axis) look very different to what I get – I get results in Db and his results are in… who knows. It looks like ratio to me.


Is there a way that I can convert my Y axis from decibels to magnitude? There is a Matlab function db2mag but I can’t figure out how to apply it to my transfer function results.

Thank you,
Ted


RE: Frequency response in MatLab question

(OP)
Thank you, IRstuff.
Would you mind if I ask you to provide a little more detail regarding finding reference level?

Thank you,
Ted

RE: Frequency response in MatLab question

(OP)
Hmm...
Let me try again. I think that my poor English didn't allow me to express my question properly.
As far as I understand (please correct me if I'm wrong). Db is a ratio of magnitude of two signals expressed as 20*Log(A2/A1). Is there a way that I can get the results of tfestimate function to be expressed on linear scale?

Thank you,
Ted

RE: Frequency response in MatLab question

(OP)
Thank you, IRstuff,
Problem is that this is my 7th day with MatLab and I know very little at this point. It seems to me that conversion to Db is sort of imbedded in tfestimate MatLab function and I can't figure out how to undo the logarithm in the result plot of this function. When I write the function it may look like this: tfestimate(y,x,[],[],[],Fs) in brackets I can input window size, overlap, etc. But I don't see Log10 in input arguments options and so I can't choose not to enter it... I suspect that I'm missing something very simple but I can't find it in help.
May I ask you to suggest how exactly I can write it in MatLab to get rid of Logarithm?

Thank you,
Ted

RE: Frequency response in MatLab question

The dB reference value is 1, so if you were testing a spring mass system with a k of 100 it would start at 20 dB at 0 Hz.

It is a good idea to do the fft of each channel separately and then divide one by the other to check exactly what is going on scaling wise. Alternatively generate two known signals (say A and 10*A) and see what tfestimate gives.

Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376: Eng-Tips.com Forum Policies http://eng-tips.com/market.cfm?

RE: Frequency response in MatLab question

Quote:

The dB reference value is 1

It's "1" of something, but the OP should already know that, since the transfer function required dividing the response by the input, so all the OP needs to answer is WHAT is the input?

TTFN (ta ta for now)
I can do absolutely anything. I'm an expert! https://www.youtube.com/watch?v=BKorP55Aqvg
FAQ731-376: Eng-Tips.com Forum Policies forum1529: Translation Assistance for Engineers Entire Forum list http://www.eng-tips.com/forumlist.cfm

RE: Frequency response in MatLab question

(OP)
Thank you Greg,
There's definately a scaling mess which I have to fix first.
Regarding your 2nd reply - do you mean plotting abs of tfestimate? I didn't find any other way. May be it's my English playing tricks with me but I couldn't find any other way of plotting other than Db on y axis.

Thank you,
Ted

RE: Frequency response in MatLab question

Yes, the default is to display dB, but if you take the output of tfestimate you can display that vector however you like, it's just a vector of complex numbers. ie

tfout=tfestimate(usual stuff in here);
plot(real(tfout),imag(tfout))

will give you a nyquist plot




Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376: Eng-Tips.com Forum Policies http://eng-tips.com/market.cfm?

RE: Frequency response in MatLab question

(OP)
Thank you Greg,
That's what I'm doing now and with some fiddling with signal(exported from Adams) scaling things are starting to make sense.

Thank you,
Ted

RE: Frequency response in MatLab question

(OP)
In fact I was doing it differently:
tfout=tfestimate(usual stiff in here)
tfMag=abs(tfout)
tfPhase=angle(tfout)*57.3 to approximate phase in deg.
Then plot the above as two separate plots.
After rescaling Adams signals this gives me reasonable results with heave.
Doing it the way you suggested returns a very strange plot (sorry).
I will see how things work with other signals (pitch and cpl in particular) and report.

Thank you,
Ted

RE: Frequency response in MatLab question

BTW: how is your porpoising problem study going ?

So: here's what I do for computing lateral handling analysis:
I have data channels from road & sim data with the usual names. They are segments of length power of 2. Beginning and ending with a zero input for a second or 2.
nfft = pow2(round(log2(length(STEER))));
sps = 100; % sample rate
nsegs = size(STEER,2);
for nseg=1 %:nsegs
[sstxy,f] = tfestimate(STEER(:,nseg), LATACC(:,nseg),rectwin(nfft),nfft/2,nfft,sps); % ayg by steer transfer function
[rolltxy,f] = tfestimate(LATACC(:,nseg),ROLL (:,nseg),rectwin(nfft),nfft/2,nfft,sps); % roll by ay transfer function
[rollms2txy,f] = tfestimate(LATACC(:,nseg)*9.806,ROLL (:,nseg),rectwin(nfft),nfft/2,nfft,sps); % roll by ay transfer function
[betatxy,f] = tfestimate(LATACC(:,nseg),SIDSLP(:,nseg),rectwin(nfft),nfft/2,nfft,sps); % beta by ay transfer function
[betasteertxy,f]= tfestimate(STEER (:,nseg),SIDSLP(:,nseg),rectwin(nfft),nfft/2,nfft,sps); % beta by steer transfer function
[yawvtxy,f] = tfestimate(STEER(:,nseg), YAWVEL(:,nseg),rectwin(nfft),nfft/2,nfft,sps); % yawv by steer transfer function
[yawatxy,f] = tfestimate(STEER(:,nseg), YAWACC(:,nseg),rectwin(nfft),nfft/2,nfft,sps); % yaw accel by steer transfer function
[rollyawvtxy,f] = tfestimate(YAWVEL(:,nseg),ROLL (:,nseg),rectwin(nfft),nfft/2,nfft,sps); % roll by yawvel transfer function
[rollswatxy,f] = tfestimate(STEER(:,nseg) ,ROLL (:,nseg),rectwin(nfft),nfft/2,nfft,sps); % roll by steer transfer function

inx=find(f > 3.0 | f < .02);
f(inx)=[];
sstxy(inx)=[];
rolltxy(inx)=[];
rollms2txy(inx)=[];
betatxy(inx)=[];
betasteertxy(inx)=[];
yawvtxy(inx)=[];
yawatxy(inx)=[];
rollyawvtxy(inx)=[];
rollswatxy(inx)=[];

delta_f=f(2)-f(1);

handles.ss_gain = 100.*abs(sstxy); % convert to Bode form for gain
handles.ss_phase = angle(sstxy).*180/pi;
handles.yawv_gain = 100.*abs(yawvtxy);
handles.yawv_phase = angle(yawvtxy).*180/pi;
handles.yawa_gain = 100.*abs(yawatxy);
handles.yawa_phase = angle(yawatxy).*180/pi;
handles.roll_gain = abs(rolltxy);
handles.rollms2_gain = abs(rollms2txy);
handles.roll_phase = angle(rolltxy).*180/pi;
handles.rollms2_phase = angle(rollms2txy).*180/pi;
handles.beta_gain = abs(betatxy);
handles.beta_phase = angle(betatxy).*180/pi;
handles.betaswa_gain = 100*abs(betasteertxy);
handles.betaswa_phase = angle(betasteertxy).*180/pi;
handles.rollswa_gain = 100*abs(rollswatxy);
handles.rollswa_phase = angle(rollswatxy).*180/pi;
handles.rollyawv_gain = abs(rollyawvtxy);
handles.rollyawv_phase = angle(rollyawvtxy).*180/pi;
handles.f = f;
end

RE: Frequency response in MatLab question

(OP)
Thanks for great post Cibachrome,
I’m traveling right now but will reply as soon as I’m back home.

Thank you,
Ted

RE: Frequency response in MatLab question

%Try this for the heave response with your data: Seems to work pretty well
load Question
[heavetxy,freq] = tfestimate(act_dz,heave,rectwin(nfft),nfft/2 ,nfft,Fs); % heave by act_dz transfer function
inx=find(freq > 20 | freq < .01);
freq(inx)=[];
heavetxy(inx)=[];
heave_gain = abs(heavetxy);
heave_phase = angle(heavetxy).*180/pi;

figure('name','T_M_S Porpoising Study','numbertitle','off','menubar','none')
subplot(2,2,1)
plot(freq,heave_gain,'o')
ylabel('Heave Gain (units ?)')
grid on
title('This is all you need to do Ted:')
subplot(2,2,3)
plot(freq,heave_phase,'o')
grid on
xlabel('Frequency (Hz.)')
ylabel('Phase (deg)')
sidetext('BillCobb zzvyb6@yahoo.com')

% Now go get the transfer function:
options = optimset('MaxFunEvals',1000000,'TolFun',.0000001,'Display','on');
FM0=[ .1 1 0.15 0.23 ]
for n=1: 5
[FM,error] = lsqcurvefit('splane2_fit',FM0 ,freq ,heave_gain,[],[],options);
heavesys =tf([FM(1:2 )],[FM(3:4) 1 ])
FM0=FM+ .001*rand(1);
end

N_1=num2str(FM(1));
N_0=num2str(FM(2));
D_2=num2str(FM(3));
D_1=num2str(FM(4));
D_0= '1';

f=.0 :.02:20;
[mag,phase]=bode(heavesys,f);
subplot(2,2,1)
hold on
plot(f,squeeze(mag),'b-')
mytexstr = ['$\frac{H(s)}{Z(s)}$ = $\frac{' N_1 's+' N_0 '}{' D_2 's^2+' D_1 's+' D_0 '}$']
text('string',mytexstr,'interpreter','latex','fontsize',15,'units','norm','pos',[.25 .75],'Fontweight','Bold');

[wn,z]=damp(heavesys);
subplot(2,2,3)
hold on
plot(f,squeeze(phase),'b-')
mytexstr = [ '$\omega_n = ' num2str(wn(1),3) ' Hz $']
text('string',mytexstr,'interpreter','latex','fontsize',12,'units','norm','pos',[.4 .7],'Fontweight','Bold');
mytexstr = [ '$\zeta = ' num2str(z(1),3) '$']
text('string',mytexstr,'interpreter','latex','fontsize',12,'units','norm','pos',[.4 .5],'Fontweight','Bold');

subplot(2,2,2)
bode(heavesys)
ylabel('Heave Response')
grid on

subplot(2,2,4)
step(heavesys,5)
ylabel('Heave Response (units ??)')
grid on

You'll need this: Note: I'm only using magnitude for fit criteria. phase floats
function [mag] = splane2_fit(FM,f)
sys = tf([FM(1:2 )],[FM(3:4) 1]);
mag = squeeze(bode(sys,f));

RE: Frequency response in MatLab question

Anybody considered that the porpoising is a ride problem, not necessarily an aero problem ? Sure, aero can attenuate it, but the root cause seems to be something else.

BTW: Your 1000 Hz. sample rate is way overkill. What would help your activity is a longer data segment for a smaller delta_F frequency increment, multiple segments for ensemble averaging, adding a starting & ending signal of a second or 2, and a power of 2 number of samples. Some nonlinear effects smear the FFT data from the Question.mat file, but hopefully, you get the idea...

This is pitch by act_dz.

RE: Frequency response in MatLab question

In Newey's book he mentions it as being a problem with an early Leyton House car. Basically as the front tire flexed the front of the floor got too close the ground , blocking the flow, or stalling the airflow, so it lost downforce, the front came up, and then the whole cycle repeated itself. Hard to fix with a passive suspension because the main compliance is undamped.

Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376: Eng-Tips.com Forum Policies http://eng-tips.com/market.cfm?

RE: Frequency response in MatLab question

(OP)
OK, I'm back home and able to reply.
Many thanks for your help Cibachrome!
I already found answers to my initial questions - it was a matter of getting the units/scaling straight:

For a frequency response I was using (Live script):

Heave_ Chassis Acceleration Response Magnitude Ratio

[h_az_tf,n]=tfestimate(act_az,heave_az,winsize,overlap,nfft,Fs);
[h_az_tf2,n]=tfestimate(act_az2,heave_az2,winsize,overlap,nfft,Fs);
plot(n,abs(h_az_tf),'r')
xlim(l)
hold on
plot(n,abs(h_az_tf2),'b')
hold off
title("Heave Chassis Acceleration FRF Ratio")
xlabel('Frequency (Hz)')
ylabel('Magnitude')
legend('d1',"d2")

Heave_ Chassis Acceleration Response Phase

plot(n,angle(h_az_tf)*57.3,'r')
xlim(l)
ylim([-200 200])
hold on
plot(n,angle(h_az_tf2)*57.3,'b')
hold off
title("Heave Chassis Displacement FRF PHASE")
xlabel('Frequency (Hz)')
ylabel('Phase Angle')
legend('d1',"d2")

d1_HeaveAcc_FRF_Max=max(abs(h_az_tf))
d2_HeaveAcc_FRF_Max=max(abs(h_az_tf2))
d1_HeaveAcc_FRF_RMS=rms(abs(h_az_tf))
d2_HeaveAcc_FRF_RMS=rms(abs(h_az_tf2))


This was to compare two different damper settings.
I fugured out that I need to increse the simulation time to get better resolution at each frequency step - this seems to be inline with your suggestion from your latest post.

Now let me study your code a little. This is a bit above my current MatLab knowledge level (remember my MatLab experience is less than a month at this point).

Thank you,
Ted

PS: Attached is Acar simulation data files

RE: Frequency response in MatLab question

There is a Matlab Controls setup panel that is supposed to allow you to choose units and display schemes.

ctrlpref

But it doesn't seem to be working for me. Set the preferences and reload Matlab.

But, the danger is ruining all your previous codes ! (as it did mine !)

RE: Frequency response in MatLab question

(OP)
Thanks Cibachrome,
Thing is that requests in A/Car are predefined and it is a bit too complicated to redefine them atm. So I set the units/scaling manually as can be seen in DataImport.mlx that I uploaded.


Thank you,
Ted

RE: Frequency response in MatLab question

(OP)
So far I'm having problem with adopting your transfer function code to my current data (see above uploaded .mat file).
This portion:

% Now go get the transfer function:
options = optimset('MaxFunEvals',1000000,'TolFun',.0000001,'Display','on');
FM0=[ .1 1 0.15 0.23 ]
for n=1: 5
[FM,error] = lsqcurvefit('splane2_fit',FM0 ,freq ,heave_gain,[],[],options);
heavesys =tf([FM(1:2 )],[FM(3:4) 1 ])
FM0=FM+ .001*rand(1);
end

I added above code to this live script:

[h_az_tf,n]=tfestimate(act_az,heave_az,winsize,overlap,nfft,Fs);
[h_az_tf2,n]=tfestimate(act_az2,heave_az2,winsize,overlap,nfft,Fs);
plot(n,abs(h_az_tf),'r')
xlim(l)
hold on
plot(n,abs(h_az_tf2),'b')
hold off
title("Heave Chassis Acceleration FRF Ratio")
xlabel('Frequency (Hz)')
ylabel('Magnitude')
legend('d1',"d2")

plot(n,angle(h_az_tf)*57.3,'r')
xlim(l)
ylim([-200 200])
hold on
plot(n,angle(h_az_tf2)*57.3,'b')
hold off
title("Heave Chassis Displacement FRF PHASE")
xlabel('Frequency (Hz)')
ylabel('Phase Angle')
legend('d1',"d2")
heave_gain=abs(h_az_tf)

Here I add your code with changes apropriate for my code above

options = optimset('MaxFunEvals',1000000,'TolFun',.0000001,'Display','on');
FM0=[ .1 1 0.15 0.23 ]
for n=1: 5
[FM,error] = lsqcurvefit('splane2_fit',FM0 ,n ,abs(h_az_tf),[],[],options);
heavesys =tf([FM(1:2 )],[FM(3:4) 1 ])
FM0=FM+ .001*rand(1);
end

N_1=num2str(FM(1));
N_0=num2str(FM(2));
D_2=num2str(FM(3));
D_1=num2str(FM(4));
D_0= '1';

Problem that it gives me an error:
Error using lsqcurvefit (line 286)
Function value and YDATA sizes are not equal.

How do I get your function splane2_fit to be the same size as my YDATA (abs(h_az_tf)?

Thank you,
Ted

RE: Frequency response in MatLab question

"Thing is that requests in A/Car are predefined and it is a bit too complicated to redefine them atm"

I nearly fell of my chair at that one. In order to make our custom scripts work with /Car, MSC (or some lucky contractor) have written a black box that takes Car req files and re-mangles them back into the form we used to have in Chassis, where REQ numbers were fixed, for example REQ/1080 was always toe caster camber. I am not aware of the gory details.

Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376: Eng-Tips.com Forum Policies http://eng-tips.com/market.cfm?

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Eng-Tips forums is a member-only feature.

Click Here to join Eng-Tips and talk with other members! Already a Member? Login



News


Close Box

Join Eng-Tips® Today!

Join your peers on the Internet's largest technical engineering professional community.
It's easy to join and it's free.

Here's Why Members Love Eng-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close