zxwflyer 发表于 2007-4-16 16:57

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

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

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

luwilero 发表于 2007-4-17 10:38

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

jeruselem 发表于 2008-4-28 19:44

地震波反应谱

本帖最后由 牛小贱 于 2015-3-18 20:54 编辑

clear;close all;
load northbridge.mat %输入地震波文件，时程数据格式为一行。
eqw=x;
amax=max(abs(eqw));
dt=0.01%采样频率100Hz
ff=eqw;
n=0;
zeta=0.05; %阻尼比
for w=200:-0.1:1; %自己设定
x=0;
v=0;
%alhp=35;
acc=0;
a=;
b=1;
=residue(b,a);
t=0:dt:length(eqw)*dt;
h=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t);
x=conv(h,ff)*dt;
v=diff(x)/dt;
acc=diff(v)/dt;
n=n+1;
beta1(n)=max(abs(acc))/amax;
beta2(n)=max(abs(v))/amax;
beta3(n)=max(abs(x))/amax;
T(n)=2*pi/w;
end;
save eqw_response_Northbridge T beta1 beta2 beta3
%%
ww=1./T;wf=fliplr(ww);
beta1f=fliplr(beta1);
beta2f=fliplr(beta2);
beta3f=fliplr(beta3);
figure;plot(wf,beta1f,'linewidth',2);grid on;xlabel('frequency (Hz) ');ylabel('acceleration (m/s^2)');title('Acceleration Response Spectrum-Northbridge 35gal ');
saveas(gcf,'Northbridge-Response-1','fig');saveas(gcf,'Northbridge-Response-1','bmp');
figure;plot(wf,beta2f,'linewidth',2);grid on;xlabel('frequency (Hz)');ylabel('velocity (m/s)');title('Velocity Response Spectrum-Northbridge 35gal')
saveas(gcf,'Northbridge-Response-2','fig');saveas(gcf,'Northbridge-Response-2','bmp');
figure;plot(wf,beta3f,'linewidth',2);grid on;xlabel('frequency (Hz)');ylabel('displacement (m)');title('Displacement Response Spectrum-Northbridge 35gal')
saveas(gcf,'Northbridge-Response-3','fig');saveas(gcf,'Northbridge-Response-3','bmp');
figure;plot(T,beta1,'linewidth',2);grid on;xlabel('period (s)');ylabel('acceleration (m/s^2)');title('Acceleration Response Spectrum-Northbridge 35gal ');
saveas(gcf,'Northbridge-Response-4','fig');saveas(gcf,'Northbridge-Response-4','bmp');
figure;plot(T,beta2,'linewidth',2);grid on;xlabel('period (s)');ylabel('velocity (m/s)');title('Velocity Response Spectrum-Northbridge 35gal ');
saveas(gcf,'Northbridge-Response-5','fig');saveas(gcf,'Northbridge-Response-5','bmp');
figure;plot(T,beta3,'linewidth',2);grid on;xlabel('period (s)');ylabel('displacement (m)');title('Displacement Response Spectrum-Northbridge 35gal ');
saveas(gcf,'Northbridge-Response-6','fig');saveas(gcf,'Northbridge-Response-6','bmp');

liujiashun 发表于 2010-3-4 12:26

回复 板凳 jeruselem 的帖子

谢谢啊 好东西

liujiashun 发表于 2010-4-3 20:56

谢谢

谢谢啊 呵呵 真是高手！

lihong0880 发表于 2010-12-8 10:51

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

cuidazhe 发表于 2010-12-17 14:14

%   ***********************************************
%   *** Response Spectrum By Duhamel Integral & ***
%   ***   Newmark Linear Acceleration Method    ***
%   ***         Author:Yahong DENG (PHD)      ***
%   ***      Email:hoverdyh@zju.edu.cn      ***
%   ***         Chang'an University         ***
%   ***            2008.12.11Xi'an            ***
%   ***********************************************
clc;
clear all;
close all hidden;
disp('***********************************************')
disp('*** Response Spectrum By Duhamel Integral & ***')
disp('***       Linear Acceleration Method      ***')
disp('***      Author:Yahong DENG (PHD)         ***')
disp('***       Email:hoverdyh@zju.edu.cn         ***')
disp('***          Chang An University            ***')
disp('***         2008.12.11XiAn            ***')
disp('***********************************************')
disp('')
disp('-------------------- ')
disp('**振动信号输入** ')
disp('-------------------- ')
disp('')
filenin=input('请输入数据文件名(*.txt/*.dat/*.*)：','s');
fileid=fopen(filenin,'r');
sampfre=input('请输入数据采样频率：');
linsign=input('请输入数据起始行：');
for line=1:linsign-1
    tline=fgetl(fileid);
