人工地震動生成程序

1、clearclccloseallhidden%fni=input(生成人工地震波-輸入數(shù)據(jù)文件名(20041012):,s);fid=fopen(fni,r);fs=fscanf(fid,%f,1);%采樣頻率%tu=fscanf(fid,%f,1);%上升時間長度%上升時間包絡(luò)線線形(1-直線、2拋物線、3-指數(shù)曲線)iu=fscanf(fid,%f,1);%上升時間包絡(luò)線線形參數(shù)(只有指數(shù)曲線需要具體參數(shù)，其均為1)cu=fscanf(fid,%f,1);ta=fscanf(fid,%f,1);%持時時間長度td=fscanf(fid,%f,1);%下降時間長度%下降時間包絡(luò)線線形(1-直

2、線、2拋物線、3-指數(shù)曲線)id=fscanf(fid,%f,1);%下降時間包絡(luò)線線形(只有拋物線，指數(shù)曲線需要具體參數(shù)，其余為1)cd=fscanf(fid,%f,1);批注MS1:定義強度包絡(luò)線%dp=fscanf(fid,%f,1);%陰尼比值p=fscanf(fid,%f,1);%概率系數(shù)(一般可取P=0.85)nn=fscanf(fid,%f,1);%迭代次數(shù)fno=fscanf(fid,%f,1);%輸岀數(shù)據(jù)文件名%對目標反應(yīng)譜取值x=fscanf(fid,%f,2,inf);%反應(yīng)譜頻率和幅值數(shù)據(jù)%tatus=fclose(fid);%計算生成地震波的數(shù)據(jù)長度tl=tu+ta+

3、td;%計算生成地震波的數(shù)據(jù)長度nt=round(fs*tl+1);%大于并最接近nt的2的幕次方為FFT長度nfft=2Anestpow2(nt)%計算頻率間隔(Hz)df=fs/nfft%定義反應(yīng)譜的離散頻率向量f=0:df:(nfft/2-1)*df%計算時間間隔(s)dt=1/fs;%定義的離散時間向量t=0:dt:(nt-1)*dt%生成0到2P啲隨機數(shù)為隨機相位g=rand(1,nfft/s)*2*pi;%建立時間包絡(luò)線%建立與地震波長度相同元素為1的向量en=ones(1,nt);%上升時間階段%確定上升時間段的長度l=round(tu*fs)+1%產(chǎn)生上升時間段的包絡(luò)線數(shù)組元素

4、switchiucase1%直線en(1:l)=linspace(0,1,1);%y=linspace(a,b,n)generatesarowvectoryofnpointslinearlyspacedbetweenandincludingaandb.case2%拋物線a=0:l-1;en(1:l)=(a/(l-1).A2;case3%指數(shù)曲線a=0:l-1;en(1:l)=1-exp(-cu*a/(l-1);end%持續(xù)時間階段%確定0時刻到持續(xù)時間結(jié)束時刻時間段的長度m=round(tu+ta)*fs)+1;%下降時間階段%產(chǎn)生下降時間段的包絡(luò)線數(shù)組元素switchidcase1%直線en

5、(m:nt)=linspace(1,0,nt-m+1);case2%拋物線a=0:nt-m;en(m:nt)=1-cd*(a/(nt-m).A2;case3%指數(shù)曲線a=0:nt-m;en(m:nt)=exp(-cd*a/(nt-m);end%按線性插值建立目標反應(yīng)譜離散數(shù)據(jù)%按目標反應(yīng)譜的長度生成元素為0的向量a0=zeros(x(1,:);%取目標反應(yīng)譜數(shù)據(jù)的長度n=length(x(1,:);%四舍五入取整求反應(yīng)譜最大頻率對應(yīng)數(shù)組元素的下標nb=round(x(1,n)/df)+1;fork=1:n-1%四舍五入取整求反應(yīng)譜前一個頻率數(shù)據(jù)對應(yīng)數(shù)組元素的下標l=round(x(1,K)/d

6、f)+1;%四舍五入取整求反應(yīng)譜后一個頻率數(shù)據(jù)對應(yīng)數(shù)組元素的下標m=round(x(1,K+1)/df)+1;%線性插值產(chǎn)生前后兩個頻率數(shù)據(jù)間的反應(yīng)譜數(shù)組元素a0(1:m)=linspace(x(2,k),x(2,k+1),m-l+1)end%根據(jù)目標反應(yīng)譜計算對應(yīng)的近似功率譜a1=a0;s=zeros(1,nfft/2);k=nb:ne;s(k)=2*dp/(pi.*(a1(k)A2)./f(k”(-2*log(-log(p)*pi/tl)./f(k);%將功率譜轉(zhuǎn)換成傅里葉幅值譜b1=sqrt(4*df*s)*nfft/2;%定義元素為0的反譜傳遞函數(shù)矩陣hf=zeros(ne,nfft)

7、;%計算加速度反應(yīng)譜傳遞函數(shù)矩陣forj-0:ne-1w=2*pi*df*jwd=w*sqrt(1-dp*dp);e=exp(-t.*W*dp);a=t.*wd;s=sin(a)*(1-2*dp*dp)/(1-dp*dp);c=cos(a).*(2*dp/sqrt(1-dp*dp);%計算加速度反應(yīng)譜的脈沖響應(yīng)函數(shù)向量h=wd*e.*(s+c)/fs;%通過FFT變換求加速度反應(yīng)譜傳遞函數(shù)向量hf(j+1,:)=fft(h,nfft);endmm=nn%進行生成人工地震波迭代計算%100為最大迭代次數(shù)fork=1:100%將幅值譜和相位譜轉(zhuǎn)化為實部和虛部c=b1.*exp(i*g);%將正負圓

8、頻率傅里葉譜向量組合成一仙向量d=c,c(nfft/2:-1:1);%IFFT變換，并取變換結(jié)果實部為生成的地震波e=ifft(d,nfft);%給生成的地震波加上強度包絡(luò)線y=en.*real(e(1:nt);%計算反應(yīng)譜%對生成的地震波進行FFT變換yf=fft(y,nfft);forj=1:ne%用地震波FFT變換結(jié)果和反應(yīng)譜傳遞函數(shù)的乘積的逆變換做卷積運算d=ifft(yf.*hf(j,:),nfft);%求各頻率對應(yīng)地震的最大響應(yīng)al(j)=max(real(d(1:nt);end%如果達到指定的迭代次數(shù)顯于圖形ifk=mmsubplot(2,1,1);%m同時顯示生成的地震波的強度包絡(luò)線plot(t,y,t,en,t,-en);xlabel(時間(s);ylabel(加速度(g);gridonsublpot(2,1,2);%同時顯示期望反應(yīng)譜與反應(yīng)譜計算譜、l=1:neplot(f(l),a0(l),:f(l),a1(l);xlabel(瀕率(Hz);ylabel(加速度(g);legend(目標譜,計算譜);gridon;ig=input(繼續(xù)迭代次數(shù)取值1-9,否則退岀:);ifig0&ig10%如果輸入數(shù)字是1-9mm=mm+igelsebreak;endendc=bl%期望譜與計算譜的比值來修改傅里葉值譜j=nb:ne;bl(j