%Appl3_1
clf
load w9497d;   %加載觀測數據
x=w9497d;
dt=1; N=length(x);  %由于數據的采樣間隔為天
n=0:N-1;t=n*dt;f=n/(N*dt);   %時間序列和頻率序列
subplot(2,1,1),plot(n/365.25,x); %以年為單位給出記錄數據隨時間的變化,1年有365.25天
xlabel('時間/年');
title('福州臺電導率數據');
ylabel('MS/CM');
y=fft(x);   
f1=0.0;f2=0.005;    %要濾去頻率的上限和下限,即去掉頻率范圍為0~0.005天-1
yy=zeros(1,length(y));  %設置與y相同元素的數組
for m=0:N-1   %將頻率落在該頻率范圍及其大于Nyquist頻率的波濾去
   if(m/(N*dt)>f1&m/(N*dt)<f2)...   %小于Nyquist頻率的濾波范圍
|(m/(N*dt)>(1/dt-f2)&m/(N*dt)<(1/dt-f1))  %大于Nyquist頻率的濾波范圍
%1/dt為一個頻率周期
      yy(m+1)=0.;       %置在此頻率范圍內的振動振幅為零
   else
      yy(m+1)=y(m+1);     %其余頻率范圍的振動振幅不變
   end
end
subplot(2,1,2),plot(n/365.25,real(ifft(yy))) %繪制濾波后的數據運用ifft變換回時間域并繪圖。
title('濾波后的逆Fourier輸出');
ylabel('MS/CM');
xlabel('時間/年')