Matlab | 基于复赛谱提取地震数据的地震子波

function w = ExtractWavelet(data, M)
% 功能：基于复赛谱提取输入地震数据data的地震子波w
% 输入：
%       x     -  地震信号,Nt×Nz
%       M     -  提取地震子波的长度
% 输出：
%       w     -  提取的地震子波,1×M
% 时间：2022/10/17if nargin < 2M = 101;
end[Nt, Nz] = size(data);
z1 = zeros(Nt, Nz);
% 多道平均对数分解法
for k = 1:Nz x = data(:,k);y = fft(x);     % 做傅里叶变换z = log(y);     % 并取其对数得到对数谱z1(:,k) = ifft(z);   % 再做傅里叶反变换
end
z2 = sum(z1,2)/Nz;
% 设计一个低通滤波时窗c
% c = [ones(1,15),zeros(1,Nt-30),ones(1,15)]';  
c = zeros(Nt, 1);
for i = 1:Ntc(i) = 0.5*(1+cos(2*pi*(i-1)/Nt));
end
wn = z2.*c;     % 截取复赛谱时间域的前面接近于0的一段数据
wn = fft(wn);   % 并进行傅里叶变换，得到对数谱
wn = exp(wn);   % 再取其指数，得到子波谱
w0 = real(ifft(wn));  % 然后做傅里叶反变换得到时间域的地震子波w1 = zeros(Nt,1);
for i = 1:M/2w1(i) = w0(floor(M/2)+2-i);
end
for i = floor(M/2)+1:Ntw1(i) = w0(i-floor(M/2));
end
pm = max(abs(w1));
w1 = w1./pm; % 归一化w = w1(1:M); % 提取的子波wend