end
filesty=input('请选择数据列类型（1：全部(单列/多列) 2：多列中一列）：');
if filesty==1
    vibsign=fscanf(fileid,'%f',inf);
    vibsign=vibsign';
else
    colnumb=input('请输入文件总列数：');
    colsign=input('请输入数据所在列数：');
    inpsign=fscanf(fileid,'%f',);
    vibsign=inpsign(colsign,:);
end
signumb=length(vibsign);
time=(0:1/sampfre:(signumb-1)/sampfre);
maxsign=0;
for n=1:signumb
    if abs(vibsign(n))>maxsign
      maxsign=abs(vibsign(n));
      maxtime=time(n);
    end
end
disp(' ');
fprintf('---* 输入信号基本信息 *---\n');
fprintf('输入信号长度为：%d\n',signumb);
fprintf('输入信号持时为：%.3f (s)\n',(signumb-1)/sampfre);
fprintf('输入信号最大值（绝对值）：%f\n',maxsign);
fprintf('输入信号最大值时间点：%.3f (s)\n',maxtime);
fprintf('--- 恭喜，成功读入数据！！！ ---\n');
filesta=fclose(fileid);
disp('')
disp('---------------------- ')
disp('**反应谱参数设置** ')
disp('---------------------- ')
disp('')
damnumb=input('请输入反应谱阻尼比个数（1-3）：');
switch damnumb
    case 1
      damprat(1)=input('请输入阻尼比值（0-1）：');
    case 2
      damprat(1)=input('请输入第一个阻尼比值（0-1）：');
      damprat(2)=input('请输入第二个阻尼比值（0-1）：');
    case 3
      damprat(1)=input('请输入第一个阻尼比值（0-1）：');
      damprat(2)=input('请输入第二个阻尼比值（0-1）：');
      damprat(3)=input('请输入第三个阻尼比值（0-1）：');
end
haxisty=input('请选择反应谱计算点类型（1：等间距频率计算点 2：等间距周期计算点）：');
switch haxisty
    case 1
      minfreq=input('请输入最小计算频率（一般可取0.1，且 > 0）：');
      maxfreq=input('请输入最大计算频率（一般可取10-15）：');
      delfreq=input('请输入计算频率间隔（一般可取0.01/0.02）：');
      spenumb=round((maxfreq-minfreq)/delfreq)+1;
      freoper=(minfreq:delfreq:(spenumb-1)*delfreq+minfreq);
    case 2
      minperi=input('请输入最小计算周期（一般可取0.1，且 > 0）：');
      maxperi=input('请输入最大计算周期（一般可取2-6）：');
      delperi=input('请输入计算周期间隔（一般可取0.01/0.02）：');
      spenumb=round((maxperi-minperi)/delperi)+1;
      freoper=(minperi:delperi:(spenumb-1)*delperi+minperi);
