登录
注册
帮助
课程部分内容涉及地质记录的时间序列（time series）和频谱分析，主要应用快速傅立叶变化（FFT），很简单也很实用的工具。也许有人现在或者以后会涉及到时间序列分析相关内容，所以把这门课有关时间序列部分的作业（附件time series computing lab）发上来，使用的数据为俄罗斯Vostok station 的Ice Core中氘同位素相对丰度（附件VostokIceCore1999.dat）数据。我写的答案贴在下面，还有一些相关的图在附件中（图顺序乱了，鼠标放在图上可以显示编号）。
时间序列可以应用在地学的很多方面. Matlab recipes for earth sciences by Martin H. Trauth是本非常不错的参考书，完全可以自学的。
有关时间序列部分的lecture notes也贴在附件中，感兴趣的可以看下。
希望会有所帮助
_________________________________________________________________
% Time Series and Power Spectral Analysis of Vostok Ice Core
function [dRMS] = VostokIceCore(M)
icecore=M.data;
% Resampling of ice core age data
Age=linspace(0,420000,2048);
% Interpolation for deuterium
H2=interp1(icecore(:,2),icecore(:,3),Age);
% Calculate sampling frequency and Nyquist frequency
Fs=1/(Age(2)-Age(1));
Fnyq=1/2*Fs;
% Remove DC component
Deu=H2-mean(H2);
% Detrending
polyf=polyfit(Age,Deu,1);
Deu_new=Deu-polyval(polyf,Age);
% Perform FFT
Deu_fft=fft(Deu_new);
F=Fnyq*linspace(0,1,2048/2+1);
p=(abs(Deu_fft/2048)).^2;
P=[p(1) 2*p(2:2048/2) p(2048/2+1)];
figure;
subplot(211),plot(Age,Deu_new);
xlabel('Age (year)');ylabel('Delta D (‰)');
subplot(212),plot(F,P);
xlabel('Frequecy (1/year)');ylabel('Power');
% Plot Period vs Power using semi-logarithmic scale
figure
semilogx(ones(1,length(F)-1)./F(2:end),P(2:end),'linewidth',2);
xlabel('Period (year)');ylabel('power');
% Calculte Welch's periodogram using two overlapping segments length of 1536
Deu_1=Deu_new(1:1536);
Deu_2=Deu_new(2048-1536+1:end);
% Apply a Hann widow to each of the segments
Hann=1/2*(1-cos(2*pi*(0:1536-1)'/1536));
Hann_1=Deu_1'.*Hann;
Hann_2=Deu_2'.*Hann;
Hann_fft_1=fft(Hann_1);
Hann_fft_2=fft(Hann_2);
p_1=(abs(Hann_fft_1/1536)).^2;
p_2=(abs(Hann_fft_2/1536)).^2;
p_avg=(p_1+p_2)/2;
p_HW=[p_avg(1);2*p_avg(2:1536/2);p_avg(1536/2+1)];
% Compute the frequency
F_HW=Fnyq*linspace(0,1,1536/2+1);
figure
semilogx(ones(1,length(F_HW)-1)./F_HW(2:end),p_HW(2:end),'k','linewidth',2);
xlabel('Period (year)');ylabel('power');
% Filtering out frequencies greater than 1.25e-5
k1=find(F>1.25E-5);K1=k1(1);
X1=Deu_fft;
X1(K1:end-K1)=0;
x1=ifft(X1);
dRMS1=1-sum((Deu_new-real(x1)).^2)/sum((Deu_new-mean(Deu_new)).^2);
figure;
subplot(211),plot(Age,Deu_new);
xlabel('Age (year)');ylabel('Delta D (‰)');
subplot(212),plot(Age,real(x1),'linewidth',2);
xlabel('Age (year)');ylabel('Delta D after filter (F>1.25E-5) applied (‰)')
% Filtering out frequencies greater than 3.125e-5
k2=find(F>3.125E-5);K2=k2(1);
X2=Deu_fft;
X2(K2:end-K2)=0;
x2=ifft(X2);
dRMS2=1-sum((Deu_new-real(x2)).^2)/sum((Deu_new-mean(Deu_new)).^2);
figure;
subplot(211),plot(Age,Deu_new);
xlabel('Age (year)');ylabel('Delta D (‰)');
subplot(212),plot(Age,real(x2),'linewidth',2);
xlabel('Age (year)');ylabel('Delta D after filter (F>3.125E-5) applied (‰)')
% Filtering out frequencies greater than 5.5e-5
k3=find(F>5.5E-5);K3=k3(1);
X3=Deu_fft;
X3(K3:end-K3)=0;
x3=ifft(X3);
dRMS3=1-sum((Deu_new-real(x3)).^2)/sum((Deu_new-mean(Deu_new)).^2);
figure;
subplot(211),plot(Age,Deu_new);
xlabel('Age (year)');ylabel('Delta D (‰)');
subplot(212),plot(Age,real(x3),'linewidth',2);
xlabel('Age (year)');ylabel('Delta D after filter (F>5.5E-5) applied (‰)')
dRMS=[dRMS1;dRMS2;dRMS3];
end
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>> M=importdata('VostokIceCore1999.dat',' ',60);
>> dRMS=VostokIceCore(M)
dRMS =
0.506207640217901
0.774088779676613
0.903638085831479
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1）The sampling frequency is 0.004876(2400/420000) year-1; the Nyquist frequency is 0.002438 year-1.
2）See Figure 1.
3）The three largest peaks on the periodogram (105.1 Kyr, 38.2 Kyr and 28.0 Kyr,) represent eccentricity (~100Kyr), obliquity (~40 Kyr) and procession (~ 26 and 23 Kyr) orbital cycle respectively. (Figure 2);
4）After apply the Hann window and Welch’s method, the periodoram in Figure 3 is improved compared to that in Figure 2, because side lobes on the periodogram have been reduced. The three largest peaks on the periodogram are 105.1, 39.4 and 28.7 Kyr respectively. They are indicative of eccentricity, obliquity and procession forcing.
5）In the first case, frequencies that are greater than 1.25e-5 are filtered out from the original Fourier vector. After applying the inverse FFT to the filtered Fourier vector, the resulting time series (Figure 4) only contains periods greater than 80 Kyr (1/1.25e-5), which approximately represents the eccentricity component of the original delta D time series from the Vostok ice core. About 50.6% variation of the original data is accounted for by this synthetic (or by eccentricity forcing). Afterwards, frequencies that are greater than 3.125e-5 are filtered out from the original Fourier vector. Therefore, the resulting time series (Figure 5) only contains periods greater than 32 Kyr (1/3.125e-5), which approximately represents the eccentricity plus obliquity components of the raw delta D time series. About 77.4% variation of the original data is accounted for by eccentricity and obliquity. Finally, frequencies that are greater than 5.5e-5 are filtered out from the original Fourier vector. The resulting time series (Figure 6) only contains periods greater than 18 Kyr (1/5.5e-5), which approximately represents the eccentricity plus obliquity plus procession components of the raw delta D time series. About 90.4% variation of the raw delta D data from the Vostok ice core is accounted.
附加文件
-
Figure1.png (5.6K)
下载次数: 3 -
VostokIceCore1999.dat (133.86K)
下载次数: 3 -
Lab Assignment 3.pdf (321.24K)
下载次数: 4 -
Figure6.png (7.75K)
下载次数: 1 -
Figure5.png (7.69K)
下载次数: 0 -
Figure2.png (3.8K)
下载次数: 0 -
Figure3.png (3.9K)
下载次数: 0 -
Figure4.png (7.58K)
下载次数: 5 -
Lecture1.pdf (735.04K)
下载次数: 3 -
Lecture2.pdf (582.34K)
下载次数: 3 -
Lecture3.pdf (531.71K)
下载次数: 3 -
Lecture4.pdf (843.47K)
下载次数: 3
多重回复
