地震记录合成

clear all;
clc;
L=71;
fm=40;%采样频率
dt=0.001;%时间间隔
R=zeros（1L）;%雷克子波长度
%雷克子波
B=（L-1）/2;
for m=1:1:L    
  
    t=dt*（m-B）;
    R（m）=（1-2*（t*pi*fm）^2）*exp（-（t*pi*fm）^2）;
end


subplot（241）;
 plot（R）;hold on;

% %卷积
% RE=zeros（1128）;
% RE（164）=0.5;%反射系数
% Result=zeros（128+L-1）;%卷积长度
% Result=conv（RRE）;
%  subplot（242）;
% plot（Result）;hold on;
% title（‘一维卷积‘）

%时深转换
c=256;
M=load（‘Model_V.txt‘）;%输入数据
V=（reshape（Mcc））‘;%转换成256*256的矩阵
dh=10;%深度间隔

t=zeros（2c）;  %时间域分界面矩阵      
T=zeros（cc）;%时间矩阵

   
        
for j=1:1:c  %按列循环
    T（1j）=0;  %第一行时间为零
    d=1;    %深度分界面矩阵行数
    for i=2:1:c
        T（ij）=2*dh/V（ij）+T（i-1j）;%时间矩阵
         if V（i-1j）          h（dj）=（i-1）*dh;%深度分界面矩阵
          d=d+1;
         end
         t（dj）=ceil（T（ij）/0.001）;
           
         
    end
    
        
end
%时间域赋速度值
for j=1:1:c
    
    n=0;
    for i=2:1:c
        for k=n+1:1:ceil（T（ij）/0.001）
            
            Vt（kj）=V（i-1j）;
            n=ceil（T（ij）/0.001）;
        
               
        end
     end
    
end
I=size（Vt1）;%取行数
for j=1:1:c
    for k=1:1:I
         if k>t（1j） +t（2j）
             Vt（kj）=V（cc）;
         end   
    end
end
%反射系数

            
for j=1:1:c  %按列循环
   
    for i=2:1:I
        
         Re（ij）=（Vt（ij）-Vt（