function w = ExtractWavelet(data, M)
% 功能:基于复赛谱提取输入地震数据data的地震子波w
% 输入:
% x - 地震信号,Nt×Nz
% M - 提取地震子波的长度
% 输出:
% w - 提取的地震子波,1×M
% 时间:2022/10/17
if nargin < 2
M = 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:Nt
c(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/2
w1(i) = w0(floor(M/2)+2-i);
end
for i = floor(M/2)+1:Nt
w1(i) = w0(i-floor(M/2));
end
pm = max(abs(w1));
w1 = w1./pm; % 归一化
w = w1(1:M); % 提取的子波w
end