end
disp('')
disp('-------------------- ')
disp('**   反应谱计算   ** ')
disp('-------------------- ')
disp('')
compmet=input('请选择反应谱计算方法（1：杜哈姆积分卷积法(速度慢) 2：杜哈姆积分FFT法 3：纽马克线性加速度法(速度快,推荐)）：');
oritime=cputime;
switch compmet
    case 1
      if haxisty==1
            for n=1:spenumb
                angfreq=2*pi*freoper(n);
                for m=1:damnumb
                  angfred(m)=angfreq*sqrt(1-damprat(m)^2);
                  termexp=exp(-damprat(m)*angfreq*time);
                  termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2));
                  termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2)));
                  uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre;
                  uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre;
                  uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre;
                  dispspe(m,n)=max(abs(conv(vibsign,uniresd)));
                  velospe(m,n)=max(abs(conv(vibsign,uniresv)));
                  accespe(m,n)=max(abs(conv(vibsign,uniresa)));
                  accspeb(m,n)=accespe(m,n)/maxsign;
                  fftnumb=2^nextpow2(2*signumb);
                  uniafft=fft(uniresa,fftnumb);
                  signfft=fft(vibsign,fftnumb);
                  accespe1=real(ifft(uniafft.*signfft,fftnumb));
                  accespe2=accespe1(1:signumb);
                  accespe3=max(accespe2);
                end
                fprintf('共 %4d 个频率计算点，正在计算第%4d个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
            end
            disp(' ')
            elatime=cputime-oritime;
            fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);   
      else
            for n=1:spenumb
                angfreq=2*pi/freoper(n);
                for m=1:damnumb
                  angfred(m)=angfreq*sqrt(1-damprat(m)^2);
                  termexp=exp(-damprat(m)*angfreq*time);
                  termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2));
                  termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2)));
                  uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre;
                  uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre;
                  uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre;
                  dispspe(m,n)=max(abs(conv(vibsign,uniresd)));
                  velospe(m,n)=max(abs(conv(vibsign,uniresv)));
                  accespe(m,n)=max(abs(conv(vibsign,uniresa)));
                  accspeb(m,n)=accespe(m,n)/maxsign;
                end   
            fprintf('共 %4d 个周期计算点，正在计算第   %4d个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
            end
            disp(' ')
            elatime=cputime-oritime;
            fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);
      end
    case 2
      fftnumb=2^nextpow2(2*signumb-1);
      if haxisty==1
            for n=1:spenumb
                angfreq=2*pi*freoper(n);
                for m=1:damnumb
                  angfred(m)=angfreq*sqrt(1-damprat(m)^2);
                  termexp=exp(-damprat(m)*angfreq*time);
                  termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2));
                  termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2)));
                  uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre;
                  uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre;
                  uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre;
                  fftunid=fft(uniresd,fftnumb);
                  fftuniv=fft(uniresv,fftnumb);
                  fftunia=fft(uniresa,fftnumb);
                  fftsign=fft(vibsign,fftnumb);
                  dispspt=real(ifft(fftunid.*fftsign,fftnumb));
                  accespe(m,n)=max(dispspt(1:signumb));
                  velospt=real(ifft(fftuniv.*fftsign,fftnumb));
                  velospe(m,n)=max(velospt(1:signumb));
                  accespt=real(ifft(fftunia.*fftsign,fftnumb));
                  accespe(m,n)=max(accespt(1:signumb));
                  accspeb(m,n)=accespe(m,n)/maxsign;
                end
                fprintf('共 %4d 个频率计算点，正在计算第%4d个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
            end
            disp(' ')
            elatime=cputime-oritime;
            fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);   
      else
            for n=1:spenumb
                angfreq=2*pi/freoper(n);
                for m=1:damnumb
                  angfred(m)=angfreq*sqrt(1-damprat(m)^2);
                  termexp=exp(-damprat(m)*angfreq*time);
                  termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2));
                  termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2)));
                  uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre;
                  uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre;
                  uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre;
                  fftunid=fft(uniresd,fftnumb);
                  fftuniv=fft(uniresv,fftnumb);
                  fftunia=fft(uniresa,fftnumb);
                  fftsign=fft(vibsign,fftnumb);
                  dispspt=real(ifft(fftunid.*fftsign,fftnumb));
                  accespe(m,n)=max(dispspt(1:signumb));
                  velospt=real(ifft(fftuniv.*fftsign,fftnumb));
                  velospe(m,n)=max(velospt(1:signumb));
                  accespt=real(ifft(fftunia.*fftsign,fftnumb));
                  accespe(m,n)=max(accespt(1:signumb));
                  accspeb(m,n)=accespe(m,n)/maxsign;
                end   
            fprintf('共 %4d 个周期计算点，正在计算第   %4d个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
            end
            disp(' ')
            elatime=cputime-oritime;
            fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);
      end
    case 3
      if haxisty==1
            for n=1:spenumb
                angfreq=2*pi*freoper(n);
                deltime=1/sampfre;
                for m=1:damnumb
                  a(m)=angfreq*sqrt(1-damprat(m)^2)*deltime;
                  b(m)=damprat(m)*angfreq*deltime;
                  d(m)=sqrt(1-damprat(m)^2);
                  e(m)=exp(-b(m));
                  a11(m)=e(m)*(damprat(m)*sin(a(m))/d(m)+cos(a(m)));
                  a12(m)=e(m)*sin(a(m))/angfreq/d(m);
                  a21(m)=e(m)*(-angfreq)*sin(a(m))/d(m);
                  a22(m)=e(m)*(-damprat(m)*sin(a(m))/d(m)+cos(a(m)));
                  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)...
                           -2*damprat(m)/angfreq^3/deltime;
                  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)...
                           -1/angfreq^2+2*damprat(m)/angfreq^3/deltime;
                  b21(m)=e(m)*((-damprat(m)-angfreq*deltime)*sin(a(m))/angfreq/a(m)-cos(a(m))/angfreq^2/deltime)+1/angfreq^2/deltime;
                  b22(m)=e(m)*(damprat(m)*sin(a(m))/angfreq/a(m)+cos(a(m))/angfreq^2/deltime)-1/angfreq^2/deltime;
                  tempdis(m,1)=0;
                  tempvel(m,1)=0;
                  tempacc(m,1)=0;
                     for k=2:signumb
                         tempdis(m,k)=a11(m)*tempdis(m,k-1)+a12(m)*tempvel(m,k-1)+b11(m)*vibsign(k-1)+b12(m)*vibsign(k);
                         tempvel(m,k)=a21(m)*tempdis(m,k-1)+a22(m)*tempvel(m,k-1)+b21(m)*vibsign(k-1)+b22(m)*vibsign(k);
                         tempacc(m,k)=-(2*damprat(m)*angfreq.*tempvel(m,k)+angfreq^2.*tempdis(m,k));
                     end
                     dispspe(m,n)=max(abs(tempdis(m,:)));
                     velospe(m,n)=max(abs(tempvel(m,:)));
                     accespe(m,n)=max(abs(tempacc(m,:)));
                     accspeb(m,n)=accespe(m,n)/maxsign;
                end
                fprintf('共 %4d 个频率计算点，正在计算第%4d个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
            end
            disp(' ')
            elatime=cputime-oritime;
            fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);
      else
            for n=1:spenumb
                angfreq=2*pi/freoper(n);
                deltime=1/sampfre;
                for m=1:damnumb
                  a(m)=angfreq*sqrt(1-damprat(m)^2)*deltime;
                  b(m)=damprat(m)*angfreq*deltime;
                  d(m)=sqrt(1-damprat(m)^2);
                  e(m)=exp(-b(m));
                  a11(m)=e(m)*(damprat(m)*sin(a(m))/d(m)+cos(a(m)));
                  a12(m)=e(m)*sin(a(m))/angfreq/d(m);
                  a21(m)=e(m)*(-angfreq)*sin(a(m))/d(m);
                  a22(m)=e(m)*(-damprat(m)*sin(a(m))/d(m)+cos(a(m)));
                  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)...
                           -2*damprat(m)/angfreq^3/deltime;
                  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)...
                           -1/angfreq^2+2*damprat(m)/angfreq^3/deltime;
                  b21(m)=e(m)*((-damprat(m)-angfreq*deltime)*sin(a(m))/angfreq/a(m)-cos(a(m))/angfreq^2/deltime)+1/angfreq^2/deltime;
                  b22(m)=e(m)*(damprat(m)*sin(a(m))/angfreq/a(m)+cos(a(m))/angfreq^2/deltime)-1/angfreq^2/deltime;
                  tempdis(m,1)=0;
                  tempvel(m,1)=0;
                  tempacc(m,1)=0;
                     for k=2:signumb
                         tempdis(m,k)=a11(m)*tempdis(m,k-1)+a12(m)*tempvel(m,k-1)+b11(m)*vibsign(k-1)+b12(m)*vibsign(k);
                         tempvel(m,k)=a21(m)*tempdis(m,k-1)+a22(m)*tempvel(m,k-1)+b21(m)*vibsign(k-1)+b22(m)*vibsign(k);
                         tempacc(m,k)=-(2*damprat(m)*angfreq*tempvel(m,k)+angfreq^2*tempdis(m,k));
                     end
                     dispspe(m,n)=max(abs(tempdis(m,:)));
                     velospe(m,n)=max(abs(tempvel(m,:)));
                     accespe(m,n)=max(abs(tempacc(m,:)));
                     accspeb(m,n)=accespe(m,n)/maxsign;
                end
                fprintf('共 %4d 个周期计算点，正在计算第%4d个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100);
            end
            disp(' ')
            elatime=cputime-oritime;
            fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime);
      end      
