Frequency response of a system operating at a fixed frequency

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.

Answer 1

[Mathieu NOE]

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');

[Mathieu NOE]

now we can improve the code . I introduced a test to see when to stop the identification - when the next point has a magnitude 3 dB below the max modulus

I believ you can easily modify that script to pass data / arguments to and from a simulink file.

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

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;
    
    % convert identified TF real and imaginary part into module / phase
    FRF_mod(k) = sqrt(FRF_real(k).^2 + FRF_imag(k).^2);
    FRF_phase(k) = atan2(FRF_imag(k),FRF_real(k));
    
    % store max value of FRF_mod
    FRF_mod_max = max(FRF_mod_max,max(FRF_mod(k)));
    
    % stop if next value is 3 dB lower than max value
    if abs(FRF_mod(k))<0.707*FRF_mod_max
        k
        break
    end
end

k = 19

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]

Thank you for your answer

The problem I still encounter is that in my case i only have a fft that should be associated to one single frequency, so how do I evaluate the amplitude of the frequency response?

What I mean is that if I excited my sistem only with a single frequency and i perform the fft of the output, what value of this fft should I consider as the response in the amplitude and phase and how should I interpret it?

[Mathieu NOE]

you have also a frequency vector associated with the fft output , so sue that vector and probably some interpolation as your excitation frequency is probably not matching the fft bins

all of this is avoided in my demo

[Alberto Moretti]

i'm sorry what do you mean? as fft output I only have a one dimensional complex vector, am i using fft in the wrong way?

[Mathieu NOE]

no, this is fine

but you want to plot the fft output (as y data) vs. a frequency vector f - like explained in the doc

Fast Fourier transform - MATLAB fft - MathWorks France

example :

Fs = 40;
T = 1/Fs;
L = 211;
t = (0:L-1)*T;
% Create a superposition of a 2 Hz sinusoidal signal and its higher harmonics. 
% The signal contains a 2 Hz cosine wave, a 4 Hz cosine wave, and a 6 Hz sine wave.
 
X = 3*cos(2*pi*2*t) + 2*cos(2*pi*4*t) + sin(2*pi*6*t);

Y = fft(X);
% Compute the single-sided amplitude spectrum of the signal.

f = Fs*(0:(L-1)/2)/L; % <= YOUR FREQUENCY VECTOR
P2 = abs(Y/L);
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);
plot(f,P1,"-o") 
title("Single-Sided Spectrum of Original Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

[Mathieu NOE]

Another solution,if you have the Signal Processing toolbox is to use tfestimate (which is a fft based function but you don't see it if you try to edit the function)

Transfer function estimate - MATLAB tfestimate - MathWorks France

here see how to use it in my last script :

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

samples_ramp_up = 1024;     % ramp up portion of the test signal - no identification is done during this period to avoid transients
samples_ID = 4096-samples_ramp_up;  % samples used for the identification process itself 
samples = samples_ramp_up + samples_ID; % total amount of samples 

amplitude_ID = 0.1; %amplitude of the test signal 

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

FRF_mod_max = 0;


for k = 1:numel(f_id)
    
    fc = f_id(k); % select frequency 
    
    % cos and sin DFT tables (coefficients)
    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 
    
    % generate / measure the process output
    signal_out = filter(nz,dz,signal_ID);
    
    [FRF,Freq] = tfestimate(signal_ID,signal_out,hanning(samples),round(0.75*samples),samples,Fs);
    % interpolation
    FRF_id(k) = interp1(Freq,FRF,fc);
    
    % convert identified TF real and imaginary part into module / phase
    FRF_mod(k) = abs(FRF_id(k));
    FRF_phase(k) = 180/pi*angle(FRF_id(k));
    
    % store max value of FRF_mod
    FRF_mod_max = max(FRF_mod_max,max(FRF_mod(k)));
    
    % stop if next value is 3 dB lower than max value
    if abs(FRF_mod(k))<0.707*FRF_mod_max
        break
        k
    end
    
end


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),'dr');grid on
xlabel('Frequency (Hz)');
ylabel(' phase (°)')
legend('model','ID');

[Alberto Moretti]

thank you very much.

I only have one last doubt

after I compute the frequency vector, due to how it is composed, i get a vector that goes up to 10^6.

I will simply ignore evrything related to frequency different from 0.3 in my case right?

[Mathieu NOE]

ok, you would say that you cannot see how fft is used in "modern" releases of tfestimate ....

now this is my own version of it , where you have access to all the secrets of it (well indeed it's a much oldier version that I changed a bit for my own purpose, but the interesting aspects are now visible).

enjoy !

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

samples_ramp_up = 1024;     % ramp up portion of the test signal - no identification is done during this period to avoid transients
samples_ID = 4096-samples_ramp_up;  % samples used for the identification process itself 
samples = samples_ramp_up + samples_ID; % total amount of samples 

amplitude_ID = 0.1; %amplitude of the test signal 

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

FRF_mod_max = 0;


for k = 1:numel(f_id)
    
    fc = f_id(k); % select frequency 
    
    % cos and sin DFT tables (coefficients)
    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 

    % generate / measure the process output
    signal_out = filter(nz,dz,signal_ID);

%     [FRF,Freq] = tfestimate(signal_ID,signal_out,hanning(samples),round(0.75*samples),samples,Fs);
    [FRF,Cxy,Freq] = mytfe_and_coh(signal_ID,signal_out,samples,Fs,hanning(samples),round(0.75*samples));
    
    % interpolation
    FRF_id(k) = interp1(Freq,FRF,fc);
 
    % convert identified TF real and imaginary part into module / phase
    FRF_mod(k) = abs(FRF_id(k));
    FRF_phase(k) = 180/pi*angle(FRF_id(k));

    % store max value of FRF_mod
    FRF_mod_max = max(FRF_mod_max,max(FRF_mod(k)));
    
    % stop if next value is 3 dB lower than max value
    if abs(FRF_mod(k))<0.707*FRF_mod_max
        break
        k
    end

end


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),'dr');grid on
xlabel('Frequency (Hz)');
ylabel(' phase (°)')
legend('model','ID');

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

