Pulling a data point from a bode graph within a state-space simulation?

E_c=10^9; % N/m^3 rho_c=5000; % kg/m^3 Ro=.01; % m Ri=.009; % m h=1; % m A_c=(pi*Ro^2-pi*Ri^2); % m^2 V_c=(pi*Ro^2-pi*Ri^2)*h; % m^3 mc=V_c*rho_c; m2=1/3*mc; % kg m3=m2; % kg m4=m3; % kg m1=2; % kg

k1=10^4; % N/m k2=(E_c*A_c)/(1/6); % N/m k3=(E_c*A_c)/(1/3); % N/m k4=(E_c*A_c)/(1/3); % N/m k5=(E_c*A_c)/(1/6); % N/m k=[k1 k2 k3 k4 k5]; % N/m

m=[m1 m2 m3 m4]; % kg M=diag(m); % kg F=50; % N c=20; % Ns/m K=diag(k); % N/m

keq=(1/k(1)+1/k(2))^(-1); % N/m

A=[0,1,0,0,0,0,0,0; -keq/m(1),0,keq/m(1),0,0,0,0,0; 0,0,0,1,0,0,0,0; keq/m(2),0,-(k(3)+keq)/m(2),0,(k(3)/m(2)),0,0,0; 0,0,0,0,0,1,0,0; 0,0,(k(3)/m(3)),0,-(k(3)+k(4))/m(3),0,k(4)/m3,0; 0,0,0,0,0,0,0,1; 0,0,0,0,(k(4)/m(4)),0,-(k(4)+k(5))/m(4),-c/m(4)];

B=[0;1/m(1);0;0;0;0;0;0];

C=[10^3 0 0 0 0 0 0 0];

D=0;

sys=ss(A,B,C,D); [wn,zeta,P] = damp(sys); w=.05*min(wn); w2=.97*min(wn); bode(sys) -I want to pull the point when the frequency is equal to w and w2 in order to find the Magnitude ratios at those points.

Answer 1

[Star Strider]

Just ask for it!

From the documentation:

• [mag,phase,wout] = bode(sys) returns magnitudes, phase values, and frequency values wout corresponding to bode(sys).

If the data are not exactly the ones you want (or close enough to them), use the interp1 function to estimate them.

w=.05*min(wn);
w2=.97*min(wn);
[mag,phase,wout] = bode(sys);

mags = squeeze(mag);
phases = squeeze(phase);

mpi = interp1(wout, [mags,phases], [w; w2], 'linear', 'extrap');
Desired_Result = [[w; w2], mpi];
fprintf(1, '\n        w          Magnitude            Phase\n')
fprintf(1, '    %8.5f       %8.5f         %13.5E\n', Desired_Result')

          w          Magnitude            Phase
       3.27042        0.11537          -2.24180E-04
      63.44606        2.02070          -8.56391E-02

Answer 2

[John BG]

Jon

I have the impression that instead of sys=ss(A,B,C,D) perhaps you should use

    [b,a]=ss2tf(A,B,C,D)

where 'b' is the numerator and 'a' is the denominator of the transfer function of the system that you have initially defined with state-space matrices A,B,C,D

command ss works the other way round: from TF to state space.

So having introduced your lines, please next time tag the code by selecting it all and clicking on '{} code' that it formats, so that readers don't have to:

E_c=10^9; % N/m^3
rho_c=5000; % kg/m^3
Ro=.01; % m
Ri=.009; % m
h=1; % m
A_c=(pi*Ro^2-pi*Ri^2); % m^2
V_c=(pi*Ro^2-pi*Ri^2)*h; % m^3
mc=V_c*rho_c;
m2=1/3*mc; % kg
m3=m2; % kg
m4=m3; % kg
m1=2; % kg
k1=10^4; % N/m
k2=(E_c*A_c)/(1/6); % N/m
k3=(E_c*A_c)/(1/3); % N/m
k4=(E_c*A_c)/(1/3); % N/m
k5=(E_c*A_c)/(1/6); % N/m
k=[k1 k2 k3 k4 k5]; % N/m
m=[m1 m2 m3 m4]; % kg
M=diag(m); % kg
F=50; % N
c=20; % Ns/m
K=diag(k); % N/m
keq=(1/k(1)+1/k(2))^(-1); % N/m
A=[0,1,0,0,0,0,0,0; -keq/m(1),0,keq/m(1),0,0,0,0,0; 0,0,0,1,0,0,0,0; keq/m(2),0,-(k(3)+keq)/m(2),0,(k(3)/m(2)),0,0,0; 0,0,0,0,0,1,0,0; 0,0,(k(3)/m(3)),0,-(k(3)+k(4))/m(3),0,k(4)/m3,0; 0,0,0,0,0,0,0,1; 0,0,0,0,(k(4)/m(4)),0,-(k(4)+k(5))/m(4),-c/m(4)];
B=[0;1/m(1);0;0;0;0;0;0];
C=[10^3 0 0 0 0 0 0 0];
D=0;

state-space to transfer function numerator and denominator

[b,a]=ss2tf(A,B,C,D)

b=    1.0e+21 *[      0              0      0.000000000000000   0.000000000000000   0.000000000005449  0.000000000552632   0.000015020048113   0.000361067798879   6.624086603893607]

a= 1.0e+22 *[   0.000000000000000   0.000000000000000   0.000000000000001   0.000000000000111   0.000000003009263 0.000000072741618   0.001339001277477   0.000316834641557   5.673582679221318 ]

Now launch the filter design tool

fdatool

and introduce 'b' and 'a'.

The transfer function is:

Tell fdatool to generate a function that you can use as an independent code block:

function Hd = ss2tf_filter1
%SS2TF_FILTER1 Returns a discrete-time filter object.

% MATLAB Code
% Generated by MATLAB(R) 8.5.1 and the Signal Processing Toolbox 7.0.
% Generated on: 08-Mar-2016 02:16:53

% Equiripple FIR Lowpass filter designed using the FIRPM function.

% All frequency values are in Hz.
Fs = 48000;  % Sampling Frequency

Fpass = 9600;            % Passband Frequency
Fstop = 12000;           % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 16;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);

% [EOF]

you can also tell fdatool to generate a SIMULINK block. I have attached this SIMULINK block for your convenience, click on the attached ss2tf_simulink_slx.txt and change file name to ss2tf_simulink.slx

This last odd step is because file extensions like fdatool exports are not currently accepted in MATLAB CENTRAL file uploads. I already asked them to enhance the range of files that can be attached, moreover fdatool being part of MATLAB, but no answer yet.

Well, hope it helps.

If you find this answer of any help solving your question, please click on the thumbs-up vote link,

thanks in advance

John

[e-mail]