系统辨识中傅里叶分析计算相位响应的问题求助

我想要编写代码输出用于系统辨识的实验频率响应。为验证效果，我采用建模正弦信号和建模系统的响应进行测试。为排除暂态过程及未来测量噪声的影响，希望通过傅里叶分析确定系统输出的幅值和相位。幅值计算正常，但相位存在问题。

我试过自己编写的代码以及Matlab文档中的代码，最终采用了某视频中的代码，但相位问题仍未解决。我分别查看输入频率下的幅值和相位，并与lsqcurvefit的结果进行对比。我的代码如下：

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

视频中的fourier函数如下：

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

内容的提问来源于stack exchange，提问作者Ilmir Gogorev