人工地震波生成程序

function man_made_seismic（）
%%%%%%%%%%%%%%%%%%%%%%%%%%生成人工地震波
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%
%fni=input（‘生成人工地震波－输入文件名：‘‘s‘）;  %%%%%%%%%%%输入文件名
fni=‘seismic.txt‘;
fid=fopen（fni‘r‘）;     %%%打开文件
fs=fscanf（fid‘%f‘1）;  %%%采样频率，人工设置，一般可以设为100，200（因为地震频率一般在0～15HZ）
%fs=100;  %%%采样频率

%%%%%××××××××××××××××××××××××××××××××××××××××
tu=fscanf（fid‘%f‘1）;  %%%上升时间长度
%tu=4.0;
%上升时间包络线线性（1直线；2-抛物线；3-指数曲线）
iu=fscanf（fid‘%f‘1）;
%iu=1; %%直线
%上升时间包络线参数（指数时候需要，其他情况均为1）
cu=fscanf（fid‘%f‘1）; %%%指数形式的变量
ta=fscanf（fid‘%f‘1）; %%%持续时间长度
td=fscanf（fid‘%f‘1）; %%%下降时间长度

%下降时间包络线线性（1直线；2-抛物线；3-指数曲线）
id=fscanf（fid‘%f‘1）;

%下降时间包络线参数（抛物线、指数时候需要，其他情况均为1）
cd=fscanf（fid‘%f‘1） ;%%%
dp=fscanf（fid‘%f‘1）; %%%阻尼比值
p=fscanf（fid‘%f‘1） ;%%%概率系数（一般取0.85）
nn=fscanf（fid‘%d‘1）;;%%迭代次数
fno=fscanf（fid‘%s‘1）;%%输出文件名
x=fscanf（fid‘%f‘[2inf]）;%反应谱频率和幅值数据  （根据设计反应谱计算）
status=fclose（fid）; %%%关闭文件
aaaaa=450.0
%%%%%××××××××××××××××××××××××××××××××××××××××
%plot（x（1:）x（2:））


%%计算地震波总时间长度
tl=tu+ta+td;   %%上升+持续时间+下降


%%%计算生成地震波的数据长度
nt=round（fs*tl+1） %%%采样频率*总持续时间+1
%寻找大于并最接近nt的2的4幂次方为傅立叶变换fft的长度
nfft=2^（nextpow2（nt））  %%%计算与数据长度最近的2的整数次幂
%%%%计算频率间隔（hZ）
df=fs/nfft; %%%%采样频率/fft点数目，没有乘以2PI
%定义反应谱的离散频率向量
f=0:df:（nfft/2.0-1.0）*df;   %f为0~奈归斯特频率的一半（圆频率）
%%计算时间间隔（s）
dt=1.0/fs;  %%%采样频率倒数，时间间隔
%%定义离散时间向量
t=0:dt:（nt-1）*dt;  %时间序列
%%生成0～2*pi的随机数作为随机相位
g=rand（1nfft/2）*2*pi;  %%%随机相位角


%%%建立时间包络线
%%（1）建立与地震波长度相同的1元素向量
en=ones（1nt）;   %%%取值为1的原因是包络线中间段值为1
%%（2）上升时间阶段
    %%%（a）确定上升时间长度
       l=round（tu*fs）+1;  %%上升时间乘采样率  －－－－上升段数据的数目
    %%%产生上升时间段的包络数组
       switch iu
           case 1 %%%直线
                en（1:l）=linspace（01l）;
           case 2 %%抛物线
                a=0:1:l-1;
                en（1:l）=（a/（l-1））.^2;
           case 3 %%指数曲线
                 a=0:1:l-1;
                en（1:l）=1-exp（-cu*a/（l-1））;
       end
       
       
%%（3）持续时间阶段
     %%%确定0时刻到持续时间结束时间段的长度
       m=round（（tu+ta）*fs）+1;
       
%%（4）下降时间阶段
    %%%产生下降阶段的包络线数组元素
 
    switch id       
           case 1 %%%直线
                en（m:nt）=linspace（10nt-m+1）;
           case 2 %%抛物线
                a=0:nt-m;
                en（m:nt）=1-cd*（a*td/（nt-m））.^2;
           case 3 %%指数曲线
                 a=（0:1:nt-m）;    %%%%%%%
                en（m:nt）=exp（-cd*a*td/（nt-m））;
        end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%包络线完成


  
  %%%按照线性插值建立目标反应谱离散数据
  %%%按照目标反应谱的长度生成0元素数组
  a0=zeros（1nfft/2）;   %%%%%%%%%%%%%%%%%%%反应谱点数目由时间总长度和采样率来确定

  %length（a0）
  %%%取目标反应谱的长度
  n=length（x（1:））;    %%%%%%目标反应谱的长度
  %%%四舍五入取整数，求反应谱最小频率对应数组元素的下标（round函数是寻找最接近该数的整数）
  nb=round（x（11）/df）+1;   %%%%反应谱内最小谱值对应的频率
  %%%四舍五入取整数，求反应谱最大频率对应数组元素的下标
  ne=round（x（1n）/df）+1;   %%%%%x存储的为设计反应谱序列，第一列频率，第二列谱值
  
    for k=1:n-1   %%%%n为设计反应谱序列点数目
        %%%（1）四舍五入取整数求反应谱前一个频率数据对应数组元素的下标
        l=round（x（1k）/df）+1;    
        %%%（2）四舍五入取整数求反应谱后一个频率数据对应数组元素的下标
        m=round（