求把地震波时程转换成反应谱的matlab 程序

zxwflyer

那位高手有没有把地震波时程转换成反应谱的matlab 程序啊？

[ 本帖最后由 eight 于 2008-4-28 20:04 编辑 ]

luwilero

与楼主共同期待。
我也找了好久了。

jeruselem

1. clear;close all;
   
2. load northbridge.mat %输入地震波文件，时程数据格式为一行。
   
3. eqw=x;
   
4. amax=max(abs(eqw));
   
5. dt=0.01  %采样频率100Hz
   
6. ff=eqw;
   
7. n=0;
   
8. zeta=0.05; %阻尼比
   
9. for w=200:-0.1:1; %自己设定
   
10. x=0;
    
11. v=0;
    
12. %alhp=35;
    
13. acc=0;
    
14. a=[1,2*zeta*w,w^2];
    
15. b=1;
    
16. [r,p]=residue(b,a);
    
17. t=0:dt:length(eqw)*dt;
    
18. h=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t);
    
19. x=conv(h,ff)*dt;
    
20. v=diff(x)/dt;
    
21. acc=diff(v)/dt;
    
22. n=n+1;
    
23. beta1(n)=max(abs(acc))/amax;
    
24. beta2(n)=max(abs(v))/amax;
    
25. beta3(n)=max(abs(x))/amax;
    
26. T(n)=2*pi/w;
    
27. end;
    
28. save eqw_response_Northbridge T beta1 beta2 beta3
    
29. %%
    
30. ww=1./T;wf=fliplr(ww);
    
31. beta1f=fliplr(beta1);
    
32. beta2f=fliplr(beta2);
    
33. beta3f=fliplr(beta3);
    
34. figure;plot(wf,beta1f,'linewidth',2);grid on;xlabel('frequency (Hz) ');ylabel('acceleration (m/s^2)');title('Acceleration Response Spectrum-Northbridge 35gal ');
    
35. saveas(gcf,'Northbridge-Response-1','fig');saveas(gcf,'Northbridge-Response-1','bmp');
    
36. figure;plot(wf,beta2f,'linewidth',2);grid on;xlabel('frequency (Hz)');ylabel('velocity (m/s)');title('Velocity Response Spectrum-Northbridge 35gal')
    
37. saveas(gcf,'Northbridge-Response-2','fig');saveas(gcf,'Northbridge-Response-2','bmp');
    
38. figure;plot(wf,beta3f,'linewidth',2);grid on;xlabel('frequency (Hz)');ylabel('displacement (m)');title('Displacement Response Spectrum-Northbridge 35gal')
    
39. saveas(gcf,'Northbridge-Response-3','fig');saveas(gcf,'Northbridge-Response-3','bmp');
    
40. figure;plot(T,beta1,'linewidth',2);grid on;xlabel('period (s)');ylabel('acceleration (m/s^2)');title('Acceleration Response Spectrum-Northbridge 35gal ');
    
41. saveas(gcf,'Northbridge-Response-4','fig');saveas(gcf,'Northbridge-Response-4','bmp');
    
42. figure;plot(T,beta2,'linewidth',2);grid on;xlabel('period (s)');ylabel('velocity (m/s)');title('Velocity Response Spectrum-Northbridge 35gal ');
    
43. saveas(gcf,'Northbridge-Response-5','fig');saveas(gcf,'Northbridge-Response-5','bmp');
    
44. figure;plot(T,beta3,'linewidth',2);grid on;xlabel('period (s)');ylabel('displacement (m)');title('Displacement Response Spectrum-Northbridge 35gal ');
    
45. saveas(gcf,'Northbridge-Response-6','fig');saveas(gcf,'Northbridge-Response-6','bmp');

复制代码

liujiashun

谢谢啊 好东西

liujiashun

谢谢啊 呵呵 真是高手！

lihong0880

不明白northbridge.mat那个文件的含义，是一行数据表示时间呢，还是两行，一行时刻，一行对应加速度幅值呢

cuidazhe

1. %     ***********************************************
   
2. %     *** Response Spectrum By Duhamel Integral & ***
   
3. %     ***   Newmark Linear Acceleration Method    ***
   
4. %     ***         Author:Yahong DENG (PHD)        ***
   
5. %     ***        Email:hoverdyh@zju.edu.cn        ***
   
6. %     ***           Chang'an University           ***
   
7. %     ***            2008.12.11  Xi'an            ***
   
8. %     ***********************************************
   
9. clc;
   
10. clear all;
    
11. close all hidden;
    
12. disp('***********************************************')
    
13. disp('*** Response Spectrum By Duhamel Integral & ***')
    
14. disp('***       Linear Acceleration Method        ***')
    
