Wyznaczanie odpowiedzi częstotliwościowej filtru FIR rzędu N-1 dla pojedynczej wartości częstotliwości cyfrowej:
>> N=8;
>> h=rand(1,N);
>> f=0.25;
>> n=0:N-1;
>> s_we=exp(j*2*pi*f*n);
>> s_wy=filter(h,1,s_we);
>> H=s_wy/s_we;
Wyznaczanie transmitancji “Z” filtru typu FIR o zadanej odpowiedzi impulsowej dla pojedynczej wartości z:
>> f=0.5;
>> z=exp(j*2*pi*f)
>> h=rand(1,8);
>> p=0:length(h)-1;
>> H=sum(h.*(z.^p));
Wyznaczanie przybliżenia transformaty D-TFT ze wzoru definicyjnego (27):
>> f=-30:0.01:30;
>> T=1/20;
>> x=rand(1,32);
>> n=-10:21;
i następnie:
>> k=0; for ff=f, k=k+1; X(k)=sum(x.*exp(-j*2*pi*f*n*T)); end;
albo:
>> k=(1:length(f))’;
>> X(k)=exp(-j*2*pi*f’*n*T)*x’;
Przykłady projektowania filtrów typu IIR
Przykład 1:
Zaprojektować pasmowoprzepustowy cyfrowy filtr eliptyczny o pasmie przepustowym w przedziale częstotliwości 100 Hz do 300 Hz, przy częstotliwości Nyquist'a równej 1000 Hz (częstotliwość próbkowania 2 000 Hz). Rząd filtru 5. Tętnienia w pasmie przepustowym nie powinny przekraczać 3 dB, a tłumienie w pasmie zaporowym powinno wynosić co najmniej 70 dB.
Odp.:
>> [b,a] = ellip( 5, 3, 70, [100/1000, 300/1000]);
w celu sprawdzenia wyników:
>> F=logspace(1,3,512); - wektor 512 wartości rozłożonych logarytmicznie na odcinku od 10 do 1000
>> H1 = freqz( b, a, F, 2000 ); - w punktach określonych przez F, dla częstotliwości próbkowania 2000 Hz
>>[H2,W]=freqz(b,a,512,’whole’); - odpowiedź częstotliwościowa określona w 512 punktach okręgu jednostkowego pulsacji cyfrowej oraz odpowiedni wektor pulsacji cyfrowych;
>> figure(1);
>> semilogx(F, 20*log10(abs(H1)) ); grid;
>> figure(2);
>> plot(W,20*log10(abs(H2))); grid;
dla porównania obu powyższych wyników można podać polecenie:
>> figure(1);
>> hold on
>> semilogx(W(1:257)*1000/pi,20*log10(abs(H2(1:257))),’r’);
natomiast po komendach:
>> figure(3);
>> semilogx( F, unwrap(angle(H1))); grid;
dostajemy wykres fazy właśnie zaprojektowanego filtru - warto ten wykres przemnożyć przez odpowiednio przeskalowaną amplitudę w celu wyeliminowania z wykresu fazy wartości odpowiadających niewielkim amplitudom:
>> hold on;
>> ind=find(floor(abs(H1)*100));
>> L=length(ind);
>> H1prog=zeros(1,512);
>> H1prog(ind)=ones(1,L);
>> FazaH1=angle(H1).*H1prog;
>> semilogx(F,unwrap(FazaH1),’m’);
Przykład 2:
Znaleźć rząd dla filtrów cyfrowych Butterwortha i eliptycznego oraz pasmo znormalizowane przy następujących wymaganiach: pasmo przepustowe między 1000Hz a 2000Hz, pasmo zaporowe zaczyna się o 500 Hz od wymienionych częstotliwości, częstotliwość próbkowania wynosi 20KHz, tętnienia w pasmie przepustowym max. 2dB, tłumienie w pasmie zaporowym - przynajmniej 70 dB.
Odp.:
>>[n,Wn] = buttord( [1000,2000]/10000, [500, 2500]/10000, 2, 70)
>> n = 16
>> Wn = 0.0983 0.2034
i dalej już projektujemy filtr typu Butterwortha stosując otrzymane wartości parametrów:
>>[b,a] = butter(n, Wn);
Analogicznie dla filtru eliptycznego dostaniemy:
>>[n,Wn] = ellipord( [1000,2000]/10000, [500, 2500]/10000, 2, 70)
>> n = 6
>> Wn = 0.1000 0.2000
>>[b,a] = ellip(n, 2, 70, Wn);
Widać, że różnica w rzędach tych dwóch typów filtrów dla tej samej specyfikacji w dziedzinie częstotliwościowej jest bardzo duża.
Przykład 3:
Przykład ma na celu wykazanie różnic w zastosowaniu różnych metod projektowania filtru cyfrowego typu IIR w oparciu o prototyp analogowy.
Projektowanie filtru analogowego typu Butterwortha rzędu N o pulsacji granicznej „1” (czyli znormalizowanego):
>> N=6;
>>[za,pa,ka]=buttap(N);
>>[ba,aa]=zp2tf(za,pa,ka);
>> Fg=1/(2*pi); - przeliczenie znormalizowanej pulsacji granicznej z radianow/sekundę na częstotliwość w Hz;
>> Fmax=2*Fg; - zakładamy, że maksymalna częstotliwość (częstotliwość Nyquista) wynosi 2 razy Fg;
>> Fs=2*Fmax; - przyjmujemy częstotliwość próbkowania (w Hz) równą podwojonej częstotliwości Nyquista;
projektowanie filtru cyfrowego za pomocą instrukcji „butter”:
>>[b1,a1]=butter(N,0.5); - częstotliwość graniczna wynosi 0.5, ponieważ przyjęliśmy, że Fmax=2*Fg;
b)projektowanie filtru cyfrowego z wykorzystaniem prototypu analogowego i transformacji dwuliniowej:
>>[b2,a2]=bilinear(ba,aa,Fs);
projektowanie filtru cyfrowego z wykorzystaniem prototypu analogowego i zachowaniem wartości próbek odpowiedzi impulsowej:
>>[b3,a3]=impinvar(ba,aa,Fs); - warto zauważyć, że otrzymane współczynniki są zespolone;
następnie korzystamy z funkcji freqz, w celu wyznaczenia próbek charakterystyki częstotliwościowej każdego z filtrów;
wektor częstotliwości pobieramy jedynie raz (256 wartości rozłożonych liniowo w zakresie od 0 do Fmax - Fmax w Hz odpowiada wartości częstotliwości cyfrowej równej „1”):
>>[H1,Fd]=freqz(b1,a1,256,Fs);
>> H2=freqz(b2,a2,256,Fs);
>> H3=freqz(b3,a3,256,Fs);
dla ćwiczenia wyznaczamy także logarytmicznie rozłożony wektor częstotliwości:
>> Flog=logspace(-2,log10(Fmax),256); - wektor 256 wartości w Hz rozłożonych logarytmicznie na odcinku od 0.01 do Fmax;
wzorcowa charakterystyka filtru analogowego (spróbkowana w punktach Flog) jest następująca:
>> Ha=freqs(ba,aa,Flog*2*pi); - częstotliwości w Hz zostały przeliczone na pulsację w radianach/sekundę, gdyż tak interpretuje wektor zmiennej niezależnej funkcja freqs;
Filtracja filtrem typu FIR z wykorzystaniem FFT (splotu kołowego)
Kompletne m-pliki do realizacji algorytmów overlap-add i overlap-save:
function y=over_add(x,h,NB);
% y=over_add(x,h,NB)
%
% x - ciag wejsciowy;
% h - odpowiedz impulsowa (Nh<<NB);
% NB - dlugosc bloku dla FFT - wtedy dlugosc wycinka "x": Nxb;
% y - wynik splatania;
% czesc inicjalizacyjna:
H=fft(h,NB);
Nh=length(h);
Nxb=NB-(Nh-1);
K=floor(length(x)/Nxb);
x=x(:)';
buf=zeros(1,Nh-1);
% algorytm wlasciwy:
for k=0:K-1,
n=1+Nxb*k:Nxb+Nxb*k;
xb=[x(n),zeros(1,Nh-1)];
XB=fft(xb);
YB=XB.*H;
yb=real(ifft(YB));
y(n)=[buf+yb(1:Nh-1),yb(Nh:Nxb)];
buf=yb(Nxb+1:NB);
end;
function y=over_save(x,h,NB);
% y=over_save(x,h,NB)
%
% x - ciag wejsciowy;
% h - odpowiedz impulsowa (Nh<<NB);
% NB - dlugosc bloku dla FFT - wtedy dlugosc wycinka "x": Nxb;
% y - wynik splatania;
% czesc inicjalizacyjna:
H=fft(h,NB);
Nh=length(h);
Nxb=NB-(Nh-1);
K=floor(length(x)/Nxb);
if NB+Nxb*(K-1)>length(x), K=K-1;end; % na wszelki wypadek
x=x(:)';
buf=zeros(1,Nh-1);
% algorytm wlasciwy:
for k=0:K-1,
n=1+Nxb*k:NB+Nxb*k;
% ny=Nh+Nxb*k:(Nh-1)+Nxb+Nxb*k; % po to by nie bylo roznicy w przesunieciu
ny=1+Nxb*k:Nxb+Nxb*k; % pojawi sie przesuniecie o "Nh-1"
xb=x(n);
XB=fft(xb);
YB=XB.*H;
yb=real(ifft(YB));
y(ny)=yb(Nh:Nh+Nxb-1);
end;