Good afternoon. I want to write a code to output experimental frequency response for system identification. For verification I use the modeled sinusoid and the response of the modeled system. To exclude the influence of transient process and in the future measurement noise, I want to determine the amplitude and phase of the system output using Fourier analysis. The amplitude is fine, but the phase is a problem.

I tried my own code and the code from Matlab doc. I decided on the code from the video: . I still have a problem with the phase. I watch amplitude and phase at the input frequency respectively. I compared with the results of lsqcurvefit. My code:

close all, clear, clc;

step_t = 0.01;

T = 5;
w = (2*pi)/T;
Am = 2;

t_real = 0:step_t:100;

sin_mass = Am*sin(w*t_real);

Tf = 10;
sys = tf([0 1],[Tf 1])

y = lsim(sys,sin_mass,t_real);

figure()
plot(t_real, sin_mass, t_real, y),grid on
legend('input','output')

[f,F_x_f]=fourier(t_real,y,'sinus');

% plot the magnitude spectrum
figure(2)
stem(f,abs(F_x_f))
grid on
xlim([0,20])
xlabel('Frequency in Hz')
ylabel('Amplitude spectrum of the force in N')

% plot the phase spectrum
figure(3)
stem(f,angle(F_x_f))
grid on
xlim([0,20])
xlabel('Frequency in Hz')
ylabel('Phase spectrum of the force in °')

[M,I] = min(abs(f-(1/T)));

out_prob = F_x_f(I);

Amplitude_est = abs(out_prob)
Phase_est = angle(out_prob)

f_out = @(x,u)x(1)*sin(w*u + x(2));

par0 = [Am,0];

par_out = lsqcurvefit(f_out, par0, t_real, y')

y_est = par_out(1)*sin(w*t_real + par_out(2));

figure()
plot(t_real, y, t_real, y_est),grid on
legend('output','output est')

Function fourier from the video:

function [f,X_f]=fourier(t,x_t,modus)
    % check if mode is set
    if nargin<3
        % set to default value
        modus='pulse';
    end
    
    % Number of values -> scalar
    N=length(t);

    % maximum time (in s) -> scalar
    t_max=t(N);
    % Time step (in s) -> scalar
    t_step=t(2);
    % maximum frequency (in Hz) -> scalar
    f_max=0.5/t_step;
        
    % perform Fourier transform -> vector of length N
    if strcmp(modus,'pulse')
        % the unit of the spectrum is 1/Hz (e.g. V/Hz, A/Hz, ...)
        XX(1:N)=t_max/(N-1)*fft(x_t);        
    elseif strcmp(modus,'sinus')
        % the unit of the spectrum is 1 (e.g. V, A, ...)
        XX(1:N)=2/N*fft(x_t);
    else
        error(['The modus ',modus,' is unknown.'])
    end
    % Meaning:
    % for even N:
    %   XX(1) contains the DC component (real)
    %   XX(2:N/2) contains the positive frequency components in ascending order (usually complex)
    %   XX(N/2+1) contains the highest frequency (real)
    %   XX(N/2+2:end) contains the negative frequency components in decreasing order (usually complex)
    % for odd N:
    %   XX(1) contains the DC component (real)
    %   XX(2:(N+1)/2) contains the positive frequency components in ascending order (usually complex)
    %   XX((N+3)/2:end) contains the negative frequency components in decreasing order (usually complex)   
    
    % check whether N is even or odd
    if mod(N,2)
        % N is odd
        % create frequency spectrum -> row vector
        X_f(1:(N+1)/2)=XX(1:(N+1)/2);
        % create frequency range (in Hz) -> row vector
        f=linspace(0,f_max*(N-1)/N,(N+1)/2);
    else
        % N is even
        % create frequency spectrum -> row vector
        X_f(1:N/2+1)=XX(1:N/2+1);
        % create frequency range (in Hz) -> row vector
        f=linspace(0,f_max,N/2+1);
    end
end