提取瑞雷面波的频散曲线

提取瑞雷面波的频散曲线，并用频散曲线模拟面波并成像，随即进行面波的逆频散

clear all;clc;
x=[02.557.51012.51517.52022.523.52526.527.528.5303132.533353637.538404142.5434547.5485052.55557.56062.56567.570];
y=[200199.5199198.5198197196.5196195194193192189186183177173.5167164154149142.5139130126121119.5113.5107.5107104101.51009998.59897.59797];
%  plot（xy）;
%ylim（[80 220]）;
% hold on;
y1=[203196.5194198.5199195196.5194193196195192187186183176173.5167163154149145.5139130126121119.5113.5107.5107104101.51009998.59897.59797];
scatter（xy1‘r‘）;
xlabel（‘f（HZ）‘）;
ylabel（‘Phase Velocity（m/s）‘）;

dt=0.001;
nfft=1024;
df=1.0/（dt*nfft）;
xx=df:df:70;
yy1=spline（xyxx）;
figure
plot（xxyy1）;

yy2=yy1（length（yy1））*ones（1nfft-length（yy1））;

yy=[yy1 yy2];
% %------------------------------------定义子波
fp=35;
dt=0.002;
t=-0.063:dt:0.064;
%wave1=（1-2*（pi*fp.*t）.^2）.*exp（-（pi*fp.*t）.^2）; 
%figure;plot（twave1）;title（‘雷克子波‘）; 

at=sin（2*pi*fp.*t）./（pi.*t）;
wave1=at;
w=wave1;
w=w./max（abs（w））;
lw=length（w）;
%---
figure;plot（tw）;
af=abs（fft（wnfft））;
df=1/（dt*nfft）;
for j=1:nfft
    f1（j）=j*df;
end

af（1）=0;af（nfft/2+1）=0;%将不对称点设为零值。
for kk=1:35
for k=1:nfft
    if k<=nfft/2
        wf（k）=af（k）*exp（-i*f1（k）*2*pi*kk/yy（k））;%*exp（-i*f1（k）*2*pi*kk/yy（k））;
    else
        wf（k）=af（k）*exp（i*f1（k）*2*pi*kk/yy（k））;%*exp（i*f1（k）*2*pi*kk/yy（k））;
    end
end
 id=real（ifft（wfnfft））;
 dd（:kk）=id（1:nfft/2）‘;
end
figure;wigb（dd0.8）;
x=1;n=1:1:35;
h=（n-1）*x;
dt=0.002;  
v1 =200;
flow=97;
fhigh=200;
[d]=my_inverse_dispertion（dddthflowfhighyyv1）
%d=d（1:200:）;
figure;wigb（d0.8）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        1695  2013-03-09 19:59  mianbo.m