Filter coefficients from frequency response with invfreqs/invfreqz

Hi,

I have a linear time invariant system (spring-mass-damper), the input of the system is a trapezoidal acceleration power spectral density applied to the base of the system, the output is the root mean square of the acceleration of the mass. To evaluate the rms of the output I would like to use lyap, but I can use it only if the input of the system is a white noise.

I want to design a filter that has as input a white noise ( a flat power spectral density) and as output the input of my LTI system: the trapezoidal acceleration power spectral density. I need only the coefficients of the filter so that I can evaluate the space-state form of the filter and of the filter + LTI (needed to use lyap).

I don't understand exactly how to use invfreqs/invfreqz. The magnitude response of the designed filter doesn't match with what I would like to have.

Here's the code:

clear all; close all; clc

%% Defining the trapezoidal APSD

df = 0.1;

f_a = 20; f_b = 100; f_c = 300; f_d = 2000; % freq points

s_l = 3; s_r = -5; % slopes

APSD_m = 1; % flat mid PSD

f_l = f_a:df:f_b; f_m = f_b:df:f_c; f_r = f_c:df:f_d;

APSD_l = APSD_m * (f_l/f_b) .^ (s_l/(10*log10(2))); % left PSD

APSD_r = APSD_m * (f_r/f_c) .^ (s_r/(10*log10(2))); % right PSD

APSD_m = APSD_m * ones(1,length(f_m)); % mid PSD

f_APSD = [f_l,f_m,f_r];

APSD = [APSD_l,APSD_m,APSD_r];

%% Filter invfreqz

desiredResponse = sqrt(APSD); % The magnitude response I would like to have

normalizedFrequency = f_APSD./f_APSD(end) * pi;

[b,a] = invfreqz(desiredResponse,normalizedFrequency,10,13,[],30); % to get filter coefficients

zplane(b,a)

[h,w] = freqz(b,a,length(f_APSD)); % to see the frequency response of the filter

figure()

loglog(f_APSD,APSD,'-k',f_APSD,abs(h).^2,'-g')

legend('Target APSD','Filter Response')

[Paul]

Hi Fabio,

What is "lyap"? If a Matlab function, can you provide a link to the doc page?

[Star Strider]

lyap, dlyap

Answer 1

[Paul]

Assuming that f_APSD is in Hz, the filter, h, below gives a reasonable match to APSD.

df = 0.1;
f_a = 20; f_b = 100; f_c = 300; f_d = 2000; % freq points
s_l = 3; s_r = -5; % slopes
APSD_m = 1; % flat mid PSD
f_l = f_a:df:f_b; f_m = f_b:df:f_c; f_r = f_c:df:f_d;
APSD_l = APSD_m * (f_l/f_b) .^ (s_l/(10*log10(2))); % left PSD
APSD_r = APSD_m * (f_r/f_c) .^ (s_r/(10*log10(2))); % right PSD
APSD_m = APSD_m * ones(1,length(f_m)); % mid PSD
f_APSD = [f_l,f_m,f_r];
APSD = [APSD_l,APSD_m,APSD_r];
%% Filter invfreqz
%{
desiredResponse =  sqrt(APSD); % The magnitude response I would like to have
normalizedFrequency = f_APSD./f_APSD(end) * pi; 
[b,a] = invfreqz(desiredResponse,normalizedFrequency,10,13,[],30); % to get filter coefficients
zplane(b,a)
[h,w] = freqz(b,a,length(f_APSD)); % to see the frequency response of the filter
figure()
loglog(f_APSD,APSD,'-k',f_APSD,abs(h).^2,'-g')
legend('Target APSD','Filter Response')
%}
g = APSD(1)/(2*pi*f_APSD(1));
h = tf(g*[1 0],1)*tf(1,[1/(2*pi*f_b) 1])*tf(1,[1/(2*pi*f_c) 1])^2;
h = tf(g*[1 0],1)*tf(1,[1/(2*pi*f_c*0.85) 1])^3;
zpk(h)

