用MATLAB拟合生成地震波程序

zntimmon

求教各位大大谁有源程序，或者帮忙指导下如何做，有的话发我邮箱里，谢谢了。zhangfanyu88@yahoo.com.cn

ch_j1985

LZ搜索一下版面：地震波
看看此帖是否有用？
http://forum.vibunion.com/thread-25956-1-1.html

[ 本帖最后由 ch_j1985 于 2008-11-20 21:50 编辑 ]

happy

1. clear
   
2. clc
   
3. close all hidden
   
4. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
5. fni=input('生成人工地震波-输入数据文件名(20041012):','s');
   
6. fid=fopen(fni,'r');
   
7. fs=fscanf(fid,'%f',1);%采样频率
   
8. tu=fscanf(fid,'%f',1);
   
9. %上升时间包络线线形(1-直线、2-抛物线、3-指数曲线)
   
10. iu=fscanf(fid,'%f',1);
    
11. cu=fscanf(fid,'%f',1);
    
12. ta=fscanf(fid,'%f',1);%持时时间长度
    
13. td=fscanf(fid,'%f',1);%下降时间长度
    
14. id= fscanf(fid,'%f',1);
    
15. cd=fscanf(fid,'%f',1);
    
16. dp=fscanf(fid,'%f',1);%阴尼比值
    
17. p=fscanf(fid,'%f',1);%概率系数（一般可取P=0.85）
    
18. nn=fscanf(fid,'%f',1);%迭代次数
    
19. fno=fscanf(fid,'%f',1);%输出数据文件名
    
20. x=fscanf(fid,'%f',[2,inf]);%反应谱频率和幅值数据
    
21. tatus=fclose(fid);
    
22. tl=tu+ta+td;
    
23. nt=round(fs*tl+1);
    
24. nfft=2^nestpow2(nt)
    
25. %计算频率间隔（Hz）
    
26. df=fs/nfft
    
27. %定义反应谱的离散频率向量
    
28. f=0:df:(nfft/2-1)*df
    
29. %计算时间间隔（s）
    
30. dt=1/fs;
    
31. %定义的离散时间向量
    
32. t=0:dt:(nt-1)*dt
    
33. %生成0到2PI的随机数为随机相位
    
34. g=rand(1,nfft/s)*2*pi;
    
35. en=ones(1,nt);
    
36. l=round(tu*fs)+1
    
37. switch iu
    
38. case 1 %直线
    
39. en(1:l)=linspace(0,1,1);% y = linspace(a,b,n) generates a row vector y of n points linearly
    
40. spaced between and including a and b.
    
41. case 2 %抛物线
    
42. a=0:l-1;
    
43. en(1:l)=(a/(l-1)).^2;
    
44. case 3 %指数曲线
    
45. a=0:l-1;
    
46. en(1:l)=1-exp(-cu*a/(l-1));
    
47. end
    
48. m=round((tu+ta)*fs)+1;
    
49. switch id
    
50. case 1 %直线
    
51. en(m:nt)=linspace(1,0,nt-m+1);
    
52. case 2 %抛物线
    
53. a=0:nt-m;
    
54. en(m:nt)=1-cd*(a/(nt-m)).^2;
    
55. case 3 %指数曲线
    
56. a=0:nt-m;
    
57. en(m:nt)=exp(-cd*a/(nt-m));
    
58. end
    
59. a0=zeros(x(1,:));
    
60. n=length(x(1,:));
    
61. nb=round(x(1,n)/df)+1;
    
62. for k=1:n-1
    
63. l=round(x(1,K)/df)+1;
    
64. m=round(x(1,K+1)/df)+1;
    
65. a0(1:m)=linspace(x(2,k),x(2,k+1),m-l+1)
    
66. end
    
67. a1=a0;
    
68. s=zeros(1,nfft/2);
    
69. k=nb:ne;
    
70. s(k)=2*dp/(pi.*(a1(k).^2)./f(k)./(-2*log(-log(p)*pi/tl)./f(k)));
    
71. %将功率谱转换成傅里叶幅值谱
    
72. b1=sqrt(4*df*s)*nfft/2;
    
73. hf=zeros(ne,nfft);
    
74. for j-0:ne-1
    
75. w=2*pi*df*j
    
76. wd=w*sqrt(1-dp*dp);
    
77. e=exp(-t.*W*dp);
    
78. a=t.*wd;
    
79. s=sin(a)>*((1-2*dp*dp)/(1-dp*dp));
    
80. c=cos(a).*(2*dp/sqrt(1-dp*dp));
    
81. h=wd*e.*(s+c)/fs;
    
82. hf(j+1,:)=fft(h,nfft);
    
83. end
    
84. mm=nn
    
85. for k=1:100
    
86. c=b1.*exp(i*g);
    
87. d=[c,c(nfft/2:-1:1)];
    
88. e=ifft(d,nfft);
    
89. y=en.*real(e(1:nt));
    
90. yf=fft(y,nfft);
    
91. for j=1:ne
    
92. d=ifft(yf.*hf(j,:),nfft);
    
93. al(j)=max(real(d(1:nt)));
    
94. end
    
95. if k==mm
    
96. subplot(2,1,1);
    
97. plot(t,y,t,en,t,-en);
    
98. xlabel('时间（s）');
    
99. ylabel('加速度（g）');
    
100. grid on
     
101. sublpot(2,1,2);
     
102. l=1:ne
     
103. plot(f(l),a0(l),':'f(l),a1(l));
     
104. xlabel('频率（Hz）');
     
105. ylabel('加速度（g）');
     
106. legend('目标谱','计算谱');
     
107. grid on;
     
108. ig=input('继续迭代次数[取值1-9,否则退出]:');
     
109. if ig>0&ig<10 %如果输入数字是1-9
     
110. mm=mm+ig
     
111. else
     
112. break;
     
113. end
     
114. end
     
115. c=bl
     
116. j=nb:ne;
     
117. bl(j)=c(j).*a0(j)./a1(j);
     
118. end
     
119. fid=fopen(fno,'W');
     
120. for K=1:nt
     
121. fprintf(fid,'%f%f\n',t(k),y(k));
     
122. end
     
123. status=fclose(fid);

复制代码