请帮我将上述代码用matlab语言写出来
时间: 2023-12-30 09:01:51 浏览: 98
MATLAB代码
以下是用Matlab语言实现的代码:
% 读取面波数据并画图
clear all
st = irisFetch.Traces('II', 'MBAR', '00', 'BHZ', '2017-08-08 08:00:00', '2017-08-08 08:10:00');
dt = st(1).sampleRate; % 采样间隔
data = [];
scale = 0.05;
dx = 2;
figure('Position', [100, 100, 800, 600]);
for i = 1:length(st)
d = st(i).data;
data = [data, d];
t = (0:length(d)-1) * dt;
plot(t, d*scale+(i-1)*dx, 'b', 'LineWidth', 1);
hold on;
end
xlabel('Time (s)');
ylabel('Offset (m)');
xlim([t(1), t(end)]);
ylim([-10, (i-1)*dx+10]);
title('Surface wave data');
set(gca, 'FontSize', 12, 'Box', 'on');
saveas(gcf, 'Surface_wave.png');
hold off;
% 二维FFT
d = data.';
n = length(d(1, :)); % 时间点数
m = length(d(:, 1)) * 5; % 空间方向采样点数,m增大可以让FK谱光滑一点,以达到插值效果
D = zeros(m, n);
D(1:length(d(:, 1)), :) = d; % 补零
% 时间采样率
fs = 1 / dt;
% 空间采样率
xs = 1 / dx;
% 频率 (赫兹)
f = (-n/2:n/2-1) * fs / n;
% 波数 (每米)
k = (-m/2:m/2-1) * 2*pi*xs / m;
% 二维FFT
fk = fft2(D);
% 作图
pmin = -10;
P = abs(fftshift(fk));
P = P / max(P(:));
P = 10 * log10(P);
P2 = abs(fk);
P2 = P2 / max(P2(:));
P2 = 10 * log10(P2);
figure('Position', [100, 100, 1000, 800]);
subplot(2, 2, 1);
pcolor(f, k, P2);
shading interp;
colormap('jet');
caxis([pmin, 0]);
xlabel('Frequency (Hz)');
ylabel('Wave number (2\pi/m)');
set(gca, 'FontSize', 12, 'Box', 'on');
subplot(2, 2, 2);
pcolor(f, k, P);
shading interp;
colormap('jet');
caxis([pmin, 0]);
hold on;
plot([f(n/2+1), f(end), f(end), f(n/2+1), f(n/2+1)], ...
[k(1), k(1), k(m/2+1), k(m/2+1), k(1)], ...
'r--', 'LineWidth', 2);
xlabel('Frequency (Hz)');
ylabel('Wave number (m^{-1})');
set(gca, 'FontSize', 12, 'Box', 'on');
subplot(2, 2, 3);
pcolor(f(n/2+1:end), k(1:m/2), P(1:m/2, n/2+1:end));
shading interp;
colormap('jet');
caxis([pmin, 0]);
xlabel('Frequency (Hz)');
ylabel('Wave number (m^{-1})');
set(gca, 'FontSize', 12, 'Box', 'on');
subplot(2, 2, 4);
pcolor(f(n/2+1:end), abs(k(m/2+1:end)), P(1:m/2, n/2+1:end).');
shading interp;
colormap('jet');
caxis([pmin, 0]);
colorbar('FontSize', 12);
xlabel('Frequency (Hz)');
ylabel('Wave number (m^{-1})');
xlim([0, 100]);
set(gca, 'FontSize', 12, 'Box', 'on');
sgtitle('FK Spectrum', 'FontSize', 14, 'FontWeight', 'bold');
saveas(gcf, 'FK_spectrum.png');
阅读全文