对此可以用滤波的方法将大的趋势项去掉。

测试的代码如下

% 测试积分对正弦信号的作用

clc

clear

close all

%% 原始正弦信号

ts = 0.001;

fs = 1/ts;

t = 0:ts:1000*ts;

f = 50;

dis = sin(2*pi*f*t); % 位移

vel = 2*pi*f.*cos(2*pi*f*t); % 速度

acc = -(2*pi*f).^2.*sin(2*pi*f*t); % 加速度

% 多个正弦波的测试

% f1 = 400;

% dis1 = sin(2*pi*f1*t); % 位移

% vel1 = 2*pi*f1.*cos(2*pi*f1*t); % 速度

% acc1 = -(2*pi*f1).^2.*sin(2*pi*f1*t); % 加速度

% dis = dis + dis1;

% vel = vel + vel1;

% acc = acc + acc1;

%

结：频域积分正常恢复信号，时域积分恢复加入的高频信息有误差

% 加噪声测试

acc = acc + (2*pi*f).^2*0.2*randn(size(acc));

% 结：噪声会使积分结果产生大的趋势项

figure

ax(1) = subplot(311);

plot(t, dis), title('位移')

ax(2) = subplot(312);

plot(t, vel), title('速度')

ax(3) = subplot(313);

plot(t, acc), title('加速度')

linkaxes(ax, 'x');

% 由加速度信号积分算位移

[disint, velint] = IntFcn(acc, t, ts, 2);

axes(ax(2)); hold on

plot(t, velint, 'r'), legend({'原始信号',

'恢复信号'})

axes(ax(1)); hold

on

plot(t, disint, 'r'), legend({'原始信号', '恢复信号'})

%% 测试积分算子的频率特性

n = 30;

amp = zeros(n, 1);

f = [5:30 40:10:480];

figure

for i = 1:length(f)

fi = f(i);

acc = -(2*pi*fi).^2.*sin(2*pi*fi*t); % 加速度

[disint, velint] = IntFcn(acc, t, ts, 2); % 积分算位移

amp(i) = sqrt(sum(disint.^2))/sqrt(sum(dis.^2));

plot(t, disint)

drawnow

% pause

end

close

figure

plot(f, amp)

title('位移积分的频率特性曲线')

xlabel('f')

ylabel('单位正弦波的积分位移幅值')

以上代码中使用IntFcn函数实现积分，它是封装之后的函数，可以实现时域积分和频域积分，其代码如下

%

积分操作由加速度求位移，可选时域积分和频域积分

function [disint, velint] = IntFcn(acc, t, ts, flag)

if flag == 1

% 时域积分

[disint, velint] = IntFcn_Time(t, acc);

velenergy = sqrt(sum(velint.^2));

velint = detrend(velint);

velreenergy = sqrt(sum(velint.^2));

velint =

velint/velreenergy*velenergy; disenergy = sqrt(sum(disint.^2));

disint = detrend(disint);

disreenergy = sqrt(sum(disint.^2));

disint = disint/disreenergy*disenergy; %

此操作是为了弥补去趋势时能量的损失

% 去除位移中的二次项

p = polyfit(t, disint, 2);

disint = disint - polyval(p, t);

else

% 频域积分

velint = iomega(acc, ts, 3, 2);

velint = detrend(velint);

disint = iomega(acc, ts, 3, 1);

% 去除位移中的二次项

p = polyfit(t, disint, 2);

disint = disint - polyval(p, t);

end

end

其中时域积分的子函数如下

% 时域内梯形积分

function [xn, vn] = IntFcn_Time(t, an)

vn = cumtrapz(t, an);

vn = vn - repmat(mean(vn), size(vn,1), 1);

xn = cumtrapz(t, vn);

xn = xn - repmat(mean(xn), size(xn,1), 1);

end

频域积分的子函数如下(此代码是一个老外编的，在频域内实现积分和微分操作)

function dataout = iomega(datain, dt,

datain_type, dataout_type)

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

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

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

%

% IOMEGA is a

MATLAB script for converting displacement, velocity, or

% acceleration

time-series to either displacement, velocity, or

% acceleration

times-series. The script takes an array of waveform data

% (datain),

transforms into the frequency-domain in order to more easily

% convert into

desired output form, and then converts back into the time

% domain

resulting in output (dataout) that is converted into the

desired

% form.

%

% Variables:

% ----------

%

% datain = input waveform

data of type datain_type

%

% dataout = output waveform

data of type dataout_type

%

% dt = time increment

(units of seconds per sample)

%

% 1 - Displacement

% datain_type = 2 -

Velocity

% 3 - Acceleration

%

% 1 - Displacement

% dataout_type

= 2 -

Velocity

% 3 - Acceleration

%

%

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

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

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

% Make sure that

datain_type and dataout_type are either 1, 2 or 3

if (datain_type < 1 || datain_type >

3)

error('Value for datain_type must be a 1, 2 or 3');

elseif (dataout_type < 1 || dataout_type

> 3)

error('Value for dataout_type must be a 1, 2 or 3');

end

% Determine

Number of points (next power of 2), frequency increment

% and Nyquist

frequency

N = 2^nextpow2(max(size(datain)));

df = 1/(N*dt);

Nyq = 1/(2*dt);

% Save frequency

array

iomega_array = 1i*2*pi*(-Nyq : df : Nyq-df);

iomega_exp = dataout_type - datain_type;

% Pad datain

array with zeros (if needed)

size1 = size(datain,1);

size2 = size(datain,2);

if (N-size1 ~= 0 && N-size2 ~=

0)

if size1 > size2

datain = vertcat(datain,zeros(N-size1,1));

else

datain = horzcat(datain,zeros(1,N-size2));

end

end

% Transform

datain into frequency domain via FFT and shift output (A)

% so that

zero-frequency amplitude is in the middle of the array

% (instead of the

beginning)

A = fft(datain);

A = fftshift(A);

% Convert datain

of type datain_type to type dataout_type

for j = 1 : N

if iomega_array(j) ~= 0

A(j) = A(j) * (iomega_array(j) ^ iomega_exp);

else

A(j) = complex(0.0,0.0);

end

end

% Shift new

frequency-amplitude array back to MATLAB format and

% transform back

into the time domain via the inverse FFT.

A = ifftshift(A);

datain = ifft(A);

% Remove zeros

that were added to datain in order to pad to next

% biggerst power

of 2 and return dataout.

if size1 > size2

dataout = real(datain(1:size1,size2));

else

dataout = real(datain(size1,1:size2));

end

return