大家好！下面是我在做复数道分析时求瞬时频率的程序，可是结果总是不正确，哪位能帮下忙，帮斧正一下，并和我联系,将不胜感激！
clc;
clear;
depth1=100; %第一层厚度 
depth2=[80];%第二层厚度
velocity_up=1500;
velocity_mid=2000;
velocity_down=2500;
sample_inter=1;%采样间隔1ms
trace_num=1;%道数
trace_length=600;%整个地层深度
wave_t=160; %define the wavelet length 160ms
fp=25; %the wavelet main frequence is 25hz
t1=2*depth1/velocity_up;%深时转换
t2=t1+(2*depth2/velocity_up);
ref1= (velocity_mid-velocity_up)/(velocity_mid+velocity_up);% ref1=0.1429
ref2= (velocity_down-velocity_mid)/(velocity_down+velocity_mid); % ref2=0.1111

sample_num=trace_length/sample_inter;%地层深度采样点个数
for i=1:trace_num
for j=1:sample_num
R(j,i)=0;
end
k1=round(t1*1000/sample_inter); % k1=133
R(k1,i)=ref1;
k2=round(t2*1000/ sample_inter); % K2=[240]
R(k2(i),i)=ref2; 
end
w_t=-wave_t/(2*1000):0.001:wave_t/(2*1000);
rick=(1-2*(pi*fp*w_t).^2).*exp(-(pi*fp*w_t).^2);
for i=1:trace_num
syn(:,i)=conv(R(:,i),rick);
figure(1);
plot(syn(:,i));
axis([0 500 -0.3 0.4]);
z=hilbert(syn(:,i));
a=real(z);
b=imag(z);
c=atan(b/a);
%c=angle(z);
%c= unwrap(c);
figure(2);
plot(c);
axis([0 500 -1 1]);
P=diff(c);
figure(3);
plot(P);
axis([0 500 -0.3 0.4]);
end

我的邮箱是wolfshouse·163.com。