%巴特沃茲Butterworth濾波器的離散濾波算法 模擬程序
%程序設計:劉輝
%pangelaca@yahoo.com.cn (email/msn)
%由于我在柏林工業大學學習,源程序的一些變量名由德語縮寫成.(但其實閱讀上沒有影響)
%各種參數都可以替換,產生不同階,不同截止頻率,不同通帶的Butterworth濾波器.
%程序會顯示其傳遞函數,振幅特性曲線,相位特性曲線,沖激響應,階躍響應等,涉及了信號處理的方方面面.
%將圖像顯示部分去掉,利用產生的分子系數和分母系數,可以實現直接對任何輸入信號的濾波.

%各參數說名:
%fd: 3DB截止頻率
%fs: 當下降到as(DB)時的通頻
%fT: 采樣頻率
%as: 下降通頻時的截止DB

fd=200;
fs=280;
fT=8000;
as=20;

wd=2*pi*fd;
ws=2*pi*fs;
As= 10^(-as/20);

N=ceil(0.5*log10(As^(-2)-1)/log10(ws/wd));

Ord=0:N-1;
Sx=wd*exp(j*pi*(0.5+1/2/N+Ord/N))

Zx=exp(Sx/fT);

KoefZx=poly (Zx);

k=abs(prod(1-Zx));
vpa(k)

Zahler=zeros(1,N);
Zahler(1)=k;

[Hz,w]=freqz(Zahler,KoefZx,2048);
subplot(3,2,1);
plot (abs(Hz));
title ('幅頻特性曲線 |H(z)|');

LogampHz=log(abs(Hz));
subplot(3,2,2);
plot(w,LogampHz);
title ('對數特性曲線 H(z)');

PhaHz=angle(Hz);
subplot(3,2,3);
plot(w,PhaHz);
title ('相位特性曲線');

ImpulsFolge=zeros(1,800);
ImpulsFolge(1)=1;

Impulsantwort=filter(Zahler,KoefZx,ImpulsFolge);
subplot(3,2,5);
bar (Impulsantwort);
title ('沖激響應 h(t)');

SprungFolge=ones(1,800);

Sprungantwort=filter(Zahler,KoefZx,SprungFolge);
subplot(3,2,6);
bar (Sprungantwort);
title ('階躍響應');