super_cycle.m

我写的一个地震单台定位的程序

function y = super_cycle();
%use recursion to get the superiority cycle of the P wave
%close all
global data;
if nargin==0
    start = 10;
    t = 10;
else if nargin == 1
        t = 10;
    else if nargin>2
            error('this function only need two parameters.');
        end
    end
end
H = get(gcf,'UserData');
start = str2num(get(H(1),'String'));
t = str2num(get(H(2),'String'));
fs = str2num(get(H(3),'String'));
alfa = str2num(get(H(4),'String'));
%alfa = 0.225 ;%alfa is a moving coefficient
%fs = 100;
t = t*fs;
start = start*fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
[name path] = uigetfile('*.*','please choose a file','please open a file');
file = strcat(path,name);
data = load(file);
%}
up = diff(data(:,1));
fft_up = abs(fft(up));
%{
figure
x = 1:1:length(fft_up);
x = x./length(fft_up);
x = x*100;
subplot(311)
plot(x,fft_up)
ew = data(:,2);
ns = data(:,3);
%}
Ap = max(abs(up(start:start+t)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ml = length(up);
N = 5;
for i = 5:ml
    R_up(i) = sum(up(i-N+1:1:i))/N;
end
%figure
%plot(up);
d_up = diff(up);
D = diff(R_up);
R_up = R_up.*R_up;
D = D.*D;
for k = start:1:t+start
    R_up(k) = alfa*R_up(k-1) + up(k)^2;
    D(k) = alfa*D(k-1) + d_up(k)^2;
    T(k) = 2*pi*sqrt(R_up(k)/D(k))/100;
end
%{
figure
plot(T(start:start+t))
%}
y = max(T);
%{
disp('superiority cycle: ');
disp(y);
%}
Ml = 6.3*log10(y)+7.1;
set(H(6),'String',Ml);
%{
disp('the grade of the event: ');
disp(Ml);
%}
epicenter = 10^(-0.51118*log10(1/y)-0.18298*log10(Ap)+1.59766);
%{
disp('the distance of the epicenter is: ');
disp(delta);
%}
set(H(5),'String',epicenter);