帮我看看这个算法的matlab描述

wuqiong

这里有公式的内容

[ 本帖最后由 wuqiong 于 2008-12-29 17:49 编辑 ]

wuqiong

这里是一个根据三份量地震波记录求解方位角的算法：
震源方位角是运用垂直及水平分量记录的互相关函数(cross-correlations)来求取的（Nakamura ，1988）。垂直及水平分量的互相关函数，可以由以下公式所求得：




其中
为第i点垂直与东西分量的互相关函数
为第i-1点垂直与东西分量的互相关函数
为第i点垂直与南北分量的互相关函数
为第i-1点垂直与南北分量的互相关函数
为第i点垂直分量记录
为第i点东西分量记录
为第i点南北分量记录
为平滑常数(0.0-1.0)
由公式(1)及(2)求得互相关函数之后，再带入下式可求得地震波入射之方位角 。

---(3)
方位角 之所在象限和 与 之正负符号有关，其对应关系如下：

象限
I	负	负
II	正	负
III	正	正
IV	负	正

我的matlab函数为：
function az = azimath_naka(record,alfa,type)
%this function is used to compute the azimth of the seismic record by using
%the method of NAKAMURA. and the alfa is a moving average parameter range
%from 0.0 to 1.0. the parameter type is the type of the seismeter:
%                     1 ------------------- vocerty
%                     2 ------------------- displacement
%                     3 ------------------- aceleration
if nargin ~= 3
    error('the parameter type and alfa must be confirmed at first');
end
if type == 3
    record(:,1) = cumtrapz(record(:,1));
    record(:,2) = cumtrapz(record(:,2));
    record(:,3) = cumtrapz(record(:,3));
end
ud = detrend(record(:,1));
ew = detrend(record(:,2));
ns = detrend(record(:,3));
% figure
% subplot(311),plot(ud);ylabel('ud')
% subplot(312),plot(ns);ylabel('ns')
% subplot(313),plot(ew);ylabel('ew')
%compute the crosscorr of the ud-ew and ud-ns
in = length(ud);
maxlag = in;samp_seg = 100;overlap = 0;flag = 'biased';
% Rue  = cum2x(ud,ew,maxlag, samp_seg, overlap, flag);
Rue = xcorr(ud,ew,'unbiased');
Rue = Rue(maxlag+1:end);
% Run = cum2x(ud,ns,maxlag, samp_seg, overlap, flag);
Run = xcorr(ud,ns,'unbiased');
Run = Run(maxlag+1:end);
figure
subplot(211),plot(Rue)
subplot(212),plot(Run)
for i = 2:in
    Rue(i) = Rue(i-1)*alfa + ud(i)*ew(i);
    Run(i) = Run(i-1)*alfa + ud(i)*ns(i);
    Tue = Rue(i);
    Tun = Run(i);
    temp = atan(Tue./Tun)*180/pi;
    if (Tue<0)&(Tun<0)
        az(i) = temp;
    elseif (Tue>0)&(Tun<0)
        az(i) = 180 - abs(temp);
    elseif (Tue>0)&(Tun>0)
        az(i) = 180 + abs(temp);
    elseif (Tue<0)&(Tun>0)
        az(i) = 360 - abs(temp);
    end
end
不知道怎么了，做得结果总是不对，不知道我的程序是什么地方错了，请来看看。郁闷我好久了，不好意思啊

[ 本帖最后由 wuqiong 于 2008-12-29 17:47 编辑 ]

wuqiong

附件中可以看到公式

yufeng

是出现运行问题 还是结果不正确 分量记录什么的是否已知

wuqiong

不是运行问题,是结果不对.