15. disp('***        Author:Yahong DENG (PHD)         ***')
    
16. disp('***       Email:hoverdyh@zju.edu.cn         ***')
    
17. disp('***          Chang An University            ***')
    
18. disp('***           2008.12.11  XiAn              ***')
    
19. disp('***********************************************')
    
20. disp('  ')
    
21. disp('-------------------- ')
    
22. disp('**  振动信号输入  ** ')
    
23. disp('-------------------- ')
    
24. disp('  ')
    
25. filenin=input('请输入数据文件名(*.txt/*.dat/*.*)：','s');
    
26. fileid=fopen(filenin,'r');
    
27. sampfre=input('请输入数据采样频率：');
    
28. linsign=input('请输入数据起始行：');
    
29. for line=1:linsign-1
    
30.     tline=fgetl(fileid);
    
31. end
    
32. filesty=input('请选择数据列类型（1：全部(单列/多列) 2：多列中一列）：');
    
33. if filesty==1
    
34.     vibsign=fscanf(fileid,'%f',inf);
    
35.     vibsign=vibsign';
    
36. else
    
37.     colnumb=input('请输入文件总列数：');
    
38.     colsign=input('请输入数据所在列数：');
    
39.     inpsign=fscanf(fileid,'%f',[colnumb,inf]);
    
40.     vibsign=inpsign(colsign,:);
    
41. end
    
42. signumb=length(vibsign);
    
43. time=(0:1/sampfre:(signumb-1)/sampfre);
    
44. maxsign=0;
    
45. for n=1:signumb
    
46.     if abs(vibsign(n))>maxsign
    
47.         maxsign=abs(vibsign(n));
    
48.         maxtime=time(n);
    
49.     end
    
50. end
    
51. disp(' ');
    
52. fprintf('---* 输入信号基本信息 *---\n');
    
53. fprintf('输入信号长度为：%d\n',signumb);
    
54. fprintf('输入信号持时为：%.3f (s)\n',(signumb-1)/sampfre);
    
55. fprintf('输入信号最大值（绝对值）：%f\n',maxsign);
    
56. fprintf('输入信号最大值时间点：%.3f (s)\n',maxtime);
    
57. fprintf('--- 恭喜，成功读入数据！！！ ---\n');
    
58. filesta=fclose(fileid);
    
59. disp('  ')
    
60. disp('---------------------- ')
    
61. disp('**  反应谱参数设置  ** ')
    
62. disp('---------------------- ')
    
63. disp('  ')
    
64. damnumb=input('请输入反应谱阻尼比个数（1-3）：');
    
65. switch damnumb
    
66.     case 1
    
67.         damprat(1)=input('请输入阻尼比值（0-1）：');
    
68.     case 2
    
69.         damprat(1)=input('请输入第一个阻尼比值（0-1）：');
    
70.         damprat(2)=input('请输入第二个阻尼比值（0-1）：');
    
71.     case 3
    
72.         damprat(1)=input('请输入第一个阻尼比值（0-1）：');
    
73.         damprat(2)=input('请输入第二个阻尼比值（0-1）：');
    
74.         damprat(3)=input('请输入第三个阻尼比值（0-1）：');
    
75. end
    
76. haxisty=input('请选择反应谱计算点类型（1：等间距频率计算点 2：等间距周期计算点）：');
    
77. switch haxisty
    
78.     case 1
    
79.         minfreq=input('请输入最小计算频率（一般可取0.1，且 > 0）：');
    
80.         maxfreq=input('请输入最大计算频率（一般可取10-15）：');
    
81.         delfreq=input('请输入计算频率间隔（一般可取0.01/0.02）：');
    
82.         spenumb=round((maxfreq-minfreq)/delfreq)+1;
    
83.         freoper=(minfreq:delfreq:(spenumb-1)*delfreq+minfreq);
    
84.     case 2
    
85.         minperi=input('请输入最小计算周期（一般可取0.1，且 > 0）：');
    
86.         maxperi=input('请输入最大计算周期（一般可取2-6）：');
    
87.         delperi=input('请输入计算周期间隔（一般可取0.01/0.02）：');
    
88.         spenumb=round((maxperi-minperi)/delperi)+1;
    
89.         freoper=(minperi:delperi:(spenumb-1)*delperi+minperi);
    
90. end
    
91. disp('  ')
    
92. disp('-------------------- ')
    
93. disp('**   反应谱计算   ** ')
    
94. disp('-------------------- ')
    
95. disp('  ')
    
96. compmet=input('请选择反应谱计算方法（1：杜哈姆积分卷积法(速度慢) 2：杜哈姆积分FFT法 3：纽马克线性加速度法(速度快,推荐)）：');
    
97. oritime=cputime;
    