%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function  and Coherence Estimate 
    
% compute PSD and CSD
window = window(:);
n = length(x);		% Number of data points
nwind = length(window); % length of window
if n < nwind    % zero-pad x , y if length is less than the window length
    x(nwind)=0;
    y(nwind)=0;  
    n=nwind;
end
x = x(:);		% Make sure x is a column vector
y = y(:);		% Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap));	% Number of windows
					% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;

Pxx = zeros(nfft,1); 
Pyy = zeros(nfft,1); 
Pxy = zeros(nfft,1); 

for i=1:k
    xw = window.*x(index);
    yw = window.*y(index);
    index = index + (nwind - noverlap);
    
    Xx = fft(xw,nfft);
    Yy = fft(yw,nfft);
    Xx2 = abs(Xx).^2;
    Yy2 = abs(Yy).^2;
    Xy2 = Yy.*conj(Xx);
    Pxx = Pxx + Xx2;
    Pyy = Pyy + Yy2;
    Pxy = Pxy + Xy2;
    
end

% Select first half
if ~any(any(imag([x y])~=0))   % if x and y are not complex
    if rem(nfft,2)    % nfft odd
        select = [1:(nfft+1)/2];
    else
        select = [1:nfft/2+1];   % include DC AND Nyquist
    end

    Pxx = Pxx(select);
    Pyy = Pyy(select);
    Pxy = Pxy(select);
    
else
    select = 1:nfft;
end
Txy = Pxy ./ Pxx;                   % transfer function estimate 
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy);    % coherence function estimate 

f = (select - 1)'*Fs/nfft;

end

[Alberto Moretti]

thank you very much.

I only have one last doubt

after I compute the frequency vector, due to how it is composed, i get a vector that goes up to 10^6.

I will simply ignore evrything related to frequency different from 0.3 in my case right?

[Mathieu NOE]

hello again

what is your sampling rate , and how can you obtain a freq vector that goes up to 10^6 (that sound enormous compared to your test signal at 0.3 Hz !!)

again, see my last code , the frequency vector (f) is given by the function mytfe_and_coh (and you pass the sampling rate Fs to the function)

[Txy,Cxy,f] = mytfe_and_coh(x,y,nfft,Fs,window,noverlap)

[Paul]

@Mathieu NOE

With respect to these lines from the original answer

%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;

How is that the phase of signal_ID(ind) can be ignored when computing the FRF? Obviously it all works; I'm trying to understand why it works as well as it does.

I modified the code to explicitly compute output/input in the frequency domain. It turns out FRF_den only has small, but not zero, imaginary part relative to its real part, which I guess is why your code can work the way it does. So I guess my real question is if the code is constructed to ensure that signal_ID(ind) has that small imaginary part. I would expect the ramp up time to be an integer number of periods of the input, but samples_ramp_up is fixed, so I must be missing something.

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
samples_ramp_up = 1024;     % ramp up portion of the test signal - no identification is done during this period to avoid transients
samples_ID = 4096-samples_ramp_up;  % samples used for the identification process itself 
samples = samples_ramp_up + samples_ID; % total amount of samples 
amplitude_ID = 0.1; %amplitude of the test signal 
dt = 1/Fs;
t = (0:samples-1)*dt;
for k = 1:numel(f_id)
    
    fc = f_id(k); % select frequency 
    
    % cos and sin DFT tables (coefficients)
    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 
    % 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;
     FRFnum = sum(signal_out(ind).*exp(-1j*2*pi*fc*t(ind)));
     FRFden = sum(signal_ID(ind) .*exp(-1j*2*pi*fc*t(ind)));
     FRF = FRFnum/FRFden;
     FRF_real(k) = real(FRF);
     FRF_imag(k) = imag(FRF);
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,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,wrapTo180(FRF_phase*180/pi),'dr');grid on
xlabel('Frequency (Hz)');
ylabel(' phase (°)')
legend('model','ID');

[Mathieu NOE]

hello Paul

hope you're doing well

so I owe you some explanations ...

so we assume the TF gain is simply a complex valued scalar at a given f_id value

we can create the test signal using sin or cos functions . To have a common ground with the dft convention (dft(x) = 2/n*sum(x(t)*exp(-jwt)) , the time signals (input and output) must be multiplied by exp(-jwt) = cos(wt) - j*sin(wt) with an integer number of periods. This is nothing new for you

so my convention, as in the dft , is to consider the real part as what is related to cos and the imaginary part to the (-sin)

to make my life easier , I considered a purely "real" input signal so it's solely based on cos function - of course you can do whatever amplitude and phase , but why making things overcomplicated ??

with that "real" only input signal , as we extract the real and imaginary parts of the output signal , we then only have to divide those numbers by the amplitude of the input signal () and here we have the TF real and imaginary parts.

so simple !

then yes you are right, the number of periods must be an integer and I was a bit too quick on that aspect (my very first version of the code was for 1 f_id value so I adapted Fs to have the right number of periods, and I forgot to change the logic when I generalized to the case f_id is now a vector ...oops

I will change that and post the corrected version very soon

[Alberto Moretti]

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.

I'll try your function, thank you very much for your help

[Mathieu NOE]

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 ...

Answer 2

[Paul]

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]

Thank you very much for your answer.

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