%Appl3_2
close all
load c818312b.dat;   %加載觀測(cè)數(shù)據(jù)
dt=1;   %采樣間隔為分鐘
x=c818312b;N=length(x);t=[0:N-1]*dt;
figure(1)  %繪出記錄數(shù)據(jù)隨時(shí)間的變化
plot(t/(24*60),x);
%以天為單位給出原始記錄，1天為24*60分鐘
xlabel('時(shí)間/天');
title('龍巖臺(tái)的地傾斜數(shù)據(jù)');
ylabel('振幅');
figure(2)
x=detrend(x);   %去掉長(zhǎng)期趨勢(shì)（這里為漸增的趨勢(shì)）
y=fft(x);   %對(duì)該數(shù)據(jù)進(jìn)行頻譜分析
plot([0:N-1]/(N*dt),abs(y)*2/N); %繪出振幅譜
xlabel('頻率/分鐘^-^1');
title('地傾斜數(shù)據(jù)的振幅譜');
ylabel('振幅');
xlim([0 1/(4*60)]);  %只繪出周期大于4小時(shí)（240分鐘）的振幅
hold on;plot([1/(26*60) 1/(26*60)],ylim, 'k:');
plot([1/(22*60) 1/(22*60)],ylim, 'k:');
plot([1/(13*60) 1/(13*60)],ylim, 'k:');
plot([1/(11*60) 1/(11*60)],ylim, 'k:');
plot([1/(8.5*60) 1/(8.5*60)],ylim, 'k:');
plot([1/(7.5*60) 1/(7.5*60)],ylim, 'k:');
text(1/(24*60), 2, '日潮', 'rotation',90)
text(1/(12*60), 2, '半日潮', 'rotation',90)
text(1/(8*60), 2, '三分之一潮', 'rotation',90);
hold off