四、以Matlab來做「頻譜分析」

以「台灣地區河川流量站資料庫」中綠水觀測站1956至2004年的資料來看，對每月平均流量的記錄對頻率加以分析（扣除1956、1959、1990、2000這四年的資料，因為當中有些月份的資料有欠缺），希望可以得到某些訊息。

----------------------------------------------------------------------------------------------------
close all
clear all
f=load('1956to2004month.txt')
x=f(1:end, 1)
y=f(1:end, 2)

figure(1); plot(x,y)
xlabel('year'); % X軸註解，以1957年為0，每隔一個月，增加1/12，1958年為2，依此類推
ylabel('流量 cms');% Y軸註解
title('1956～2004年綠水觀測站月平均流量圖');
[Pxx,f]=periodogram(y,[],1024,1);
magnitude=abs(Pxx);
figure(2); plot(f,magnitude), grid


----------------------------------------------------------------------------------------------------

圖一為每月的流量分佈

可以在頻譜分析圖中，看到兩個峰值（peak），第一個峰值出現在頻率＝0處，可以忽略不管。第二個峰值出現在頻率約為0.08之處，換算成週期為12.5年。
意即每月的平均流量多寡，約每隔12.5年，會有一次循環。

----------------------------------------------------------------------------------------------------
曾於課堂上聽老師介紹過－自然界的諸多訊號可以用sin函數來檢視，
因此，底下嘗試用sin函數來fit上述的資料。

假設 y= a*sin(b*x) 為時間與流量的關係函數
則以matlab來算迴歸函數
----------------------------------------------------------------------------------------------------

clear all ; close all ;
y0 = load('1956to2004month.txt');
xdata = y0(1:end, 1);
ydata = y0(1:end, 2);
f = fittype('a*sin(b*x)');
fit1 = fit(xdata,ydata,f,'StartPoint',[1 1]);
fdata = feval(fit1,xdata);
I = abs(fdata - ydata) > 1.5*std(ydata);
outliers = excludedata(xdata,ydata,'indices',I);
figure(1); plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*') %數值超過1.5個標準差以上，以粉紅色的*標示
text(4.5,-30,'Y = 5.244 × Sin ( 0.8625X )') %標記所算出的回歸方程式
----------------------------------------------------------------------------------------------------

所得之Sin函數為 Y = 5.244 × Sin ( 0.8625X )
算出週期約為7.6年( 2*pi/0.8625≒7.6)

意即每月的平均流量多寡，約每隔7.6年，會有一次循環。

推測與上一個頻譜分析所做出的結果不一樣的原因：（若有錯，歡迎指正）
1.兩個不一樣的方法，在資料處理上本來就有誤差。
2.用Sin函數去做迴歸時，信心程度只有95％，Coefficients (with 95% confidence bounds)
3.資料數仍不夠多，所找出的規則性不夠明確。
4.筆者本身的數學有待加強................




如果以不同的迴歸方式來計算迴歸線時，以下為例：
-------------------------------------------------------------------------------------------------
clear all ; close all ;
y0 = load('1956to2004month.txt');
xdata = y0(1:end, 1);
ydata = y0(1:end, 2);
f = fittype('a*sin(b*x)');
fit1 = fit(xdata,ydata,f,'StartPoint',[1 1]);
fdata = feval(fit1,xdata);
I = abs(fdata - ydata) > 1.5*std(ydata);
outliers = excludedata(xdata,ydata,'indices',I);
figure(1); plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*') %數值超過1.5個標準差以上，以粉紅色的*標示
text(4.5,225,'紅色線 ： Y = 5.244 × Sin ( 0.8625 X )') %標記所算出的回歸方程式
hold on ;
fit2 = fit(xdata,ydata,f,'StartPoint',[1 1],'Exclude',outliers);
fit3 = fit(xdata,ydata,f,'StartPoint',[1 1],'Robust','on');
plot(fit2,'g-')
plot(fit3,'b:')
text(4.5,200,'綠色線 ： Y = 1.879 × Sin ( 1.059 X )') %標記所算出的回歸方程式
text(4.5,175,'藍色線 ： Y = 4.459 × Sin ( 0.7496 X )') %標記所算出的回歸方程式

--------------------------------------------------------------------------------------------------
所得到的圖如下：
紅色線為之前第一種所用的迴歸方式；
綠色線為Exclude Dada所算得的迴歸線；
藍色線為以Robust迴歸（強韌迴歸）方式來計算迴歸線，所算得的週期（約為8.4年）略大於第一種方式（7.6年）。