ans =
 
  6.5822e06 s
  -----------
  (s+1602)^3
 
Continuous-time zero/pole/gain model.

figure
semilogx(f_APSD,20*log10(APSD),f_APSD,db(squeeze(bode(h,2*pi*f_APSD)))),grid
xlabel('Hz');ylabel('dB')

Not sure how good the match needs to be, or if it's really needed to match sqrt(APSD) as I don't fully understand why this approximation helps to evaluate the RMS of the output in response to some input. If the input and the system are known, shouldn't it be possible to find the RMS of the output from that information alone?

[Fabio Cavagna]

I don't understand exactly what you've done.

I'll try to explain better my problem. I have a spring-mass-damper system which is excited at its base by an acceleration. This acceleration is expressed through a power spectral density (APSD), I don't have the time history. I want to evaluate the root mean square (rms) of the acceleration of the mass. One possible solution is: I integrate the APSD multiplied by the power of the absolute value of the transmissibility and then I can take the square root. Another solution is the lyapunov equation, from Matlab help: A*X + X*A' + Q = 0; where A is the state-space form matrix A of the system, Q is B*W*B' with B the ss form matrix B and W is the power spectral density of the input which has to be a white noise (so with a flat power spectral density). The solution of this equation is the variance of the state of the system, with the ss matrix C I can evaluate the variance of the output as C*X*C' and then taking the square root, the rms ot the output. To use this method, as I said, I must have a system in which the input is a white noise and the output can be the rms of the mass acceleration. I can imagine this system as a series of a shaping filter (input white noise, output APSD) and the spring-mass-damper system (input APSD, output mass rms). So that the augmented system has input white noise and output mass rms. What I need is to design the shaping filter so that its magnitude response is the square root of the power spectral density (S_uu(om) = abs(H(j*om))^2 * 1, where 1 is assumed to be the PSD of the white noise, S_uu the APSD and H the filter transfer function), and I need its coefficients to write the transfer function and its state-space form (needed to write the state-space form of the augmented system, whose matrices will be used in the lyapunov equation).

[Paul]

I just constructed a transfer function with frequency response equivalent to APSD.

But, to use this Lyapunov equation approach, I can see why you'd need a transfer function fit to sqrt(APSD).

Can you provided a link to a reference for the Lyapunov equation approach?

Are you obligated to use the Lyapunov equation approach? If not, we know that the PSD of the output is the product of the PSD of the input and the magnitude squared of the transfer function, both of which you have. Then the RMS of the output should be the sqrt of the integral of the PSD of the output, shouldn't it?

I'm going to try an example.

f_APSD and APSD as defined above

df = 0.1;
f_a = 20; f_b = 100; f_c = 300; f_d = 2000; % freq points
s_l = 3; s_r = -5; % slopes
APSD_m = 1; % flat mid PSD
f_l = f_a:df:f_b; f_m = f_b+df:df:f_c; f_r = f_c+df:df:f_d;
APSD_l = APSD_m * (f_l/f_b) .^ (s_l/(10*log10(2))); % left PSD
APSD_r = APSD_m * (f_r/f_c) .^ (s_r/(10*log10(2))); % right PSD
APSD_m = APSD_m * ones(1,length(f_m)); % mid PSD
f_APSD = [f_l,f_m,f_r];
APSD = [APSD_l,APSD_m,APSD_r];

Define simple second order model

wn    = 2*pi*200;
zeta  = 0.3;

Second order system

h = tf(wn^2,[1 2*zeta*wn wn^2]);
sys = ss(h);

Desired figure of merit, which I think is RMS^2, based on the Lyapunov solution using the directions above

