基于FFT由加速度数据实时计算速度与位移的在线实现问题

离线实现（已验证可行） #

我通过FFT方法从加速度数据计算速度和位移，核心逻辑为：

• 对加速度数据做FFT变换
• 构造频率数组Omega
• 应用高通滤波截断低频分量
• 频域内通过除以jω得到速度频域信号，除以-ω²得到位移频域信号
• 逆FFT转换回时域的速度和位移

对应的Matlab离线计算函数：

function [velocity,position]=acc2velpos(acc,dt,f_low)
nsamples = length(acc(:,1));
domega = 2*pi/(dt*nsamples);      
acchat = fft(acc);
% Make frequency array:
Omega = zeros(nsamples,1);

% Create Omega:
if rem(nsamples,2) == 0
nh = nsamples/2;
Omega(1:nh+1) = domega*(0:nh);
Omega(nh+2:nsamples) = domega*(-nh+1:-1);
else
nh = (nsamples-1)/2;
Omega(1:nh+1) = domega*(0:nh) ;
Omega(nh+2:nsamples) = domega*(-nh:-1);
end

% High-pass filter:
n_low=floor(2*pi*f_low/domega);
acchat(1:n_low)=0;
acchat(nsamples-n_low+1:nsamples)=0;

% Multiply by omega^2:
% Skip the measurement in the first time step since it might be inaccurate
poshat(2:nsamples) = -acchat(2:nsamples) ./ Omega(2:nsamples).^2;
velhat(2:nsamples) = acchat(2:nsamples) ./ (Omega(2:nsamples)*1j);

% Inverse Fourier transform:
velocity = real(ifft(velhat));
position=real(ifft(poshat));
% --- End of function ---
end

输入参数：采样频率Fs=256Hz，采样点数N=2560（时长10s），截止频率1Hz。离线计算结果与文献《Identification of a Mechanism's Vibration Velocity and Displacement Based on the Acceleration Measurement》（Gu Mingkun，Lü Zhenhua,2010）一致。

在线实现问题（结果失真） #

我需要将该算法改为在线计算，搭建了Simulink模型，其中Matlab Function模块代码如下：

function [vel,pos, veldumcheck]=acc2velpossimulink(acc)
    dt=10/2560;
    f_low=1;
    %velstore=zeros(2560,2560);
    persistent inputArray;

    % Initialize the array on the first function call
    if isempty(inputArray)
        inputArray = zeros(1,1);
    end

    % Store the input in the array
    inputArray = [inputArray; acc];
    
    if length(inputArray)<1000
        vel=0;
        pos=0;
        veldumcheck=0;
        
    else
        nsamples = length(inputArray)-1;
        domega = 2*pi/(dt*2560);
        acchat = fft(inputArray(2:end));
        
    
    % Make frequency array:
        Omega = zeros(nsamples,1);
    % Create Omega:
        if rem(nsamples,2) == 0
            nh = nsamples/2;
            Omega(1:nh+1) = domega*(0:nh);
            Omega(nh+2:nsamples) = domega*(-nh+1:-1);
        else 
            nh = (nsamples-1)/2;
            Omega(1:nh+1) = domega*(0:nh) ;
            Omega(nh+2:nsamples) = domega*(-nh:-1);
        end
        
        % High-pass filter:
        n_low=floor(2*pi*f_low/domega);
        acchat(1:n_low)=0;
        acchat(nsamples-n_low+1:nsamples)=0;
        
        velhat=zeros(1,nsamples);
        velhat = complex(velhat, zeros(size(velhat), class(velhat)) );
        poshat =zeros(1,nsamples);
        poshat = complex(poshat, zeros(size(poshat), class(poshat)) );
        % Multiply by omega^2:
        % Skip the measurement in the first time step since it might be inaccurate
        poshat(2:nsamples) = -acchat(2:nsamples) ./ Omega(2:nsamples).^2;
        velhat(2:nsamples) = acchat(2:nsamples) ./ (Omega(2:nsamples)*1j);
        

        % Inverse Fourier transform for velocity:
        veldum=zeros(1,nsamples);
        veldum(:) = real(ifft((velhat)));
        veldumcheck=veldum(:);
        vel=veldum(end);
        
        % Inverse Fourier transform for position:
        posdum=zeros(1,nsamples);
        posdum(:) = real(ifft(poshat));
              
        pos=posdum(end);
% --- End of function ---
    end
end

在线实现逻辑：逐时间步采集加速度数据并存储，对当前已采集的所有数据做FFT计算，输出最后一个时间步的速度和位移结果。但结果与离线计算差异极大，且存在异常现象：

• 离线计算10秒数据时，t=4.5s至4.85s的速度为正值；但仅计算前4.75秒数据时，该时间段末尾的速度为负值
• 在线实现输出结果失真且持续为负值

寻求该算法在线计算的可行解决方案。

内容的提问来源于stack exchange，提问作者Ali Caglar Sonmez