% 數字濾波器matlab的原代碼  
 　　 
function lvbo(Ua,Ub,choise)
%參考指令：lvbo(2*pi,10*pi,1/0/-1)
U1=min(Ua,Ub); 
U2=max(Ua,Ub);
Us=16*U2; 
T=2*pi/Us;
T_sum=4*max(2*pi/Ua,2*pi/Ub);
sum=T_sum/T;
t=T:T:T_sum;
x=sin(U1*t)+0.8*sin(U2*t);
X=DFT(x);
figure(1); subplot(221)
U=Us/sum:Us/sum:Us;
stem(U,abs(X));grid on
axis([Us/sum,Us/2,0,1.2*max(abs(X))])
title('原模擬信號采樣頻譜圖')
Ucd=U1+(U2-U1)*1/5;Usd=U2-(U2-U1)*1/5;
switch choise
case 1
Hz_ejw=IIR_DF_BW(Ucd,1,Usd,30,T,sum);
case -1
Hz_ejw=IIR_DF_CF(Ucd,1,Usd,30,T,sum);
case 0
Hz_ejw=FIR_DF_HM(U1,U2,T,sum);
otherwise
Hz_ejw=IIR_DF_BW(Ucd,1,Usd,30,T,sum);
end
Y=X.*Hz_ejw;
y=1/sum*conj(DFT(conj(Y)));
figure(1); subplot(224)
plot(t,real(y)); title('模擬信號濾波后');grid on
axis([0,T_sum,-max(real(y))*1.5,max(real(y))*1.5])
subplot(222); 
plot(t,x); hold on
axis([0,T_sum,-max(x)*1.2,max(x)*1.2])
x=sin(U1*t);
plot(t,x,':r');grid on
title('模擬信號濾波前')
function Hz_ejw=IIR_DF_BW(Ucd,Ap,Usd,As,t,sum)
% 巴特沃思濾波器
E=(10^(0.1*Ap)-1)^0.5;
V=(10^(0.1*As)-1)^0.5;
Wc=Ucd*t; Ws=Usd*t;
Ucd=Wc/t; Usd=Ws/t;
Uca=(2/t)*tan(Ucd*t/2); Usa=(2/t)*tan(Usd*t/2);
N=ceil(log10(V/E)/log10(Usa/Uca));
k=[1:2*N];
Spk=exp(j*(pi/2+(2*k-1)/(2*N)*pi));
i=find(real(Spk)<0);
Sk(1:N)=Spk(i);
den=real(poly(Sk'));
k0=polyval(den,0);
disp('模擬巴特沃思濾波器的歸一化統函數 Ha(s) 為')
tf(k0,den)
syms s z T;
den_jU=1;
s=s/Uca;
for i=1:N
den_jU=s^(N-i+1)*den(i)+den_jU;
end 
Ha_s=simple(1/den_jU);
H_z=subs(Ha_s,'s',(2/T)*((1-1/z)/(1+1/z)));
k=1:sum;w=(2*pi/sum)*k;
ejw=exp(j*w);
Hz_ejw=subs(H_z,{z,T},{ejw,t*ones(1,length(ejw))});
figure(1); subplot(223)
plot(w,abs(Hz_ejw)); grid on
title('巴特沃思低通濾波器')
axis([2*pi/sum,pi,-0.2,1.2*max(abs(Hz_ejw))])
function Hz_ejw=IIR_DF_CF(Ucd,Ap,Usd,As,t,sum)
% 切比雪夫低通濾波器
E=(10^(0.1*Ap)-1)^0.5;
V=(10^(0.1*As)-1)^0.5;
Wc=Ucd*t; Ws=Usd*t;
Ucd=Wc/t; Usd=Ws/t;
Uca=(2/t)*tan(Ucd*t/2); Usa=(2/t)*tan(Usd*t/2);
N=ceil(acosh(V/E)/acosh(Usa/Uca));;
A=1/E+(1/E^2+1)^0.5;
a=1/2*(A^(1/N)-A^(-1/N));
b=1/2*(A^(1/N)+A^(-1/N));
k=1:2*N;
Spk=-a*sin((2*k-1)/(2*N)*pi)+j*b*...
cos((2*k-1)/(2*N)*pi);
i=find(real(Spk)<0);
Sk(1:N)=Spk(i);
den=real(poly(Sk'));
k0=1;
disp('模擬切比雪夫低通濾波器的歸一化統函數 Ha(s) 為')
tf(k0,den)
if (rem(N,2)==1)
for i=1:N
k0=k0*(-Sk(i));
end
elseif ((rem(N,2))==0)
k0=1;
for i=1:N
k0=k0*(-Sk(i));
end
end
if (rem(N,2)==0)
k0=10^(-0.05*Ap)*k0;
end
k0=real(k0);
syms s z T;
den_jU=1;
s=s/Uca;
for i=1:N
den_jU=s^(N-i+1)*den(i)+den_jU;
end
Ha_s=simple(1/den_jU);
H_z=subs(Ha_s,'s',(2/T)*((1-1/z)/(1+1/z)));
k=1:sum;w=(2*pi/sum)*k;
ejw=exp(j*w);
Hz_ejw=subs(H_z,{z,T},{ejw,t*ones(1,length(ejw))});
figure(1); subplot(223)
plot(w,abs(Hz_ejw));grid on
title('切比雪夫低通濾波器')
axis([2*pi/sum,pi,-0.5,max(abs(Hz_ejw))])
function Hz_ejw=FIR_DF_HM(U1,U2,T,sum)
wp=U1*T;ws=U2*T;
kuan=ws-wp;
M=sum;
n=[0:1:M-1];
wc=(ws+wp)/2;
hd=H_D(wc,M);
window=hamming_m(M);
h_z=hd.*window;
Hz_ejw=DFT(h_z);
k=1:sum;w=(2*pi/sum)*k;
figure(1); subplot(223)
plot(w,abs(Hz_ejw));grid on
axis([2*pi/sum,pi,-0.2,1.2*max(abs(Hz_ejw))]);
title('海明窗函數低通濾波器')
function hd=H_D(wc,N)
M=(N-1)/2;
for k=-M:M
if k==0
hd(k+M+1)=wc/pi;
else
hd(k+M+1)=sin(wc*k)/(pi*k);
end
end
function wn=hamming_m(M)
n=0:M-1;
wn(n+1)=0.54-0.46*cos((2*pi*n)/(M-1));
function Xk=DFT(xn)
% 離散傅立葉變換,xn為原序列,Xk為DFT變換后的序列
N=length(xn);
n=0:N-1;
k=0:N-1;
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;