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