function [f,v,A_g]=HRLRT_final(seis,dt,x,fmin_plot,fmax_plot,vmin_plot,vmax_plot)
seis=seis';%將地震記錄轉(zhuǎn)置(本程序分析的數(shù)據(jù)縱軸為道數(shù)，橫軸為采樣點數(shù))
n=size(seis,1);%道數(shù)
l=size(seis,2);%采樣點數(shù)
T=dt*l;%采樣時長
df=1/T;%傅里葉變換后的頻率間隔
fmin = 0;%最小頻率
fmax = fmax_plot;%最大頻率
f=[fmin:df:fmax];%掃描頻率點組
fn=length(f);%掃描頻率點數(shù)
vmin = vmin_plot;%搜索最小速度
vmax = vmax_plot;%搜索最大速度
dv = 1;%搜索速度間隔
v=[vmin:dv:vmax];%掃描速度點組
vn=length(v);%掃描速度點數(shù)
p=1./v;%掃描慢度點組
pn=length(p);%掃描慢度點數(shù)
clear fx fv lo lof
for i=1:n
? ? fx(i,:)=fft(seis(i,:));
end
lo=zeros(n,pn); % 給復數(shù)矩陣預分配空間
for k=1:n
? ? for j=1:pn
? ? ? ? lo(k,j)=-1*sqrt(-1)*2*pi*p(j)*x(k);
? ? end
end
stacknum = 2; %預加權迭代次數(shù)
clear fp pcg
for i=1:fn
? ? pcg=diag(ones(1,pn));
? ? lof=exp(lo*df*i); %正變換算子
? ? for stack=1:stacknum
? ? ? ? fv(:,i)=lsqr(lof,fx(:,i),1e-6,1,pcg);
? ? ? ? resmean=mean(fv(:,i))/pn;
? ? ? ? resfact=sum((fv(:,i)-resmean).*conj(fv(:,i)-resmean))/(pn-1);
? ? ? ? for j=1:pn
? ? ? ? ? ? pcg(j,j)=1.0/(resfact+fv(j,i)*conj(fv(j,i)));
? ? ? ? end
? ? end
end
A = abs(fv);%求取復數(shù)的模
%【歸一化】
A_g = A;
% for i = 1:size(A_g,2)
%? ? A_g(:,i)=A_g(:,i)/norm(A_g(:,i));
% end
for i = 1:size(A_g,2)
? ? A_g(:,i)=A_g(:,i)./max(A_g(:,i));
end
end