98. switch compmet
    
99.     case 1
    
100.         if haxisty==1
     
101.             for n=1:spenumb
     
102.                 angfreq=2*pi*freoper(n);
     
103.                 for m=1:damnumb
     
104.                     angfred(m)=angfreq*sqrt(1-damprat(m)^2);
     
105.                     termexp=exp(-damprat(m)*angfreq*time);
     
106.                     termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2));
     
107.                     termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2)));
     
108.                     uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre;
     
109.                     uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre;
     
110.                     uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre;
     
111.                     dispspe(m,n)=max(abs(conv(vibsign,uniresd)));
     
112.                     velospe(m,n)=max(abs(conv(vibsign,uniresv)));
     
113.                     accespe(m,n)=max(abs(conv(vibsign,uniresa)));
     
114.                     accspeb(m,n)=accespe(m,n)/maxsign;
     
115.                     fftnumb=2^nextpow2(2*signumb);
     
116.                     uniafft=fft(uniresa,fftnumb);
     
117.                     signfft=fft(vibsign,fftnumb);
     
118.                     accespe1=real(ifft(uniafft.*signfft,fftnumb));
     
119.                     accespe2=accespe1(1:signumb);
     
120.                     accespe3=max(accespe2);
     
121.                 end
     
122.                 fprintf('共 %4d 个频率计算点，正在计算第  %4d  个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
     
123.             end
     
124.             disp(' ')
     
125.             elatime=cputime-oritime;
     
126.             fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);     
     
127.         else
     
128.             for n=1:spenumb
     
129.                 angfreq=2*pi/freoper(n);
     
130.                 for m=1:damnumb
     
131.                     angfred(m)=angfreq*sqrt(1-damprat(m)^2);
     
132.                     termexp=exp(-damprat(m)*angfreq*time);
     
133.                     termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2));
     
134.                     termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2)));
     
135.                     uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre;
     
136.                     uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre;
     
137.                     uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre;
     
138.                     dispspe(m,n)=max(abs(conv(vibsign,uniresd)));
     
139.                     velospe(m,n)=max(abs(conv(vibsign,uniresv)));
     
140.                     accespe(m,n)=max(abs(conv(vibsign,uniresa)));
     
141.                     accspeb(m,n)=accespe(m,n)/maxsign;
     
142.                 end   
     
143.             fprintf('共 %4d 个周期计算点，正在计算第   %4d  个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
     
144.             end
     
145.             disp(' ')
     
146.             elatime=cputime-oritime;
     
147.             fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);
     
148.         end
     
149.     case 2
     
150.         fftnumb=2^nextpow2(2*signumb-1);
     
151.         if haxisty==1
     
152.             for n=1:spenumb
     
153.                 angfreq=2*pi*freoper(n);
     
154.                 for m=1:damnumb
     
155.                     angfred(m)=angfreq*sqrt(1-damprat(m)^2);
     
156.                     termexp=exp(-damprat(m)*angfreq*time);
     
157.                     termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2));
     
158.                     termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2)));
     
159.                     uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre;
     
160.                     uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre;
     
161.                     uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre;
     
162.                     fftunid=fft(uniresd,fftnumb);
     
163.                     fftuniv=fft(uniresv,fftnumb);
     
164.                     fftunia=fft(uniresa,fftnumb);
     
165.                     fftsign=fft(vibsign,fftnumb);
     
166.                     dispspt=real(ifft(fftunid.*fftsign,fftnumb));
     
167.                     accespe(m,n)=max(dispspt(1:signumb));
     
168.                     velospt=real(ifft(fftuniv.*fftsign,fftnumb));
     
169.                     velospe(m,n)=max(velospt(1:signumb));
     
170.                     accespt=real(ifft(fftunia.*fftsign,fftnumb));
     
171.                     accespe(m,n)=max(accespt(1:signumb));
     
172.                     accspeb(m,n)=accespe(m,n)/maxsign;
     
173.                 end
     
174.                 fprintf('共 %4d 个频率计算点，正在计算第  %4d  个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
     
175.             end
     
176.             disp(' ')
     
177.             elatime=cputime-oritime;
     
178.             fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);     
     
179.         else
     
180.             for n=1:spenumb
     
181.                 angfreq=2*pi/freoper(n);
     
182.                 for m=1:damnumb
     
183.                     angfred(m)=angfreq*sqrt(1-damprat(m)^2);
     
184.                     termexp=exp(-damprat(m)*angfreq*time);
     
185.                     termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2));
     
186.                     termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2)));
     
187.                     uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre;
     
188.                     uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre;
     
189.                     uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre;
     
190.                     fftunid=fft(uniresd,fftnumb);
     
