%fliter_example.m
%單位沖擊信號通過帶阻濾波器
t = 1:40;                                  %信號的時間
x =zeros(size(t));                         
x(1) = 1;                                  %產生單位沖擊信號
[b,a] = butter(10,[0.3 0.7],'stop');       %設計帶阻濾波器
y = filter(b,a,x);                         %濾波                          
hold on;                                   %畫圖濾波前后的時域數據
stem(t,x,'marker','o');
stem(t,y,'marker','.');
xlabel('時間/ts(s)');
legend('原始數據','濾波后數據'); 
fx=fft(x);                                 %對單位沖擊信號進行傅立葉變換  
fx=fftshift(fx);
fy=fft(y);                                 %濾波后的信號進行傅立葉變換
fy=fftshift(fy);
figure;
subplot(2,1,1);                            %畫圖顯示濾波前后的數據頻譜
f=(t-20)/20;
plot(f,abs(fx),'b-',f,abs(fy),'b-.');      %畫信號的幅頻曲線
xlabel('數字頻率/\pi(rad)');
title('幅頻曲線');
legend('原始數據的幅頻曲線譜','濾波后數據的幅頻曲線');
subplot(2,1,2);  
plot(f,angle(fx),'b-',f,angle(fy),'b-.');      %畫信號的相頻曲線
xlabel('數字頻率/\pi(rad)');
title('相頻曲線');
legend('原始數據的相頻曲線','濾波后數據的相頻曲線');