W = 1;
X = lyap(sys.A,sys.B*W*sys.B');
V = sys.C*X*sys.C.'

V = 1.0472e+03

The PSD of the output in response to an input with constant, unit, PSD is abs(H(f))^2. The figure of merit is computed by integrating the PSD of the output and mutiplying by 2, becasue we are only integrating over positive frequencies.

Vs = 2*integral(@(f) reshape(abs(freqresp(sys,2*pi*f)).^2,size(f)),0,inf)

Vs = 1.0472e+03

For the real problem at hand, define the transfer function from acceleration of the base to acceleration of the mass

h = tf([2*zeta*wn wn^2],[1 2*zeta*wn wn^2]);
sys = ss(h);

Overlay the frequency response of the system with the PSD of the input

figure
semilogx(f_APSD,20*log10(APSD),f_APSD,db(squeeze(bode(sys,2*pi*f_APSD)))),grid
legend('abs(H(f))','APSD')

Compute the integral of product of APSD and abs(H(f))^2, only over the limits of frequency as defined for APSD. Linear interpolation seems reasonable given that df seems small.

2*integral(@(f) interp1(f_APSD,APSD,f,'linear').*reshape(abs(freqresp(sys,2*pi*f)).^2,size(f)),f_APSD(1),f_APSD(end))

ans = 1.2210e+03

Extend APSD by asusming it's zero at dc and rolls down to zero linearly at f = 10000. Of course, you can model this more realistically if desired

apsd = @(f) interp1([0 f_APSD 1e4],[0 APSD 0],f,'linear',0);
2*integral(@(f) apsd(f).*reshape(abs(freqresp(sys,2*pi*f)).^2,size(f)),0,inf)

ans = 1.2254e+03

[Fabio Cavagna]

Thanks for the feedback,

As I wrote previously

"One possible solution is: I integrate the APSD multiplied by the power of the absolute value of the transmissibility and then I can take the square root."

Which is essentialy the code you wrote, but I already had this solution.

My problem is: let's imagine that the frequency of the spring-mass-damper system is not a single value but a vector of 10000 values, I should do a for loop and then perform the integration 10000 times! This is why I want to try the lyapunov approach in which basically I define the shaping filter once and then in the for loop I have only lyapunov.

[Paul]

I don't understand this statement: "the frequency of the spring-mass-damper system is not a single value but a vector of 10000 values." Which frequency? The natural frequency of the spring-mass-damper? Does this mean you actually have 10000 different spring-mass-damper systems under consideration?

[Fabio Cavagna]

Yes. 10000 is a fictitious number, it could be 100 or 1000

[Paul]

Specified PSD of base motion

df = 0.1;
f_a = 20; f_b = 100; f_c = 300; f_d = 2000; % freq points
s_l = 3; s_r = -5; % slopes
APSD_m = 1; % flat mid PSD
f_l = f_a:df:f_b; f_m = f_b+df:df:f_c; f_r = f_c+df:df:f_d;
APSD_l = APSD_m * (f_l/f_b) .^ (s_l/(10*log10(2))); % left PSD
APSD_r = APSD_m * (f_r/f_c) .^ (s_r/(10*log10(2))); % right PSD
APSD_m = APSD_m * ones(1,length(f_m)); % mid PSD
f_APSD = [f_l,f_m,f_r];
APSD = [APSD_l,APSD_m,APSD_r];

LTI filter that has a magnitude response of approximately the sqrt(APSD). Might be able to do a little better with a higher order model.

sqrtapsdsys = zpk([-117.9 , -1.899],[-23.9 , -493.4, -2040],2589.6); 

Compare APSD with magnitude squared of the approximating filter

figure
semilogx(f_APSD,20*log10(APSD),f_APSD,db(squeeze(bode(sqrtapsdsys,2*pi*f_APSD)).^2)),grid
legend('APSD','sqrtapsd^2')

Plant model for testing. Base acceleration is the input, mass acceleration is the output.

wn    = 2*pi*200;
zeta  = 0.3;

h = tf([2*zeta*wn wn^2],[1 2*zeta*wn wn^2]);
sys = ss(h);

Augment the plant with the pre-filter.

augsys = sys*sqrtapsdsys;

Compute the firgure of merit using lyap

W = 1;
X = lyap(augsys.A,augsys.B*W*augsys.B');
V = augsys.C*X*augsys.C.'

V = 1.2229e+03

Compute the figure of merit by direct integration of the product of APSD and abs(sys)^2. APSD approximated past the specified boundaries. Would be better to properly extend APSD.

apsd = @(f) interp1([0 f_APSD 1e4],[0 APSD 0],f,'linear',0);
2*integral(@(f) apsd(f).*reshape(abs(freqresp(sys,2*pi*f)).^2,size(f)),0,inf)

ans = 1.2254e+03

Resuls are pretty close.

[Fabio Cavagna]

That's amazing, thank you very much. Can I ask you how you defined the filter coefficients? Are there any special criteria? So if I wanted to increase the order of the filter or try with another APSD specification I would know how to do it.

[Paul]

Probably could have got something close to those coefficients by hand with a bit of trial and error. What I actually did was:

1) extend the APSD to 0 on the left and inf on the right by setting f_a = 0 and f_d = inf.

2) estimated the phase of sqrt(APSD) using Bode's gain-phase relationship.