191.                     fftuniv=fft(uniresv,fftnumb);
     
192.                     fftunia=fft(uniresa,fftnumb);
     
193.                     fftsign=fft(vibsign,fftnumb);
     
194.                     dispspt=real(ifft(fftunid.*fftsign,fftnumb));
     
195.                     accespe(m,n)=max(dispspt(1:signumb));
     
196.                     velospt=real(ifft(fftuniv.*fftsign,fftnumb));
     
197.                     velospe(m,n)=max(velospt(1:signumb));
     
198.                     accespt=real(ifft(fftunia.*fftsign,fftnumb));
     
199.                     accespe(m,n)=max(accespt(1:signumb));
     
200.                     accspeb(m,n)=accespe(m,n)/maxsign;
     
201.                 end   
     
202.             fprintf('共 %4d 个周期计算点，正在计算第   %4d  个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
     
203.             end
     
204.             disp(' ')
     
205.             elatime=cputime-oritime;
     
206.             fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);
     
207.         end
     
208.     case 3
     
209.         if haxisty==1
     
210.             for n=1:spenumb
     
211.                 angfreq=2*pi*freoper(n);
     
212.                 deltime=1/sampfre;
     
213.                 for m=1:damnumb
     
214.                     a(m)=angfreq*sqrt(1-damprat(m)^2)*deltime;
     
215.                     b(m)=damprat(m)*angfreq*deltime;
     
216.                     d(m)=sqrt(1-damprat(m)^2);
     
217.                     e(m)=exp(-b(m));
     
218.                     a11(m)=e(m)*(damprat(m)*sin(a(m))/d(m)+cos(a(m)));
     
219.                     a12(m)=e(m)*sin(a(m))/angfreq/d(m);
     
220.                     a21(m)=e(m)*(-angfreq)*sin(a(m))/d(m);
     
221.                     a22(m)=e(m)*(-damprat(m)*sin(a(m))/d(m)+cos(a(m)));
     
222.                     b11(m)=e(m)*((2*damprat(m)^2-1+b(m))*sin(a(m))/angfreq^2/a(m)+(2*damprat(m)+angfreq*deltime)*cos(a(m))/angfreq^3/deltime)...
     
223.                            -2*damprat(m)/angfreq^3/deltime;
     
224.                     b12(m)=e(m)*((1-2*damprat(m)^2)*sin(a(m))/angfreq^2/a(m)-2*damprat(m)*cos(a(m))/angfreq^3/deltime)...
     
225.                            -1/angfreq^2+2*damprat(m)/angfreq^3/deltime;
     
226.                     b21(m)=e(m)*((-damprat(m)-angfreq*deltime)*sin(a(m))/angfreq/a(m)-cos(a(m))/angfreq^2/deltime)+1/angfreq^2/deltime;
     
227.                     b22(m)=e(m)*(damprat(m)*sin(a(m))/angfreq/a(m)+cos(a(m))/angfreq^2/deltime)-1/angfreq^2/deltime;
     
228.                     tempdis(m,1)=0;
     
229.                     tempvel(m,1)=0;
     
230.                     tempacc(m,1)=0;
     
231.                      for k=2:signumb
     
232.                          tempdis(m,k)=a11(m)*tempdis(m,k-1)+a12(m)*tempvel(m,k-1)+b11(m)*vibsign(k-1)+b12(m)*vibsign(k);
     
233.                          tempvel(m,k)=a21(m)*tempdis(m,k-1)+a22(m)*tempvel(m,k-1)+b21(m)*vibsign(k-1)+b22(m)*vibsign(k);
     
234.                          tempacc(m,k)=-(2*damprat(m)*angfreq.*tempvel(m,k)+angfreq^2.*tempdis(m,k));
     
235.                      end
     
236.                      dispspe(m,n)=max(abs(tempdis(m,:)));
     
237.                      velospe(m,n)=max(abs(tempvel(m,:)));
     
238.                      accespe(m,n)=max(abs(tempacc(m,:)));
     
239.                      accspeb(m,n)=accespe(m,n)/maxsign;
     
240.                 end
     
241.                 fprintf('共 %4d 个频率计算点，正在计算第  %4d  个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
     
242.             end
     
243.             disp(' ')
     
244.             elatime=cputime-oritime;
     
245.             fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);
     
246.         else
     
247.             for n=1:spenumb
     
248.                 angfreq=2*pi/freoper(n);
     
249.                 deltime=1/sampfre;
     
250.                 for m=1:damnumb
     
251.                     a(m)=angfreq*sqrt(1-damprat(m)^2)*deltime;
     
252.                     b(m)=damprat(m)*angfreq*deltime;
     
253.                     d(m)=sqrt(1-damprat(m)^2);
     