end
disp('')
disp('------------------------ ')
disp('**数据 & 图形 输出** ')
disp('------------------------ ')
disp('')
tempdat=date;
datyorn=input('请选择数据输出选项（0：不输出 1：仅输出位移谱 2：仅输出速度谱 3：仅输出加速度谱(含β谱) 4：输出所有谱）：');
if datyorn==0
   disp('--- 没有生成数据文件！！！---')
elseif datyorn==1
   fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s');
   fileid=fopen(fnout,'wt');
   switch damnumb
   case 1
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'%8s\t%15s\n','# F/T','# D_Spectrum');
         fprintf(fileid,'* ------------------------------------\n');
         for n=1:spenumb
             fprintf(fileid,'%8.4f\t%15.6E\n',freoper(n),dispspe(1,n));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   case 2
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* --------------------------------------------------------\n');
         fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2');
         fprintf(fileid,'* --------------------------------------------------------\n');
         for n=1:spenumb
             fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   case 3
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* -------------------------------------------------------------------------------\n');
         fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2','# D_Spectrum_3');
         fprintf(fileid,'* -------------------------------------------------------------------------------\n');
         for n=1:spenumb
             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));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   end
elseif datyorn==2
   fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s');
   fileid=fopen(fnout,'wt');
   switch damnumb
   case 1
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'%8s\t%15s\n','# F/T','# V_Spectrum');
         fprintf(fileid,'* ------------------------------------\n');
         for n=1:spenumb
             fprintf(fileid,'%8.4f\t%15.6E\n',freoper(n),velospe(1,n));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   case 2
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* ---------------------------------------------------------\n');
         fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# V_Spectrum_1','# V_Spectrum_2');
         fprintf(fileid,'* ---------------------------------------------------------\n');
         for n=1:spenumb
             fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),velospe(1,n),velospe(2,n));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   case 3
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* --------------------------------------------------------------------------------\n');
         fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# V_Spectrum_1','# V_Spectrum_2','# V_Spectrum_3');
         fprintf(fileid,'* --------------------------------------------------------------------------------\n');
         for n=1:spenumb
             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));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   end
