# Frequency response of a system operating at a fixed frequency

38 views (last 30 days)

Show older comments

Alberto Moretti
on 10 Jun 2024

Commented: Mathieu NOE
on 12 Jun 2024

I am trying to simulate the behavior of a system in simulink. I am providing a sinusoidal input with a fixed frequency (e.g. 0.3 Hz) and I am exporting the results in matlab. I want to verify that at this operating frequency, the frequency response is flat. Once done that I'll increase the frequency and so on untill I reach a decrease in the frequency response of -3dB.

However I am facing dissiculties in evaluating this.

I tried doing something like this:

Y_measured = fft(measured);

Y_target = fft(target);

Y = Y_measured./Y_target;

amplitude = abs(Y);

phase = unwrap(angle(Y(1:N/2+1)));

I then plotted the amplitude and the phase over time.

The fact is that being the system evaluated at a fixed frequency, the frequency response diagram loses a lot of significance and I become only intrested in what would be a specific point of the diagram (I guess?).

So in conclusion:

Is this a correct way of evaluating the frequency response of my system?

What would be an efficien way of plotting and showing the results?

Thank you.

##### 0 Comments

### Accepted Answer

Mathieu NOE
on 10 Jun 2024

Edited: Mathieu NOE
on 11 Jun 2024

hello

we could use the general fft approach , but you could also use the stepped sine identification method. We only need to generate a cos/sin waveforms at the desired frequency, generate a test signal that must be fed to the process being identified then project (à la DFT) the output on the cos/sin signals and from there you get the real and imaginary parts of your process

here a little demo - the process is here a simple butterworth filter

EDIT : updated the code so the ID process is operated on a fixed an round number of periods (as it should be);

the fact that the previous version would give correct results is because the samples amount was relatively high so that the last non complete period would have minimal impact on the result.

But this new version is the one and only correct implementation

ThemeCopy

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% stepped sinus identification

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%% inputs %%%%%%%%%%%%%%%

Fs = 500;

freq = linspace(0,Fs/2.56,300);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%tranfer function (process model)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% process is simulated here with a bandpass butterworth filter

fl = 10; % lower corner freq

fh = 100; % upper corner freq

N = 2; % order

[nz,dz] = butter(N,[fl fh]/(Fs/2)); % "process"

[g,p] = dbode(nz,dz,1/Fs,2*pi*freq);

%% Stepped sine main loop

f_id = (10:5:150); % frequencies at which we want to identify our process

period_ramp_up = 10;

period_ID = 20;

amplitude_ID = 0.1; %amplitude of the test signal

dt = 1/Fs;

FRF_mod_max = 0;

for k = 1:numel(f_id)

fc = f_id(k); % select frequency

% compute samples needed

samples_ramp_up = round(period_ramp_up*Fs/fc); % ramp up portion of the test signal - no identification is done during this period to avoid transients

samples_ID = round(period_ID*Fs/fc); % samples used for the identification process itself

samples = samples_ramp_up + samples_ID; % total amount of samples

% cos and sin DFT tables (coefficients)

t = (0:samples-1)*dt;

tcos = cos(2*pi*fc*t); %

tsin = -sin(2*pi*fc*t); %

w = [linspace(0,1,samples_ramp_up) ones(1,samples_ID)]; % ramp up "window"

signal_ID = amplitude_ID*tcos.*w; % signal with ramp up

% uncomment the line below if you want to have a look at the ID signal

% plot(t(1:samples_ramp_up),signal_ID(1:samples_ramp_up),'r')

% hold on

% plot(t(1+samples_ramp_up:samples),signal_ID(1+samples_ramp_up:samples),'b')

% hold off

% generate / measure the process output

signal_out = filter(nz,dz,signal_ID);

% TF real and imag DFT computation (once ramp up is finished)

ind = (samples_ramp_up+1:samples);

FRF_real(k) = sum(signal_out(ind).*tcos(ind))*2/samples_ID/amplitude_ID;

FRF_imag(k) = sum(signal_out(ind).*tsin(ind))*2/samples_ID/amplitude_ID;

end

% convert identified TF real and imaginary part into module / phase

FRF_mod = sqrt(FRF_real.^2 + FRF_imag.^2);

FRF_phase = atan2(FRF_imag,FRF_real);

figure(1);

subplot(2,1,1),semilogx(freq,20*log10(g),'b',f_id(1:k),20*log10(FRF_mod),'dr');grid on

title(' Stepped sine Identification');

ylabel('dB ');

legend('model','ID');

subplot(2,1,2),semilogx(freq,wrapTo180(p),'b',f_id(1:k),wrapTo180(FRF_phase*180/pi),'dr');grid on

xlabel('Frequency (Hz)');

ylabel(' phase (°)')

legend('model','ID');

##### 14 Comments

Mathieu NOE
on 12 Jun 2024

### More Answers (1)

Paul
on 10 Jun 2024

Hi Alberto,

Conceptually, what you're doing is correct (though there are other alternatives to the input signal). Keep in mind that the only element of Y that you care about is the element that corresponds as closely as possible to the target frequency (0.3 Hz in this example). Or you could interpolate Y to that frequency, or you could zero-pad the FFTs to ensure you get a frequency sample exactly at the target frequency. Another thing to consider is the effect of the startup transient; you may want to only do the FFT processing on the data after the system reaches steady state. And you should make sure that you have enough steady state cycles of data. If the system is nonlinear, you may want to consider checking sensitivity of the results to the amplitude of the input. And there may be other complications is you're system is unstable.

Also may be interested in frestimate, for which there is is also a GUI interface (I think in the Model Linearizer).

Whatever path you take, I suggest first working with a model for which you know what the frequency response should be to help iron out the workflow.

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!