254.                     e(m)=exp(-b(m));
     
255.                     a11(m)=e(m)*(damprat(m)*sin(a(m))/d(m)+cos(a(m)));
     
256.                     a12(m)=e(m)*sin(a(m))/angfreq/d(m);
     
257.                     a21(m)=e(m)*(-angfreq)*sin(a(m))/d(m);
     
258.                     a22(m)=e(m)*(-damprat(m)*sin(a(m))/d(m)+cos(a(m)));
     
259.                     b11(m)=e(m)*((2*damprat(m)^2-1+b(m))*sin(a(m))/angfreq^2/a(m)+(2*damprat(m)+angfreq*deltime)*cos(a(m))/angfreq^3/deltime)...
     
260.                            -2*damprat(m)/angfreq^3/deltime;
     
261.                     b12(m)=e(m)*((1-2*damprat(m)^2)*sin(a(m))/angfreq^2/a(m)-2*damprat(m)*cos(a(m))/angfreq^3/deltime)...
     
262.                            -1/angfreq^2+2*damprat(m)/angfreq^3/deltime;
     
263.                     b21(m)=e(m)*((-damprat(m)-angfreq*deltime)*sin(a(m))/angfreq/a(m)-cos(a(m))/angfreq^2/deltime)+1/angfreq^2/deltime;
     
264.                     b22(m)=e(m)*(damprat(m)*sin(a(m))/angfreq/a(m)+cos(a(m))/angfreq^2/deltime)-1/angfreq^2/deltime;
     
265.                     tempdis(m,1)=0;
     
266.                     tempvel(m,1)=0;
     
267.                     tempacc(m,1)=0;
     
268.                      for k=2:signumb
     
269.                          tempdis(m,k)=a11(m)*tempdis(m,k-1)+a12(m)*tempvel(m,k-1)+b11(m)*vibsign(k-1)+b12(m)*vibsign(k);
     
270.                          tempvel(m,k)=a21(m)*tempdis(m,k-1)+a22(m)*tempvel(m,k-1)+b21(m)*vibsign(k-1)+b22(m)*vibsign(k);
     
271.                          tempacc(m,k)=-(2*damprat(m)*angfreq*tempvel(m,k)+angfreq^2*tempdis(m,k));
     
272.                      end
     
273.                      dispspe(m,n)=max(abs(tempdis(m,:)));
     
274.                      velospe(m,n)=max(abs(tempvel(m,:)));
     
275.                      accespe(m,n)=max(abs(tempacc(m,:)));
     
276.                      accspeb(m,n)=accespe(m,n)/maxsign;
     
277.                 end
     
278.                 fprintf('共 %4d 个周期计算点，正在计算第  %4d  个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
     
279.             end
     
280.             disp(' ')
     
281.             elatime=cputime-oritime;
     
282.             fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);
     
283.         end      
     
284. end
     
285. disp('  ')
     
286. disp('------------------------ ')
     
287. disp('**  数据 & 图形 输出  ** ')
     
288. disp('------------------------ ')
     
289. disp('  ')
     
290. tempdat=date;
     
291. datyorn=input('请选择数据输出选项（0：不输出 1：仅输出位移谱 2：仅输出速度谱 3：仅输出加速度谱(含β谱) 4：输出所有谱）：');
     
292. if datyorn==0
     
293.      disp('--- 没有生成数据文件！！！---')
     
294. elseif datyorn==1
     
295.      fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s');
     
296.      fileid=fopen(fnout,'wt');
     
297.      switch damnumb
     
298.      case 1
     
299.          fprintf(fileid,'* ------------------------------------\n');
     
300.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
301.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
302.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
303.          fprintf(fileid,'* %s\n',tempdat);
     
304.          fprintf(fileid,'* DATA\n');
     
305.          fprintf(fileid,'* ------------------------------------\n');
     
306.          fprintf(fileid,'%8s\t%15s\n','# F/T','# D_Spectrum');
     
307.          fprintf(fileid,'* ------------------------------------\n');
     
308.          for n=1:spenumb
     
309.              fprintf(fileid,'%8.4f\t%15.6E\n',freoper(n),dispspe(1,n));
     
310.          end
     
311.          fprintf(fileid,'*DATAEND ');
     
312.          disp('--- 成功生成数据文件！！！---');
     
313.          filesta=fclose(fileid);
     
314.      case 2
     
315.          fprintf(fileid,'* ------------------------------------\n');
     
316.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
317.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
318.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
319.          fprintf(fileid,'* %s\n',tempdat);
     
320.          fprintf(fileid,'* DATA\n');
     
321.          fprintf(fileid,'* --------------------------------------------------------\n');
     
322.          fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2');
     