3) with the amplitude and phase of sqrt(APSD), I then used invfreqs with iter=30 and unit weights to get a transfer function with two zeros and four poles. I used a bit of trial and error on m and n to get the lowest order transfer function with magnitude that still seemed like a good fit to sqrt(APSD). I got a lot of warnings from invfreqs

4. The result from invfreqs had a very high frequency pole, so I used freqsep to get rid of it, leaving me with the 2-zero, 3-pole form above.

[Fabio Cavagna]

Could you provide the code for it?

[Paul]

If asking for the code for step 2, I'm afraid not. It's still relatively fresh, etc.

The form of Bode's gain-phase theorem that I used is:

For a minimum phase transfer function G(s), the phase of G(j*w) at frequency w0 is

where

I used a two-sided difference to estimate the the derivative of M(u)

[Paul]

Here is a complete code. It's not the same code that I used in the comment above, but gives a very similar result.

The key piece is to esimate the phase of a SISO system based only on the gain. This code uses the methodology from Relationship of magnitude response to phase response. It seems to work for this problem and simple cases, but I certainly haven't tested it a lot.

Here's an example of estimating the phase from the magnitude for a simple problem.

warning('off'); % turn warnings off to not clutter up the output.  Lot's of warnings!
h = tf(1,[1 1]);h = h*tf([1 5],[1 100]);
w = logspace(-2,3,1000);
Hw = squeeze(freqresp(h,w));
w0 = w(1:20:end);
figure
semilogx(w,180/pi*angle(Hw),w0,180/pi*phaseest(@(w) reshape(abs(freqresp(h,w)),size(w)),w0),'o')
xlabel('rad/sec'); ylabel('deg')

Define the desired APSD

df = 0.1;
f_a = 20; f_b = 100; f_c = 300; f_d = 2000; % freq points
s_l = 3; s_r = -5; % slopes
APSD_m = 1; % flat mid PSD
f_l = f_a:df:f_b; f_m = f_b+df:df:f_c; f_r = f_c+df:df:f_d;
APSD_l = APSD_m * (f_l/f_b) .^ (s_l/(10*log10(2))); % left PSD
APSD_r = APSD_m * (f_r/f_c) .^ (s_r/(10*log10(2))); % right PSD
APSD_m = APSD_m * ones(1,length(f_m)); % mid PSD
f_APSD = [f_l,f_m,f_r];
APSD = [APSD_l,APSD_m,APSD_r];

Estimate the phase of the desired pre-filter, which has magnitude equal to the sqrt of the magnitude of APSD.

Extend f_APSD the left to get a better estimate at low frequency.

