# Frequency response of a system operating at a fixed frequency

33 views (last 30 days)
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.

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');
Alberto Moretti on 12 Jun 2024
My sampling period is 2e-6, this is because in the simulink model i need to compute the electronics of an inverter therefore i need very short steps.
Mathieu NOE on 12 Jun 2024
maybe then you could downsample the data to a much more appropriate rate (use "rate converter" or a "ZOH" block to transition to a lower sampling rate) - otherwise you gonna handle very large time data and make a lot of unnecessary computations ...

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).
And might be interested in this discussion.
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.
Alberto Moretti on 10 Jun 2024
I'd like to clarify a few things:
At the moment i am using a sampling rate of 1/(2e-6) = 500000 Hz, and since I am passing the datas from simulink with the same sampling rate, i have vectors of length n*500000 where n is the number of seconds of the simulation. So my frequency resolution is basically 1/n and you are suggesting to make my data vector bigger (by zero padding but could also be done by increasing the time of simulation am i right?) in order to cover exactly the frequency i am intrested in?
I'd also like to be sure on how i would then define the frequency vector that I will associate to the FFT (the Y of my code). Then once plotted the Y amplitude (that I can compute in dB without committing any mistake (?)) and phase with the corresponding frequency vector I will only look at the value in f = 0.3 right? And if the system can follow the input without problems should I expect an amplitude of 0? whilst for the phase?
Thank you again