elseif datyorn==3
   fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s');
   fileid=fopen(fnout,'wt');
   switch damnumb
   case 1
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* -----------------------------------------\n');
         fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# A_Spectrum','# A_Spectrum_β');
         fprintf(fileid,'* -----------------------------------------\n');
         for n=1:spenumb
             fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),accespe(1,n),accspeb(1,n));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   case 2
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* ------------------------------------------------------------------------------\n');
         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β');
         fprintf(fileid,'* ------------------------------------------------------------------------------\n');
         for n=1:spenumb
             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));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   case 3
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* ---------------------------------------------------------------------------------------------------------------------\n');
         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β',...
                        '# A_Spectrum_3','# A_Spectrum_3β');
         fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------\n');
         for n=1:spenumb
             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),...
                     accespe(2,n),accspeb(2,n),accespe(3,n),accspeb(3,n));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
      end
else
   fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s');
   fileid=fopen(fnout,'wt');
   switch damnumb
   case 1
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* --------------------------------------------------------------------------\n');
         fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum','# V_Spectrum','# A_Spectrum');
         fprintf(fileid,'* --------------------------------------------------------------------------\n');
         for n=1:spenumb
             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));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   case 2
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------------------------------------\n');
         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');
         fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------------------------------------\n');
         for n=1:spenumb
             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));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
   case 3
         fprintf(fileid,'* ------------------------------------\n');
         fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n');
         fprintf(fileid,'* YAHONG DENG (PHD)\n');
         fprintf(fileid,'* CHANG_AN UNIVERSITY\n');
         fprintf(fileid,'* %s\n',tempdat);
         fprintf(fileid,'* DATA\n');
         fprintf(fileid,'* -------------------------------------------------------------------------------------------------------------------------------------------------------------\n');
         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',...
                        '# V_Spectrum_1','# V_Spectrum_2','# V_Spectrum_3','# A_Spectrum_1','# A_Spectrum_2','# A_Spectrum_3');
         fprintf(fileid,'* -------------------------------------------------------------------------------------------------------------------------------------------------------------\n');
         for n=1:spenumb
             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),...
                            velospe(1,n),velospe(2,n),velospe(3,n),accespe(1,n),accespe(2,n),accespe(3,n));
         end
         fprintf(fileid,'*DATAEND ');
         disp('--- 成功生成数据文件！！！---');
         filesta=fclose(fileid);
    end