323.          fprintf(fileid,'* --------------------------------------------------------\n');
     
324.          for n=1:spenumb
     
325.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n));
     
326.          end
     
327.          fprintf(fileid,'*DATAEND ');
     
328.          disp('--- 成功生成数据文件！！！---');
     
329.          filesta=fclose(fileid);
     
330.      case 3
     
331.          fprintf(fileid,'* ------------------------------------\n');
     
332.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
333.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
334.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
335.          fprintf(fileid,'* %s\n',tempdat);
     
336.          fprintf(fileid,'* DATA\n');
     
337.          fprintf(fileid,'* -------------------------------------------------------------------------------\n');
     
338.          fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2','# D_Spectrum_3');
     
339.          fprintf(fileid,'* -------------------------------------------------------------------------------\n');
     
340.          for n=1:spenumb
     
341.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n),dispspe(3,n));
     
342.          end
     
343.          fprintf(fileid,'*DATAEND ');
     
344.          disp('--- 成功生成数据文件！！！---');
     
345.          filesta=fclose(fileid);
     
346.      end
     
347. elseif datyorn==2
     
348.      fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s');
     
349.      fileid=fopen(fnout,'wt');
     
350.      switch damnumb
     
351.      case 1
     
352.          fprintf(fileid,'* ------------------------------------\n');
     
353.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
354.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
355.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
356.          fprintf(fileid,'* %s\n',tempdat);
     
357.          fprintf(fileid,'* DATA\n');
     
358.          fprintf(fileid,'* ------------------------------------\n');
     
359.          fprintf(fileid,'%8s\t%15s\n','# F/T','# V_Spectrum');
     
360.          fprintf(fileid,'* ------------------------------------\n');
     
361.          for n=1:spenumb
     
362.              fprintf(fileid,'%8.4f\t%15.6E\n',freoper(n),velospe(1,n));
     
363.          end
     
364.          fprintf(fileid,'*DATAEND ');
     
365.          disp('--- 成功生成数据文件！！！---');
     
366.          filesta=fclose(fileid);
     
367.      case 2
     
368.          fprintf(fileid,'* ------------------------------------\n');
     
369.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
370.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
371.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
372.          fprintf(fileid,'* %s\n',tempdat);
     
373.          fprintf(fileid,'* DATA\n');
     
374.          fprintf(fileid,'* ---------------------------------------------------------\n');
     
375.          fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# V_Spectrum_1','# V_Spectrum_2');
     
376.          fprintf(fileid,'* ---------------------------------------------------------\n');
     
377.          for n=1:spenumb
     
378.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),velospe(1,n),velospe(2,n));
     
379.          end
     
380.          fprintf(fileid,'*DATAEND ');
     
381.          disp('--- 成功生成数据文件！！！---');
     
382.          filesta=fclose(fileid);
     
383.      case 3
     
384.          fprintf(fileid,'* ------------------------------------\n');
     
385.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
386.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
387.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
388.          fprintf(fileid,'* %s\n',tempdat);
     
389.          fprintf(fileid,'* DATA\n');
     
390.          fprintf(fileid,'* --------------------------------------------------------------------------------\n');
     
391.          fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# V_Spectrum_1','# V_Spectrum_2','# V_Spectrum_3');
     
392.          fprintf(fileid,'* --------------------------------------------------------------------------------\n');
     
393.          for n=1:spenumb
     
394.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),velospe(1,n),velospe(2,n),velospe(3,n));
     
395.          end
     
396.          fprintf(fileid,'*DATAEND ');
     
397.          disp('--- 成功生成数据文件！！！---');
     
398.          filesta=fclose(fileid);
     
399.      end
     
400. elseif datyorn==3
     
401.      fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s');
     
402.      fileid=fopen(fnout,'wt');
     
403.      switch damnumb
     
404.      case 1
     
405.          fprintf(fileid,'* ------------------------------------\n');
     
406.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
407.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
408.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
409.          fprintf(fileid,'* %s\n',tempdat);
     
410.          fprintf(fileid,'* DATA\n');
     
411.          fprintf(fileid,'* -----------------------------------------\n');
     
412.          fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# A_Spectrum','# A_Spectrum_β');
     
413.          fprintf(fileid,'* -----------------------------------------\n');
     
414.          for n=1:spenumb
     
415.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),accespe(1,n),accspeb(1,n));
     
416.          end
     
417.          fprintf(fileid,'*DATAEND ');
     
418.          disp('--- 成功生成数据文件！！！---');
     
419.          filesta=fclose(fileid);
     
420.      case 2
     
421.          fprintf(fileid,'* ------------------------------------\n');
     
422.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
423.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
424.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
425.          fprintf(fileid,'* %s\n',tempdat);
     
