f_duffing.asv

地震系统数据分析（加速度、位移、速度）输入信号的时域分析、输入信号时频变换、输出信号的时域分析、输出信号时频变换等

function f_duffing()
%%%%%%%%%%%%%%%%%%%%%%%%%%%输入信号%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid=fopen('ele1.dat','r');%打开数据
ug=fscanf(fid,'%f',inf);   %读取数据
[n,m]=size(ug)      %读取ele1的行数和列数
t=0:0.02:(n-1)*0.02            %离散的时间点
N=30/0.02      %离散时间点的个数  
figure(1)         
plot(t,ug)
xlabel('Time /s');
ylabel('acceleration');
title('时间加速度曲线');
%%%%%%%%%%%%%%%%%%%%%%%%%%%输入信号的时域分析%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%分析内容（1=自相关）
y=xcorr(ug);    %求自相关函数
figure(2);
m=(length(y)-1)/2    
plot([-m:m],y);
xlabel('Time index n');
ylabel('y(n)=Rxx(ug)');
title('自相关函数');
%%%%%%%%%%%%%%%%%%%%%%%%%%%输入信号时频变换%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=1500;%点数
n=0:N-1;
omega=2*pi/N*n
Ug=fft(ug,N)
figure(3);
plot(omega,abs(Ug));
xlabel('omega(rad)');
ylabel('Ug');
title('输入信号的频谱');
%%%%%%%%%%%%%%%%%%%%%%%%%%%赋初值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=1;c=0.3;q=0.6;%c=0.3;q=0.600;%参数赋值
alpha=0.5;beta=0.25;%Newmark参数
s=0;tn=30
x=zeros(n,1);dx=zeros(n,1);ddx=zeros(n,1);%初始化向量分配存储空间,空向量
x(1)=0.4;dx(1)=0.2;ddx(1)=0;%位移及速度赋初值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%循环计算部分%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%梯形求积%%%%%%%%%%%%%%%%%%%%
    if i==2;s=0;
    elseif i==3;s=dx(1)/tn^q+dx(i-1)/0.02^q;
    else s=dx(1)/tn^q+dx(i-1)/0.02^q;
        for j=2:i-2;
            s=s+2*dx(j)/(tn-j*0.02)^q;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    A=1/beta/0.02^2;
    B=-A*x(i-1)-dx(i-1)/beta/0.02-ddx(i-1)*(1/2/beta-1);
    C=c/gamma(1-q);
    D=x(1)/tn^q+0.02*s/2+dx(i-1)*0.02^(1-q)/(1-q)+(1-alpha)*ddx(i-1)*0.02^(2-q)/(1-q)/(2-q)+alpha*B*0.02^(2-q)/(1-q)/(2-q);
    E=alpha*A*0.02^(2-q)/(1-q)/(2-q);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    f=ug(i)
    p1=1;%三次方项系数
    p2=0;%二次方项系数
    p3=A+C*E-1;%一次项系数
    p4=B+C*D-m*f;
    p=[p1 p2 p3 p4];
    r=roots(p);
    x(i)=r(3);   %位移
    ddx(i)=x(i)*A-x(i-1)*A-dx(i-1)/beta/0.02-(0.5/beta-1)*ddx(i-1); %加速度
    dx(i)=dx(i-1)+(1-alpha)*0.02*ddx(i-1)+alpha*0.02*ddx(i);  
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%输出信号的绘图部分%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%时间历程图%%%%%%%%%%%%%%%%%%%%%%%%
figure(4);
plot(i,ddx);
xlabel('Time /s');
ylabel('Displacement');
title('输出信号时间历程图');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%输出信号的时域分析%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%分析内容（1=自相关，2=互相关）
y1=xcorr(ddx);    %求自相关函数
figure(5);
m=(length(y1)-1)/2    
plot([-m:m],y1);
xlabel('Time index n');
ylabel('y(n)=Rxx(y1)');
title('自相关函数');
y2=xcorr(ddx,ug);  %求互相关函数
figure(6);
m=(length(y2)-1)/2    
plot([-m:m],y2);
xlabel('Time index n');
ylabel('y(n)=Rxx(y2)');
title('互相关函数');
%%%%%%%%%%%%%%%%%%%%%%%%%%%输出信号时频变换%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=1500;%点数
n=0:N-1;
omega=2*pi/N*n
DDx=fft(ddx,N)
figure(7);
plot(omega,real(DDx));
xlabel('omega(rad)');
ylabel('DDx');
title('输出信号的频谱');
%%%%%%%%%%%%%%%%%%%%%%%计算功率谱密度%%%%%%%%%%%%%
%x=fft(DDx,n);
%Pyy=x.*conj(x)/n;f=100*(0:511)/n;
%figure(8);
%semilogy(f,Pyy(1:512));
%title('Power Spectral Density of x(t)');
%xlabel('Frequency (Hz)');
%ylabel('Power Spectral Density');%hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%