end
disp('')
figyorn=input('请选择图形显示选项（0：不输出 1：仅输出位移谱 2：仅输出速度谱 3：仅输出加速度谱(含β谱) 4：输出所有谱）：');
if figyorn==0
   disp(' ')
   disp('*******************************')
   disp('***   程序结束，谢谢使用！***')
   disp('*******************************')
   disp(' ')
elseif figyorn==1
   for m=1:damnumb
         plot(freoper,dispspe(m,:));
         grid on;
         hold on;
   end
   xlabel('频率(Hz)/周期(Sec)');
   ylabel('幅值');
   title('位移反应谱曲线');
   disp('--- 成功生成曲线图形！！！---');
   disp(' ')
   disp('*******************************')
   disp('***   程序结束，谢谢使用！***')
   disp('*******************************')
   disp(' ')
elseif figyorn==2
   for m=1:damnumb
         plot(freoper,velospe(m,:));
         grid on;
         hold on;
   end
   xlabel('频率(Hz)/周期(Sec)');
   ylabel('幅值');
   title('速度反应谱曲线');
   disp('--- 成功生成曲线图形！！！---');
   disp(' ')
   disp('*******************************')
   disp('***   程序结束，谢谢使用！***')
   disp('*******************************')
   disp(' ')
elseif figyorn==3
   subplot(2,1,1)
   for m=1:damnumb
         plot(freoper,accespe(m,:),'r');
         grid on;
         hold on;
   end
   xlabel('频率(Hz)/周期(Sec)');
   ylabel('幅值');
   title('加速度反应谱曲线');
   legend('幅值谱 ');
   subplot(2,1,2)
   for m=1:damnumb
         plot(freoper,accspeb(m,:));
         grid on;
         hold on;
   end
   xlabel('频率(Hz)/周期(Sec)');
   ylabel('动力系数(β)');
   legend('动力放大系数(β)谱 ');
   disp('--- 成功生成曲线图形！！！---');
   disp(' ')
   disp('*******************************')
   disp('***   程序结束，谢谢使用！***')
   disp('*******************************')
   disp(' ')
else
   subplot(3,1,1);
      for m=1:damnumb
          plot(freoper,dispspe(m,:));
          grid on;
          hold on;
      end
      xlabel('频率(Hz)/周期(Sec)');
      ylabel('位移谱值');
      title('位移/速度/加速度反应谱曲线');
   subplot(3,1,2)
      for m=1:damnumb
          plot(freoper,velospe(m,:));
          grid on;
          hold on;
      end
      xlabel('频率(Hz)/周期(Sec)');
      ylabel('速度谱值');
   subplot(3,1,3)
      for m=1:damnumb
          plot(freoper,accespe(m,:),'r');
          grid on;
          hold on;
      end
      xlabel('频率(Hz)/周期(Sec)');
      ylabel('加速度谱值');
      disp('--- 成功生成曲线图形！！！---');
   disp(' ')
   disp('*******************************')
   disp('***   程序结束，谢谢使用！***')
   disp('*******************************')
   disp(' ')
end

ccb137731 发表于 2010-12-17 14:24

果真高手。

tulong_tulong 发表于 2011-3-29 11:43

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

beauty33333 发表于 2011-11-17 16:35

好东东呀

fudy10 发表于 2011-11-19 13:01

有图有真相

manguoyong 发表于 2011-11-20 22:32

很好

xukyle 发表于 2011-12-20 10:00

下来看看那

xukyle 发表于 2011-12-20 10:05

怎么下载分扣了后没有反应呀{:3_50:}

xukyle 发表于 2011-12-20 10:16

把northbridge.mat 文件放上看看吧
页: [1] 2
查看完整版本: 求把地震波时程转换成反应谱的matlab 程序