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（