%% 設(shè)計(jì)通帶為30HZ-60Hz的6階巴特沃斯帶通濾波器
% [filtered_signal,filtb,filta] = bandpass_butterworth(inputsignal,cutoff_freq,Fs,order)
% INPUTS:
% inputsignal = input time series
% cutoff_freqs = filter corner frequencies in the form [f1 f2]
% Fs = data sampling frequency
% order = order of Butterworth filter
% OUTPUTS: 
% filtered_signal = the filtered time series
% filtb, filta = filter numerator and denominator (optional)
clear;clc;
t=0:0.005:5;
T=0.005;%時(shí)間間隔
Fs=1/T;%采樣頻率
y1=5*sin(2*pi*10*t)+3*sin(2*pi*40*t)+6*sin(2*pi*80*t);
y2=3*sin(2*pi*40*t);
[filtered_signal,filtb,filta]=bandpass_butterworth(y1,[30 60],Fs,6);
figure(1)
subplot(3,1,1)
plot(t,y1)
title('原始信號(hào)')
subplot(3,1,2)
plot(t,filtered_signal)
title('濾波后信號(hào)')
subplot(3,1,3)
plot(t,y2)
title('3*sin(2*pi*40*t)')
%濾波器的幅頻相頻曲線
figure(2)
[h,w] = freqz(filtb,filta,256);
%[h,f] = freqz(filtb,filta,256,Fs);
subplot(2,1,1)
plot(w/pi,20*log10(abs(h)))
%plot(f,20*log10(abs(h)))
xlabel('\times\pi')
%xlabel('Hz')
ylabel('gain(dB)')
subplot(2,1,2)
plot(w/pi,angle(h))
%plot(f,angle(h))
xlabel('\times\pi')
%xlabel('Hz')
ylabel('phase (rad)')

%% 設(shè)計(jì)阻帶為30HZ-60Hz的6階巴特沃斯帶阻濾波器
% [filtered_signal,filtb,filta] = bandstop_butterworth(inputsignal,cutoff_freqs,Fs,order)
% INPUTS:
% inputsignal = input time series
% cutoff_freqs = filter corner frequencies in the form [f1 f2]
% Fs = data sampling frequency
% order = order of Butterworth filter
% OUTPUTS: 
% filtered_signal = the filtered time series
% filtb, filta = filter numerator and denominator (optional)
clear;clc;
t=0:0.005:5;
T=0.005;%時(shí)間間隔
Fs=1/T;%采樣頻率
y1=5*sin(2*pi*10*t)+3*sin(2*pi*40*t)+6*sin(2*pi*80*t);
y2=5*sin(2*pi*10*t)+6*sin(2*pi*80*t);
[filtered_signal,filtb,filta]=bandstop_butterworth(y1,[30 60],Fs,6);
figure(1)
subplot(3,1,1)
plot(t,y1)
title('原始信號(hào)')
subplot(3,1,2)
plot(t,filtered_signal)
title('濾波后信號(hào)')
subplot(3,1,3)
plot(t,y2)
title('5*sin(2*pi*10*t)+6*sin(2*pi*80*t)')
%濾波器的幅頻相頻曲線
figure(2)
[h,w] = freqz(filtb,filta,256);
%[h,f] = freqz(filtb,filta,256,Fs);
subplot(2,1,1)
plot(w/pi,20*log10(abs(h)))
%plot(f,20*log10(abs(h)))
xlabel('\times\pi')
%xlabel('Hz')
ylabel('gain(dB)')
subplot(2,1,2)
plot(w/pi,angle(h))
%plot(f,angle(h))
xlabel('\times\pi')
%xlabel('Hz')
ylabel('phase (rad)')