function IntFcn_Time2(acc,fs)
n = length(acc);
time = 0:1./fs:(n-1)./fs;
nfft = 2^nextpow2(n); %FFT变换长度,大于并最接近n的2的幂次方为FFT的长度
y = fft(acc,nfft); %FFT变换
% 绘制原始加速度信号的时域波形
figure(1)
subplot(311);
plot(time,acc,'k-');
xlabel('时间(s)');
ylabel('加速度幅值(g)');
title('加速度时域波形');
grid on;
%% 绘制功率密度频谱图,计算机床振动频率
P = y.*conj(y)/nfft; %求取功率密度,其中conj(y)是求y的共轭
f = fs*(0:(nfft/2-1))/nfft; %设定频率变化范围
figure(2)
% subplot(422);
plot(f,P(1:nfft/2),'b-'); %绘制功率频谱图
xlabel('频率(Hz)');
ylabel('功率密度');
title('加速度功率密度频谱图');
grid on;
% pmax = max(P(1:nfft/2));
% id = find(P(1:nfft/2)==pmax);
% fm = f(id)
% subplot(422);
% freq_series=fs*(0:nfft/2-1)/nfft;
% plot(freq_series,abs(y(1:nfft/2))*2/nfft);
%% 梯形积分法位移信号disint和速度信号velint
velint = cumtrapz(time, acc); %利用梯形积分法求数值积分,每一个数对应原来数列之前所用数的积分,与trapz不同
velint = velint - repmat(mean(velint), size(velint,1), 1); %去除直流成分 B=repmat(A,m,n):以A的内容堆叠在(MxN)的矩阵B中
disint = cumtrapz(time, velint);
disint = disint - repmat(mean(disint), size(disint,1), 1);
% 去除速度信号的趋势项及弥补能量损失
velenergy = sqrt(sum(velint.^2));
velint = detrend(velint); % detrend 去除信号中的均值或线性趋势
velreenergy = sqrt(sum(velint.^2));
velint = velint/velreenergy*velenergy; % 此操作是为了弥补去趋势时能量的损失
% 绘制速度时域波形
figure(1)
subplot(312);
plot(time,velint,'b-');
xlabel('时间(s)');
ylabel('速度幅值(m/s)');
title('速度时域波形');
grid on;
% 去除位移信号的趋势项及弥补能量损失
disenergy = sqrt(sum(disint.^2));
disint = detrend(disint);
disreenergy = sqrt(sum(disint.^2));
disint = disint/disreenergy*disenergy;
% 去除位移中的二次项
p = polyfit(time, disint, 2); % polyfit(x,y,n) 用多项式求过已知点的表达式 x——横坐标,y——纵坐标,n——拟合的阶数
disint = disint - polyval(p, time); % polyval(p,x) 返回n次多项式p在x处的值
c2 = 1.0e+003; % 单位变换系数,位移单位为mm
disint = disint*c2;
% 绘制位移时域波形
figure(1)
subplot(313);
plot(time,disint,'b-');
xlabel('时间(s)');
ylabel('位移幅值(mm)');
title('位移时域波形');
grid on;
end
大神,看下这个程序有没有错,我用的就是这个程序。 |