function?man_made_seismic（）
%%%%%%%%%%%%%%%%%%%%%%%%%%生成人工地震波
clear
clc
close?all?hidden
%%%%%%%%%%%%%%%%%%%%%%%%%
%fni=input（‘生成人工地震波－輸入文件名：‘‘s‘）;??%%%%%%%%%%%輸入文件名
fni=‘seismic.txt‘;
fid=fopen（fni‘r‘）;?????%%%打開文件
fs=fscanf（fid‘%f‘1）;??%%%采樣頻率，人工設(shè)置，一般可以設(shè)為100，200（因為地震頻率一般在0～15HZ）
%fs=100;??%%%采樣頻率

%%%%%××××××××××××××××××××××××××××××××××××××××
tu=fscanf（fid‘%f‘1）;??%%%上升時間長度
%tu=4.0;
%上升時間包絡(luò)線線性（1直線；2-拋物線；3-指數(shù)曲線）
iu=fscanf（fid‘%f‘1）;
%iu=1;?%%直線
%上升時間包絡(luò)線參數(shù)（指數(shù)時候需要，其他情況均為1）
cu=fscanf（fid‘%f‘1）;?%%%指數(shù)形式的變量
ta=fscanf（fid‘%f‘1）;?%%%持續(xù)時間長度
td=fscanf（fid‘%f‘1）;?%%%下降時間長度

%下降時間包絡(luò)線線性（1直線；2-拋物線；3-指數(shù)曲線）
id=fscanf（fid‘%f‘1）;

%下降時間包絡(luò)線參數(shù)（拋物線、指數(shù)時候需要，其他情況均為1）
cd=fscanf（fid‘%f‘1）?;%%%
dp=fscanf（fid‘%f‘1）;?%%%阻尼比值
p=fscanf（fid‘%f‘1）?;%%%概率系數(shù)（一般取0.85）
nn=fscanf（fid‘%d‘1）;;%%迭代次數(shù)
fno=fscanf（fid‘%s‘1）;%%輸出文件名
x=fscanf（fid‘%f‘[2inf]）;%反應(yīng)譜頻率和幅值數(shù)據(jù)??（根據(jù)設(shè)計反應(yīng)譜計算）
status=fclose（fid）;?%%%關(guān)閉文件
aaaaa=450.0
%%%%%××××××××××××××××××××××××××××××××××××××××
%plot（x（1:）x（2:））


%%計算地震波總時間長度
tl=tu+ta+td;???%%上升+持續(xù)時間+下降


%%%計算生成地震波的數(shù)據(jù)長度
nt=round（fs*tl+1）?%%%采樣頻率*總持續(xù)時間+1
%尋找大于并最接近nt的2的4冪次方為傅立葉變換fft的長度
nfft=2^（nextpow2（nt））??%%%計算與數(shù)據(jù)長度最近的2的整數(shù)次冪
%%%%計算頻率間隔（hZ）
df=fs/nfft;?%%%%采樣頻率/fft點數(shù)目，沒有乘以2PI
%定義反應(yīng)譜的離散頻率向量
f=0:df:（nfft/2.0-1.0）*df;???%f為0~奈歸斯特頻率的一半（圓頻率）
%%計算時間間隔（s）
dt=1.0/fs;??%%%采樣頻率倒數(shù)，時間間隔
%%定義離散時間向量
t=0:dt:（nt-1）*dt;??%時間序列
%%生成0～2*pi的隨機數(shù)作為隨機相位
g=rand（1nfft/2）*2*pi;??%%%隨機相位角


%%%建立時間包絡(luò)線
%%（1）建立與地震波長度相同的1元素向量
en=ones（1nt）;???%%%取值為1的原因是包絡(luò)線中間段值為1
%%（2）上升時間階段
????%%%（a）確定上升時間長度
???????l=round（tu*fs）+1;??%%上升時間乘采樣率??－－－－上升段數(shù)據(jù)的數(shù)目
????%%%產(chǎn)生上升時間段的包絡(luò)數(shù)組
???????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?%%指數(shù)曲線
?????????????????a=0:1:l-1;
????????????????en（1:l）=1-exp（-cu*a/（l-1））;
???????end
???????
???????
%%（3）持續(xù)時間階段
?????%%%確定0時刻到持續(xù)時間結(jié)束時間段的長度
???????m=round（（tu+ta）*fs）+1;
???????
%%（4）下降時間階段
????%%%產(chǎn)生下降階段的包絡(luò)線數(shù)組元素
?
????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?%%指數(shù)曲線
?????????????????a=（0:1:nt-m）;????%%%%%%%
????????????????en（m:nt）=exp（-cd*a*td/（nt-m））;
????????end
??%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%包絡(luò)線完成


??
??%%%按照線性插值建立目標(biāo)反應(yīng)譜離散數(shù)據(jù)
??%%%按照目標(biāo)反應(yīng)譜的長度生成0元素數(shù)組
??a0=zeros（1nfft/2）;???%%%%%%%%%%%%%%%%%%%反應(yīng)譜點數(shù)目由時間總長度和采樣率來確定

??%length（a0）
??%%%取目標(biāo)反應(yīng)譜的長度
??n=length（x（1:））;????%%%%%%目標(biāo)反應(yīng)譜的長度
??%%%四舍五入取整數(shù)，求反應(yīng)譜最小頻率對應(yīng)數(shù)組元素的下標(biāo)（round函數(shù)是尋找最接近該數(shù)的整數(shù)）
??nb=round（x（11）/df）+1;???%%%%反應(yīng)譜內(nèi)最小譜值對應(yīng)的頻率
??%%%四舍五入取整數(shù)，求反應(yīng)譜最大頻率對應(yīng)數(shù)組元素的下標(biāo)
??ne=round（x（1n）/df）+1;???%%%%%x存儲的為設(shè)計反應(yīng)譜序列，第一列頻率，第二列譜值
??
????for?k=1:n-1???%%%%n為設(shè)計反應(yīng)譜序列點數(shù)目
????????%%%（1）四舍五入取整數(shù)求反應(yīng)譜前一個頻率數(shù)據(jù)對應(yīng)數(shù)組元素的下標(biāo)
????????l=round（x（1k）/df）+1;????
????????%%%（2）四舍五入取整數(shù)求反應(yīng)譜后一個頻率數(shù)據(jù)對應(yīng)數(shù)組元素的下標(biāo)
????????m=round（