426.          fprintf(fileid,'* DATA\n');
     
427.          fprintf(fileid,'* ------------------------------------------------------------------------------\n');
     
428.          fprintf(fileid,'%8s\t%15s\t%15s\t%15s\t%15s\n','# F/T','# A_Spectrum_1','# A_Spectrum_1β','# A_Spectrum_2','# A_Spectrum_2β');
     
429.          fprintf(fileid,'* ------------------------------------------------------------------------------\n');
     
430.          for n=1:spenumb
     
431.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),accespe(1,n),accspeb(1,n),accespe(2,n),accspeb(2,n));
     
432.          end
     
433.          fprintf(fileid,'*DATAEND ');
     
434.          disp('--- 成功生成数据文件！！！---');
     
435.          filesta=fclose(fileid);
     
436.      case 3
     
437.          fprintf(fileid,'* ------------------------------------\n');
     
438.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
439.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
440.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
441.          fprintf(fileid,'* %s\n',tempdat);
     
442.          fprintf(fileid,'* DATA\n');
     
443.          fprintf(fileid,'* ---------------------------------------------------------------------------------------------------------------------\n');
     
444.          fprintf(fileid,'%8s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\n','# F/T','# A_Spectrum_1','# A_Spectrum_1β','# A_Spectrum_2','# A_Spectrum_2β',...
     
445.                         '# A_Spectrum_3','# A_Spectrum_3β');
     
446.          fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------\n');
     
447.          for n=1:spenumb
     
448.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),accespe(1,n),accspeb(1,n),...
     
449.                      accespe(2,n),accspeb(2,n),accespe(3,n),accspeb(3,n));
     
450.          end
     
451.          fprintf(fileid,'*DATAEND ');
     
452.          disp('--- 成功生成数据文件！！！---');
     
453.          filesta=fclose(fileid);
     
454.       end
     
455.   else
     
456.      fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s');
     
457.      fileid=fopen(fnout,'wt');
     
458.      switch damnumb
     
459.      case 1
     
460.          fprintf(fileid,'* ------------------------------------\n');
     
461.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
462.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
463.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
464.          fprintf(fileid,'* %s\n',tempdat);
     
465.          fprintf(fileid,'* DATA\n');
     
466.          fprintf(fileid,'* --------------------------------------------------------------------------\n');
     
467.          fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum','# V_Spectrum','# A_Spectrum');
     
468.          fprintf(fileid,'* --------------------------------------------------------------------------\n');
     
469.          for n=1:spenumb
     
470.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),velospe(1,n),accespe(1,n));
     
471.          end
     
472.          fprintf(fileid,'*DATAEND ');
     
473.          disp('--- 成功生成数据文件！！！---');
     
474.          filesta=fclose(fileid);
     
475.      case 2
     
476.          fprintf(fileid,'* ------------------------------------\n');
     
477.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
478.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
479.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
480.          fprintf(fileid,'* %s\n',tempdat);
     
481.          fprintf(fileid,'* DATA\n');
     
482.          fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------------------------------------\n');
     
483.          fprintf(fileid,'%8s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2','# V_Spectrum_1','# V_Spectrum_2','# A_Spectrum_1','# A_Spectrum_2');
     
484.          fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------------------------------------\n');
     
485.          for n=1:spenumb
     
486.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n),velospe(1,n),velospe(2,n),accespe(1,n),accespe(2,n));
     
487.          end
     
488.          fprintf(fileid,'*DATAEND ');
     
489.          disp('--- 成功生成数据文件！！！---');
     
490.          filesta=fclose(fileid);
     
491.      case 3
     
492.          fprintf(fileid,'* ------------------------------------\n');
     
493.          fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
     
494.          fprintf(fileid,'* YAHONG DENG (PHD)\n');
     
495.          fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
     
496.          fprintf(fileid,'* %s\n',tempdat);
     
497.          fprintf(fileid,'* DATA\n');
     
498.          fprintf(fileid,'* -------------------------------------------------------------------------------------------------------------------------------------------------------------\n');
     
499.          fprintf(fileid,'%8s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2','# D_Spectrum_3',...
     
500.                         '# V_Spectrum_1','# V_Spectrum_2','# V_Spectrum_3','# A_Spectrum_1','# A_Spectrum_2','# A_Spectrum_3');
     
501.          fprintf(fileid,'* -------------------------------------------------------------------------------------------------------------------------------------------------------------\n');
     
502.          for n=1:spenumb
     
503.              fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n),dispspe(3,n),...
     
504.                             velospe(1,n),velospe(2,n),velospe(3,n),accespe(1,n),accespe(2,n),accespe(3,n));
     
505.          end
     
