%Appl7_3
%用FIR數字濾波器實現對數字信號的濾波
load grbx3.txt                %讀取數據序列
Xt=grbx3;                     %把數據賦值給變量
Fs=50;                        %設定采樣率
dt=1/Fs;                      %計算采樣間隔
n=1:length(Xt);
Nn=length(Xt);                 %序列長度
t=n/Fs;                        %時間序列
Fcs1=3;  Fcs2=23;   %設置的通帶和阻帶邊界頻率，
Ws1=Fcs1/(Fs/2);  Ws2=Fcs2/(Fs/2);              %轉換為標準頻率
Wn=[Ws1 Ws2];                    %通帶、阻帶頻率(為標準頻率)
N=400;           %濾波器階數為400
pa=(N-1)/2/Fs;                 %計算相位延遲
b=fir1(N,Wn,'stop');           % 設計FIR帶阻濾波器
figure(1);
[H,f]=freqz(b,1,Nn,Fs);    %求出濾波器幅頻相頻特性
subplot(2,1,1),plot(f,20*log10(abs(H)))
xlabel('頻率/Hz');ylabel('振幅/dB');grid on;
subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))
xlabel('頻率/Hz');ylabel('相位/^o');grid on;
figure(2);
subplot(2,1,1),plot(t,Xt);          %繪制原始波形序列
xlabel('時間/s');ylabel('振幅');title('濾波前信號') 
grid on;

Yt=filter(b,1,Xt);                  %對信號進行濾波
t=t-pa;                            %計算相位延遲
subplot(2,1,2),plot(t,Yt);           %繪制濾波后波形序列
ylim([-600 200])
xlabel('時間/s');ylabel('振幅');title('濾波后信號');
xlim([0 max(t)+1]);grid on;