# Filter coefficients from frequency response with invfreqs/invfreqz

11 vues (au cours des 30 derniers jours)
Fabio Cavagna le 16 Mai 2023
Modifié(e) : Paul le 21 Mai 2023
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')
##### 2 commentairesAfficher 1 commentaire plus ancienMasquer 1 commentaire plus ancien
Star Strider le 16 Mai 2023

Connectez-vous pour commenter.

### Réponse acceptée

Paul le 16 Mai 2023
Modifié(e) : Paul le 17 Mai 2023
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?
##### 12 commentairesAfficher 11 commentaires plus anciensMasquer 11 commentaires plus anciens
Fabio Cavagna le 21 Mai 2023
Thank you very much!

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Vibration Analysis dans Help Center et File Exchange

R2019b

### Community Treasure Hunt

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

Start Hunting!