506.          fprintf(fileid,'*DATAEND ');
     
507.          disp('--- 成功生成数据文件！！！---');
     
508.          filesta=fclose(fileid);
     
509.     end
     
510. end
     
511. disp('')
     
512. figyorn=input('请选择图形显示选项（0：不输出 1：仅输出位移谱 2：仅输出速度谱 3：仅输出加速度谱(含β谱) 4：输出所有谱）：');
     
513. if figyorn==0
     
514.      disp(' ')
     
515.      disp('*******************************')
     
516.      disp('***   程序结束，谢谢使用！  ***')
     
517.      disp('*******************************')
     
518.      disp(' ')
     
519. elseif figyorn==1
     
520.      for m=1:damnumb
     
521.          plot(freoper,dispspe(m,:));
     
522.          grid on;
     
523.          hold on;
     
524.      end
     
525.      xlabel('频率(Hz)/周期(Sec)');
     
526.      ylabel('幅值');
     
527.      title('位移反应谱曲线');
     
528.      disp('--- 成功生成曲线图形！！！---');
     
529.      disp(' ')
     
530.      disp('*******************************')
     
531.      disp('***   程序结束，谢谢使用！  ***')
     
532.      disp('*******************************')
     
533.      disp(' ')
     
534. elseif figyorn==2
     
535.      for m=1:damnumb
     
536.          plot(freoper,velospe(m,:));
     
537.          grid on;
     
538.          hold on;
     
539.      end
     
540.      xlabel('频率(Hz)/周期(Sec)');
     
541.      ylabel('幅值');
     
542.      title('速度反应谱曲线');
     
543.      disp('--- 成功生成曲线图形！！！---');
     
544.      disp(' ')
     
545.      disp('*******************************')
     
546.      disp('***   程序结束，谢谢使用！  ***')
     
547.      disp('*******************************')
     
548.      disp(' ')
     
549. elseif figyorn==3
     
550.      subplot(2,1,1)
     
551.      for m=1:damnumb
     
552.          plot(freoper,accespe(m,:),'r');
     
553.          grid on;
     
554.          hold on;
     
555.      end
     
556.      xlabel('频率(Hz)/周期(Sec)');
     
557.      ylabel('幅值');
     
558.      title('加速度反应谱曲线');
     
559.      legend('幅值谱 ');
     
560.      subplot(2,1,2)
     
561.      for m=1:damnumb
     
562.          plot(freoper,accspeb(m,:));
     
563.          grid on;
     
564.          hold on;
     
565.      end
     
566.      xlabel('频率(Hz)/周期(Sec)');
     
567.      ylabel('动力系数(β)');
     
568.      legend('动力放大系数(β)谱 ');
     
569.      disp('--- 成功生成曲线图形！！！---');
     
570.      disp(' ')
     
571.      disp('*******************************')
     
572.      disp('***   程序结束，谢谢使用！  ***')
     
573.      disp('*******************************')
     
574.      disp(' ')
     
575. else
     
576.      subplot(3,1,1);
     
577.       for m=1:damnumb
     
578.           plot(freoper,dispspe(m,:));
     
579.           grid on;
     
580.           hold on;
     
581.       end
     
582.       xlabel('频率(Hz)/周期(Sec)');
     
583.       ylabel('位移谱值');
     
584.       title('位移/速度/加速度反应谱曲线');
     
585.      subplot(3,1,2)
     
586.       for m=1:damnumb
     
587.           plot(freoper,velospe(m,:));
     
588.           grid on;
     
589.           hold on;
     
590.       end
     
591.       xlabel('频率(Hz)/周期(Sec)');
     
592.       ylabel('速度谱值');
     
593.      subplot(3,1,3)
     
594.       for m=1:damnumb
     
595.           plot(freoper,accespe(m,:),'r');
     
596.           grid on;
     
597.           hold on;
     
598.       end
     
599.       xlabel('频率(Hz)/周期(Sec)');
     
600.       ylabel('加速度谱值');
     
601.       disp('--- 成功生成曲线图形！！！---');
     
602.      disp(' ')
     
603.      disp('*******************************')
     
604.      disp('***   程序结束，谢谢使用！  ***')
     
605.      disp('*******************************')
     
606.      disp(' ')
     
607. end

复制代码

ccb137731

果真高手。

tulong_tulong

个人觉得还是应该把这个转换过程的理论搞清楚，然后自己编个比较通用的基于windows平台的.exe，这样就能省却很多事了。每次只需要给出几个参数就能生成地震波数据

beauty33333

好东东呀

fudy10

有图有真相

manguoyong

很好

xukyle

下来看看那

xukyle

怎么下载分扣了后没有反应呀

xukyle

把northbridge.mat 文件放上看看吧