%滤波器的指标
wp=2*pi*12;ws=2*pi*18;
Ap=1;As=15;
Fs=200; %抽样频率200Hz
%数字角频率
Wp=wp/Fs;
Ws=ws/Fs;
%确定阶数、截频及系数
[N1,wc1]=buttord(wp,ws,Ap,As,'s')
[num1,den1]=butter(N1,wc1,'s');
o1=linspace(0,wp,500);
o2=linspace(wp,ws,200);
o3=linspace(ws,50*2*pi,500);
H1=20*log10(abs(freqs(num1,den1,o1)));
H2=20*log10(abs(freqs(num1,den1,o2)));
H3=20*log10(abs(freqs(num1,den1,o3)));
figure(3)
%subplot(1,2,1)
plot([o1 o2 o3]/(2*pi),[H1 H2 H3])
xlabel('频率(Hz)')
title('6阶Butterworth低通滤波器的增益响应')
%模拟滤波器转数字滤波器
[numd1,dend1]=impinvar(num1,den1,Fs);
w=linspace(0,pi,512);
h1=freqz(numd1,dend1,w);
norm1=max(abs(h1));
numd1=numd1/norm1;
figure(1);
subplot(2,1,1)
plot(w/pi,20*log10(abs(h1)/norm1))
title('BW 幅度归一化的DF的幅度响应')
%计算Ap,As
w=[Wp Ws];
[h,ww]=freqz(numd1,dend1,w);
fprintf('Ap=%.4f\n',-20*log10(abs(h(1))));
fprintf('As=%.4f\n',-20*log10(abs(h(2))));
%将输入信号引入并进行滤波处理
[NUM]=xlsread('EEG signal (project 1)');
x=NUM(:,2);
L=length(x);
z1=filter(numd1,dend1,x);
figure(2);
plot(x);
title('输入信号的波形')
xlim([0,2010])
figure(5)
plot(z1)
title('经过BW滤波器后的输出波形')
xlim([0,2010])
y=fft(x);
y=fftshift(y);
figure(7)
plot((-L/2:L/2-1)*Fs/L,abs(y))
title('原信号的频谱幅度输出波形')
xlabel('频率/(Hz)');
y1=fft(z1);
y1=fftshift(y1);
L1=length(z1);
figure(8)
plot((-L1/2:L1/2-1)*Fs/L1,abs(y1));
title('经过BW滤波器的频谱幅度输出波形')
xlabel('频率/(Hz)');
%设计Chebyshev I 型滤波器
[N2,wc2]=cheb1ord(wp,ws,Ap,As,'s')
[num2,den2]=cheby1(N2,Ap,wc2,'s');
HH1=20*log10(abs(freqs(num2,den2,o1)));
HH2=20*log10(abs(freqs(num2,den2,o2)));
HH3=20*log10(abs(freqs(num2,den2,o3)));
figure(4)
%subplot(1,2,2)
plot([o1 o2 o3]/(2*pi),[HH1 HH2 HH3])
xlabel('频率(Hz)')
title('4阶Chebyshev I型低通滤波器的增益响应')
%模拟滤波器转数字滤波器
[numd2,dend2]=impinvar(num2,den2,Fs);
w=linspace(0,pi,512);
h2=freqz(numd2,dend2,w);
norm2=max(abs(h2));
numd2=numd2/norm2;
figure(1);
subplot(2,1,2)
plot(w/pi,20*log10(abs(h2)/norm1))
title('CB 幅度归一化的DF的幅度响应')
%计算Ap,As
w=[Wp Ws];
[hh,ww2]=freqz(numd2,dend2,w);
fprintf('Ap=%.4f\n',-20*log10(abs(hh(1))));
fprintf('As=%.4f\n',-20*log10(abs(hh(2))));
z2=filter(numd2,dend2,x);
figure(6);
plot(z2)
xlim([0,2010])
title('经过Chebyshev I滤波器的输出波形')
y2=fft(z2);
y2=fftshift(y2);
L2=length(z2);
figure(9)
plot((-L2/2:L2/2-1)*Fs/L2,abs(y2))
title('经过CB I滤波器的频谱幅度输出波形')
xlabel('频率/(Hz)');
% N=ceil(6.2*pi/(Ws-Wp));
% N=mod(N+1,2)+N;
% M=N-1;
% W=hanning(N)';
% Wc=(Wp+Ws)/2;
% %k=0:M;
% % hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);
% % hd(0.5*M+1)=hd(0.5*M+1)+1;
% % h=hd.*W;
% omega=linspace(0,pi,512);
% % mag=freqz(h,[1],omega);
% % figure(10)
% % plot(omega/pi,20*log10(abs(mag)))
% window2=hanning(N);
% B2=fir1(M,Wc/pi,'low',window2);
% H=freqz(B2,1,512);
% figure(10)
% plot(omega/pi,20*log10(abs(H)))
% h=B2.*W;
% z3=filter(h,1,x);
% figure(11)
% plot(z3)
% title('经过FIR滤波器后的输出波形')
% y3=fft(z3);
% figure(12)
% plot(y3)
% title('经过FIR滤波后的幅度谱')