用MATLAB拟合生成地震波程序
求教各位大大谁有源程序，或者帮忙指导下如何做，有的话发我邮箱里，谢谢了。zhangfanyu88@yahoo.com.cn回复 楼主 zntimmon 的帖子
LZ搜索一下版面：地震波看看此帖是否有用？
http://forum.vibunion.com/thread-25956-1-1.html
[ 本帖最后由 ch_j1985 于 2008-11-20 21:50 编辑 ] clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('生成人工地震波-输入数据文件名(20041012):','s');
fid=fopen(fni,'r');
fs=fscanf(fid,'%f',1);%采样频率
tu=fscanf(fid,'%f',1);
%上升时间包络线线形(1-直线、2-抛物线、3-指数曲线)
iu=fscanf(fid,'%f',1);
cu=fscanf(fid,'%f',1);
ta=fscanf(fid,'%f',1);%持时时间长度
td=fscanf(fid,'%f',1);%下降时间长度
id= fscanf(fid,'%f',1);
cd=fscanf(fid,'%f',1);
dp=fscanf(fid,'%f',1);%阴尼比值
p=fscanf(fid,'%f',1);%概率系数（一般可取P=0.85）
nn=fscanf(fid,'%f',1);%迭代次数
fno=fscanf(fid,'%f',1);%输出数据文件名
x=fscanf(fid,'%f',);%反应谱频率和幅值数据
tatus=fclose(fid);
tl=tu+ta+td;
nt=round(fs*tl+1);
nfft=2^nestpow2(nt)
%计算频率间隔（Hz）
df=fs/nfft
%定义反应谱的离散频率向量
f=0:df:(nfft/2-1)*df
%计算时间间隔（s）
dt=1/fs;
%定义的离散时间向量
t=0:dt:(nt-1)*dt
%生成0到2PI的随机数为随机相位
g=rand(1,nfft/s)*2*pi;
en=ones(1,nt);
l=round(tu*fs)+1
switch iu
case 1 %直线
en(1:l)=linspace(0,1,1);% y = linspace(a,b,n) generates a row vector y of n points linearly
spaced between and including a and b.
case 2 %抛物线
a=0:l-1;
en(1:l)=(a/(l-1)).^2;
case 3 %指数曲线
a=0:l-1;
en(1:l)=1-exp(-cu*a/(l-1));
end
m=round((tu+ta)*fs)+1;
switch id
case 1 %直线
en(m:nt)=linspace(1,0,nt-m+1);
case 2 %抛物线
a=0:nt-m;
en(m:nt)=1-cd*(a/(nt-m)).^2;
case 3 %指数曲线
a=0:nt-m;
en(m:nt)=exp(-cd*a/(nt-m));
end
a0=zeros(x(1,:));
n=length(x(1,:));
nb=round(x(1,n)/df)+1;
for k=1:n-1
l=round(x(1,K)/df)+1;
m=round(x(1,K+1)/df)+1;
a0(1:m)=linspace(x(2,k),x(2,k+1),m-l+1)
end
a1=a0;
s=zeros(1,nfft/2);
k=nb:ne;
s(k)=2*dp/(pi.*(a1(k).^2)./f(k)./(-2*log(-log(p)*pi/tl)./f(k)));
%将功率谱转换成傅里叶幅值谱
b1=sqrt(4*df*s)*nfft/2;
hf=zeros(ne,nfft);
for j-0:ne-1
w=2*pi*df*j
wd=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));
h=wd*e.*(s+c)/fs;
hf(j+1,:)=fft(h,nfft);
end
mm=nn
for k=1:100
c=b1.*exp(i*g);
d=;
e=ifft(d,nfft);
y=en.*real(e(1:nt));
yf=fft(y,nfft);
for j=1:ne
d=ifft(yf.*hf(j,:),nfft);
al(j)=max(real(d(1:nt)));
end
if k==mm
subplot(2,1,1);
plot(t,y,t,en,t,-en);
xlabel('时间（s）');
ylabel('加速度（g）');
grid on
sublpot(2,1,2);
l=1:ne
plot(f(l),a0(l),':'f(l),a1(l));
xlabel('频率（Hz）');
ylabel('加速度（g）');
legend('目标谱','计算谱');
grid on;
ig=input('继续迭代次数[取值1-9,否则退出]:');
if ig>0&ig<10 %如果输入数字是1-9
mm=mm+ig
else
break;
end
end
c=bl
j=nb:ne;
bl(j)=c(j).*a0(j)./a1(j);
end
fid=fopen(fno,'W');
for K=1:nt
fprintf(fid,'%f%f\n',t(k),y(k));
end
status=fclose(fid);
页:
[1]