f0 = [0.1 f_APSD(1:100:end)];
phasesqrtapsd2  = phaseest(@(w) (sqrt(apsdcalc(w/2/pi))),2*pi*f0);
magsqrtapsd2    = sqrt(apsdcalc(f0));

Use invfreqs to estimate the transfer function.

[b2,a2] = invfreqs(magsqrtapsd2.*exp(1j*phasesqrtapsd2),2*pi*f0,2,4,0*f0+1,30);

Convert to zpk form. There is a poles at very high frequency

sqrtapsdsys = zpk(tf(b2,a2))

sqrtapsdsys =
 
        1.3257e20 (s+201.6) (s+2.942)
  -----------------------------------------
  (s+5.274e16) (s+1885) (s+610.3) (s+61.26)
 
Continuous-time zero/pole/gain model.

Eliminate the pole at high frequency. Maybe should use invfreqs to specify third order denominator.

sqrtapsdsys = freqsep(sqrtapsdsys,2000)

sqrtapsdsys =
 
   2513.4 (s+201.6) (s+2.942)
  ----------------------------
  (s+61.26) (s+610.3) (s+1885)
 
Continuous-time zero/pole/gain model.

Compare APSD with the squared magnitude of the filter, over the frequency range of interest.

figure
semilogx( ...
    f_APSD, db(APSD), ...
    f_APSD, db(abs(freqs(b2,a2,2*pi*f_APSD)).^2),           ...
    f_APSD, db(abs(squeeze(freqresp(sqrtapsdsys,2*pi*f_APSD))).^2) ...
    ) ,grid
xlabel('Hz')
legend('APSD','(b2/a2)^2','sqrtapsdsys^2')

Maybe it would be good to adjust the gain on sqrtapsdsys so that the area under APSD and the area under abs(sqrtapsdsys)^2 over the frequency range of interest is the same?

Plant model for testing. Base acceleration is the input, mass acceleration is the output.

wn    = 2*pi*200;
zeta  = 0.3;
h = tf([2*zeta*wn wn^2],[1 2*zeta*wn wn^2]);
sys = ss(h);

Augment the plant with the pre-filter.

augsys = sys*sqrtapsdsys;

Compute the figure of merit using lyap

W = 1;
X = lyap(augsys.A,augsys.B*W*augsys.B');
V = augsys.C*X*augsys.C.'

V = 1.2081e+03

Compute the figure of merit by direct integration of the product of APSD and abs(sys)^2.

2*integral(@(f) apsdcalc(f).*reshape(abs(freqresp(sys,2*pi*f)).^2,size(f)),0,inf)

ans = 1.2253e+03

warning('on')


function apsd = apsdcalc(f)
% Compute APSD at arbitray frequency. Assumes that APSD rolls off to zero
% using the same equation for f < f_a and f > f_d.
f_a = 20; f_b = 100; f_c = 300; f_d = 2000; % freq points
s_l = 3; s_r = -5; % slopes
APSD_m = 1; 

apsd = 0*f + APSD_m;

APSD_m = 1; % flat mid PSD
apsd(abs(f)<=f_b) = APSD_m * (abs(f(abs(f)<=f_b))/f_b) .^ (s_l/(10*log10(2))); % left PSD
apsd(abs(f)>=f_c) = APSD_m * (abs(f(abs(f)>=f_c)/f_c)) .^ (s_r/(10*log10(2))); % right PSD

apsd(apsd <= 1e-10) = 1e-10;

end

function phase = phaseest(mag,w0)

phase = nan*w0;

for ii = 1:numel(w0)

    phase(ii) = -myhilbert(@(w) log(mag(w)),w0(ii));

end

end

function hi = myhilbert(f,w0)

x = @(w) f(w)./(w0-w);
% hi = integral(x,-inf,inf);
hi = ( integral(x,-inf,w0-eps) + integral(x,w0+eps,inf) )/pi;

end

[Fabio Cavagna]

Thank you very much!