zuixiaoercheng.m

利用matlab采用最小二乘法很好的实现了振动阻尼的数值求解

%leastdamp 
%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 最小二乘求模态阻尼 
%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
%fs 采样频率 
%N——汉宁窗的长度 
%elliplow——低通椭圆滤波器 
%ellipband——带通椭圆滤波器 
%Data——原始信号 
%y——Data低通滤波后得到 
%yy——y加窗卷积，用于去掉零飘 
%y2——由(y-yy)重采样 
%y3——由y2椭圆滤波器滤波 
%ymax——衰减信号的最大值 
%y5,-y5——最小二乘法拟合的曲线 
%eta——对数阻尼（模态阻尼系数） 
% 
%————————————————————— 
%需要变化的系数 
%% f 信号卓越频率（拉索被激励的频率） 
%load() 
%startp——信号开始衰减的起点 
%j—— 信号衰减过程持续的周期数 
%—————————————————————— 
clear 
N=800; 
fs=200; 
f=0.5*3; 
%________________ 
%load('x1'); 
%Data=xs'; 
%_________________ 
load('F:\Matlab.dat.21{post.content}1-2.mat'); 
i=15800:27900;Data=Data(i); 
[L,M]=size(Data); 
%Data=Data-mean(Data(1:2000)); 
%将原始信号与低通滤波的信号画出 
figure(1) 
plot((1:L)/200,Data) 
hold on 
grid on 
y=elliplow(Data); %低通滤波后的信号 
plot((1:L)/200,y,'r') 
grid on 
title('原始信号及低通滤波后的信号') 
xlabel('时间 /s') 
% 
windom=hanning(N)/(N+1)*2; %加汉宁窗 
yy=conv(y,windom); 
y=y-yy(N/2:N/2+L-1); 
figure(2) 
plot((1:L)/200,y) 
grid on 
% xlabel('Time(sec)'); 
% ylabel('Acceleration(m/s^2)') 
%hold off 
title('加汉宁窗后的信号') 
xlabel('时间 /s') 
% 重抽样，重抽样后的低通滤波信号的采样频率为20，即按1/10重抽样 
i=1:fix(L/10); 
y2=y(10*(i-1)+1)*1.15; 
figure(3) 
plot((i-1)/fs*10,y2) 
title('重抽样后的信号') 
xlabel('时间 /s') 
grid on 
%带通滤波（第一阶低通滤波），采用椭圆滤波器 
y3=ellipband(y2,fs/10,f); 
figure(4) 
plot((i-1)/fs*10,y3) 
title('带通滤波后的信号') 
xlabel('时间 /s') 
grid on 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%求每一周的最大值及相应的位置 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
b=fix(fs/10/f); %重采样后信号一个周期的点数 f=1/T,f为运动的频率，单位Hz 
w=2*pi*f; % w为运动的圆频率或角速度，单位：弧度/秒，f=w/2pi 
startp=30; %%信号开始衰减的起点 
[ymax(1),I]=max(y3(startp:startp+80)); 
x(1)=I+startp-1; 
for j=2:61; %%信号衰减过程持续的周期数 
[ymax(j),I]=max(y3(x(j-1)+b-5:x(j-1)+b+5)); 
x(j)=x(j-1)+b-6+I; 
if x(j)+b+10>L; 
break; 
end 
end 
figure(5) 
hold on 
grid on 
plot(i/20,y3,x/20,ymax,x/20,-ymax) 
xlabel('时间 s') 
ylabel('振幅') 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%最小二乘法求模态阻尼系数 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
[LL,MM]=size(ymax); 
y4=log(ymax); 
t=x/20; 

a(1,1)=MM; 
a(1,2)=sum(t); 
a(2,1)=a(1,2); 
a(2,2)=sum(t.*t); 
d(1)=sum(y4); 
d(2)=sum(y4.*t); 
bb=a\d'; 
A0=exp(bb(1)) 
eta=-bb(2)/w %模态阻尼比 
y5=A0*exp(bb(2)*t); 
figure(6); 
plot(i/20,y3,t,y5,t,-y5) 
xlabel('Time(sec)'); 
ylabel('magnitude(m/s^2 \mm)') 
title(['模态阻尼比\zeta=',num2str(eta*100),'%']) 
[pxx,ff]=psd(y3,2048*8,20); 
figure(7) 
plot(ff,pxx)