przyklad z zielinskiego

% Ćwiczenie: Projektowanie nierekursywnych filtrów cyfrowych metodą okien z zastosowaniem okna
% Kaisera
% Podaj parametry filtra (np. pasmowozaporowego)
 fpr = 1000; % częstotliwość próbkowania [Hz]
 fd1 = 150; % częstotliwość dolna 1 [Hz]
 fd2 = 200; % częstotliwość dolna 2 [Hz]
 fg1 = 300; % częstotliwość górna 1 [Hz]
 fg2 = 350; % częstotliwość górna 2 [Hz]
 dp = 0.001; % oscylacje w paśmie przepustowym np. 0.1, 0.01, 0.001
 ds = 0.0001; % oscylacje w paśmie zaporowym np. 0.001, 0.001, 0.0001
 typ = 'bs'; % lp=LowPass, hp=HighPass, bp=BandPass, bs=BandStop
% Oblicz parametry okna
 if (typ=='lp') % Filtr LP
 df=fd2-fd1; fc=((fd1+fd2)/2)/fpr; wc=2*pi*fc;
 end
 if (typ=='hp') % Filtr HP
 df=fg2-fg1; fc=((fg1+fg2)/2)/fpr; wc=2*pi*fc; end
 if (typ=='bp' | typ=='bs') % Filtr BP lub BS
 df1=fd2-fd1; df2=fg2-fg1; df=min(df1,df2);
 f1=(fd1+(df/2))/fpr; f2=(fg2-(df/2))/fpr; w1=2*pi*f1; w2=2*pi*f2;
 end
% dp=((10^(Ap/20))-1)/((10^(Ap/20))+1); % opcjonalne przeliczenie Ap i As w [dB]
% ds=10^(-As/20); % na dp i ds
 d=min(dp,ds); A=-20*log10(d);
 if (A>=50) beta=0.1102*(A-8.7); end
 if (A>21 & A<50) beta=(0.5842*(A-21)^0.4)+0.07886*(A-21); end
 if (A<=21) beta=0; end
 if (A>21) D=(A-7.95)/14.36; end
 if (A<=21) D=0.922; end
 N=(D*fpr/df)+1; N=ceil(N); if(rem(N,2)==0) N=N+1; end
 N
 pause
 M = (N-1)/2; m = 1 : M; n = 1 : N;
% Wygeneruj okno
 w = besseli( 0, beta * sqrt(1-((n-1)-M).^2/M^2) ) / besseli(0,beta);
 plot(n,w); title('Funkcja okna'); grid; pause
% Wygeneruj odpowiedź impulsową filtra
 if (typ=='lp') h=2*fc*sin(wc*m)./(wc*m); h=[ fliplr(h) 2*fc h]; end % filtr LP
 if (typ=='hp') h=-2*fc*sin(wc*m)./(wc*m); h=[ fliplr(h) 1-2*fc h]; end % filtr HP
 if (typ=='bp') % filtr BP
 h = 2*f2*sin(w2*m)./(w2*m) - 2*f1*sin(w1*m)./(w1*m);
 h = [ fliplr(h) 2*(f2-f1) h];
 end
 if (typ=='bs') % filtr BS
 h = 2*f1*sin(w1*m)./(w1*m) - 2*f2*sin(w2*m)./(w2*m);
 h = [ fliplr(h) 1+2*(f1-f2) h];
 end
 plot(n,h); title('Odp impulsowa filtra'); grid; pause
% Wymnóż odp impulsową filtra z funkcją okna
 hw = h.* w;
 plot(n,hw,'b'); title('Odpowiedz impulsowa filtra z oknem'); grid; pause
% Charakterystyka częstotliwościowa
 NF = 1000; fmin = 0; fmax = fpr/2; % wartości parametrów charakterystyki
 f = fmin : (fmax-fmin)/(NF-1) : fmax; % częstotliwość
 w = 2*pi*f/fpr; % pulsacja
 HW = freqz(hw,1,w); % widmo Fouriera
 plot(f,abs(HW)); grid; title('Moduł'); xlabel('freq [Hz]'); pause
 plot(f,20*log10(abs(HW))); grid; title('Moduł dB'); xlabel('[Hz]'); ylabel('dB');
pause
 plot(f,unwrap(angle(HW))); grid; title('FAZA'); xlabel('[Hz]'); ylabel('[rad]'); pause