%Appl3_4_1
load zdsud.dat;                          %讀取數據序列
x=zdsud;
dt=0.02; N=length(x);  %由于數據的采樣間隔0.02秒
n=0:N-1;t=n*dt;f=n/(N*dt);   %時間序列和頻率序列
y=fft(x);   
f1=0.6;f2=1.8;    %要濾去頻率的上限和下限
yy=zeros(size(y));  %設置與y相同元素的數組
for m=0:N-1   %將未在頻率范圍內的振動濾去
   if(m/(N*dt)>f1&m/(N*dt)<f2)...   %小于Nyquist頻率的保留頻率范圍
|(m/(N*dt)>(1/dt-f2)&m/(N*dt)<(1/dt-f1))  %大于Nyquist頻率的保留頻率范圍
      yy(m+1)= y(m+1);       %保留頻率范圍內的振動振幅不變
   else
      yy(m+1)=0.;     %其余頻率范圍的振動振幅為零
   end
 end
subplot(2,1,1),plot(t,zdsud);        %繪制原始值序列
xlabel('時間/s'),title('原始信號')
grid on; 
subplot(2,1,2),plot(t,real(ifft(yy))) %繪制濾波后的數據運用ifft變換回時間域并繪圖。
title('FFT濾波后的信號圖');
xlabel('時間/s')
grid on;