Projektowanie filtrów
Autorzy: Krzysztof Duda, Przemysław Korohoda
Pliki: rozdzial_04.zip
Przykład 4.1
Obliczanie charakterystyki amplitudowej oraz położenia zer i biegunów zadanej transmitancji H(s) i H(z) (plik przyklad_4_1.m).% Przykład 4.1:
% m-plik: przyklad_4_1.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./KD.
close all, clear all, clc
Ba = [1]; %współczynniki wielomianu licznika H(s)
Aa = [1 sqrt(2) 1]; %współczynniki wielomianu mianownika H(s)
Bd = [0.12858 0]; %współczynniki wielomianu licznika H(z)
Ad = [1 -1.4225 0.55301]; %współczynniki wielomianu mianownika H(z)
%% odpowiedź częstotliwościowa H(s) w Matlabie funkcja freqs
w = 0:0.1:5; %zakres częstotliwosci w [rad/s]
s = j*w;
HBa = polyval(Ba,s);
HAa = polyval(Aa,s);
Ha = HBa./HAa;
%% odpowiedź częstotliwościowa H(z) w Matlabie funkcja freqz
Om = 0:0.05:pi; %zakres częstotliwosci w [rad]
z = exp(j*Om);
HBd = polyval(Bd,z);
HAd = polyval(Ad,z);
Hd = HBd./HAd;
%% zera i bieguny
[za,pa,ka] = tf2zp(Ba,Aa);
[zd,pd,kd] = tf2zp(Bd,Ad);
%% wykresy
figure
subplot(2,2,1)
plot(w, abs(Ha),'.-'),
xlabel('\omega (rad/s)'), ylabel('|H(s)|')
axis tight,box on
subplot(2,2,2)
plot(Om, abs(Hd),'.-')
xlabel('\Omega (rad)'), ylabel('|H(z)|')
axis tight,box on
subplot(2,2,3), hold on
v = [-1.3 1.3 -1 1]; %zakres osi
plot(real(pa), imag(pa),'xb')
legend('bieguny H(s)')
plot([0 0],[v(3) v(4)],'-k') %oś OY
plot([v(1) v(2)],[0 0],'-k') %oś OX
xlabel('Re'), ylabel('Im')
axis(v), box on, axis equal
subplot(2,2,4), hold on
plot(real(zd), imag(zd),'ob')
plot(real(pd), imag(pd),'xb')
legend('zera H(z)','bieguny H(s)')
plot([0 0],[v(3) v(4)],'-k') %oś OY
plot([v(1) v(2)],[0 0],'-k') %oś OX
O = -pi:0.01:pi; O = exp(j*O);
plot(real(O),imag(O),'-k') %okrąg
xlabel('Re'), ylabel('Im')
axis(v), box on, axis equal
% KONIEC PRZYKŁADU 4.1.
Przykład 4.2
Optymalny w sensie średniokwadratowym, dolnoprzepustowy filtr FIR. Przykładowe wykorzystanie funkcji w „Przykład 4.4.A” (plik fir_mse.m).function [h,a]=fir_mse(M,Omc);
% [h,a]=fir_mse(M,Omc);
%
% projektowanie optymalnego w sensie średniokwadratowym
% dolnoprzepustowego filtra FIR o symetrycznej odpowiedzi
% impulsowej typu I i długości N=2M+1
% Omc - krawędź pasma w radianach
%
% Przykład 4.2.
% m-plik (funkcyjny): fir_mse.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./KD.
dOm= pi/1e3; %krok całkowania
Om = -pi: dOm : pi-dOm;
NOm= length(Om);
iR = zeros(M+1,M+1); % macierz iR nie zależy od współczynników filtra
iR(1,1)=1;
for k=2:M+1
iR(k, k) = 2;
end
OmD = ones(size(Om)); %żądana charakterystyka amplitudowa
ind = find(abs(Om)>Omc);
OmD(ind) = 0;
c = zeros(M+1,1);
for k=0:M
c(k+1,1) = sum( OmD.*cos(k*Om) ) / NOm;
end
a = iR*c;
for k=1:M
h(k) = a(k+1)/2;
end
h = [fliplr(h) a(1) h];
% KONIEC PRZYKŁADU 4.2.
Przykład 4.3
Optymalny w sensie średniokwadratowym, dolnoprzepustowy filtr FIR z ważonym błędem. Przykładowe wykorzystanie funkcji w „Przykład 4.4.A” (plik fir_mse_W.m).function [h, a]=fir_mse_W(M,Omc,Wp,Ws);
% [h,a]=fir_mse_W(M,Omc,Wp,Ws);
%
% projektowanie optymalnego w sensie średniokwadratowym
% dolnoprzepustowego filtra FIR o symetrycznej odpowiedzi
% impulsowej typu I i długości N=2M+1
% z ważoną charakterystyką amplitudową
% Omc - krawędź pasma w radianach
% Wp, Ws - wagi w paśmie przepustowym i zaporowym
%
% Przykład 4.3.
% m-plik (funkcyjny): fir_mse_W.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./KD.
dOm= pi/1e3; %krok całkowania
Om = -pi: dOm : pi-dOm;
NOm= length(Om);
OmD= ones(size(Om)); %żądana charakterystyka amplitudowa
ind= find(abs(Om)>Omc);
OmD(ind) = 0;
OmW = Wp*ones(size(Om)); %funkcja ważąca
OmW(ind) = Ws*ones(size(ind));
R = zeros(M+1,M+1);
for n=0:M
for m=0:M
R(n+1,m+1)= sum( OmW.*cos(Om*n).*cos(Om*m) )/NOm;
end
end
iR = inv(R);
c = zeros(M+1,1);
for k=0:M
c(k+1,1) = sum( OmW.*OmD.*cos(k*Om) )/NOm;
end
a = iR*c;
for k=1:M
h(k) = a(k+1)/2;
end
h = [fliplr(h) a(1) h];
% KONIEC PRZYKŁADU 4.3.
Przykład 4.4
Optymalny w sensie minimum błędów maksymalnych, dolnoprzepustowy filtr FIR z ważonym błędem. Przykładowe wykorzystanie funkcji w „Przykład 4.4.A” (plik fir_min_max.m).function [h, a, delt_teor]=fir_min_max(M,Omp,Oms,Wp,Ws);
% [h,a,delt_teor]=fir_min_max(M,Omp,Oms,Wp,Ws);
%
% projektowanie optymalnego w sensie minimum błedów maksymalnych
% dolnoprzepustowego filtra FIR o symetrycznej odpowiedzi
% impulsowej typu I i długości N=2M+1
% z ważoną charakterystyką amplitudową
% Omp, Oms - pasmo przepustowe i zaporowe w radianach
% Wp, Ws - wagi w paśmie przepustowym i zaporowym
%
% Przykład 4.4.
% m-plik (funkcyjny): fir_min_max.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./KD.
N_eks = M+2; %liczba ekstremów
dOm = (Omp(2)+Oms(2)-Oms(1))/N_eks;
Om = 0:dOm:pi;
Np=round(Omp(2)/dOm); %liczba ekstremów w pasmie przepustowym
Ns=N_eks-Np; %liczba ekstremów w pasmie zaporowym
Omk = [Om(1:Np) Om(end-Ns+1:end)]; %częstotliwosci ekstremów
for KK=1:1e2 %pętla może być przerwana funkcją break w linni 48
R = zeros(N_eks,N_eks-1);
c = ones(N_eks,1);
ind_c= find(Omk>=Oms(1));
c(ind_c) = 0;
znak = (-1).^([1:N_eks]);
w = znak/Wp;
w(ind_c) = znak(ind_c)/Ws;
w=w(:);
for k=0:N_eks-1
for m=0:M
R(k+1,m+1)=cos(m*Omk(k+1));
end
end
a_delt = inv([R w])*c;
%% ch-ka filtra
Om = 0:pi/1e3:pi;
Hw = zeros(size(Om));
for n=0:M
Hw = Hw + a_delt(n+1)*cos(n*Om );
end
pochodna = diff(Hw);
zmiana_znaku = [0 pochodna].*[pochodna 0];
ind1= find(pochodna==0);
ind2= find(zmiana_znaku<0);
ind = [ind2 ind1]; %indeksy do maksimów
%% obliczenie błędu
delt_teor = abs(a_delt(N_eks));
Om_eks = Om(ind);
Hw_eks = Hw(ind);
ind_eks= find(Om_eks>=Oms(1));
delt = Wp*(1-Hw_eks);
delt(ind_eks)= Ws*(0-Hw_eks(ind_eks)); %pasmo zaporowe
delt_max=max(abs(abs(delt)));
if abs(delt_max-delt_teor)<1e-9;
disp(['Wykonano ' num2str(KK) ' iteracji'])
break %wyjscie z pętli
end
if Wp>=Ws
Omk = [Om(ind) 0 pi Omp(2) Oms(1)]; %nowe częstotliwości
else
Omk = [Om(ind) 0 pi Oms(1) Omp(2)]; %nowe częstotliwości
end
% Omk = sort( Omk(1:N_eks),'ascend');
Omk = sort(Omk(1:N_eks));
Omk=Omk(end:-1:1);
end
for k=1:M
h(k) = a_delt(M-k+2)/2;
end
h = [h a_delt(1) fliplr(h)];
% KONIEC PRZYKŁADU 4.4.
Przykład 4.4.A
Przykład ten (plik przyklad_4_4_A.m) pokazuje jak można skorzystać z funkcji zawartych w przykładach 4.2, 4.3 oraz 4.4 i prezentuje porównanie otrzymanych filtrów w postaci ich charakterystyk amplitudowych.% Przykład 4.4.A:
% m-plik: przyklad_4_4_A.m
%
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./KD.
close all, clear all, clc
FontName = 'Times New Roman';
FontSize = 24;
Resolution = 600;
M = 10
N = 2*M+1
Omp = [0 0.45*pi ]; %pasmo przepustowe
Oms = [0.55*pi pi ]; %pasmo zaporowe
Omc = (Omp(2)+Oms(1))/2;
Wp = 1;
Ws = 100;
%% 1. MSE
h_1 = fir_mse(M,Omc);
%% 2. MSE wazony
Wp1 = 1;
Ws1 = 175;
h_2=fir_mse_W(M,Omc,Wp1,Ws1);
%% 3. MIN-MAX M=[10]
Wp2=1;
Ws2=10;
h_3=fir_min_max(M,Omp,Oms,Wp2,Ws2);
%% Ch-ki amplitudowe
[Hw_1, w_1] = freqz(h_1,1);
[Hw_2, w_2] = freqz(h_2,1);
[Hw_3, w_3] = freqz(h_3,1);
kreska=['-k';'-b';'-r'];
FontSize = 16;
figure, hold on
subplot(2,1,1), hold on
plot(w_1, 20*log10(abs(Hw_1)), kreska(1,:),'LineWidth',1)
plot(w_2, 20*log10(abs(Hw_2)), kreska(2,:),'LineWidth',1)
plot(w_3, 20*log10(abs(Hw_3)), kreska(3,:),'LineWidth',1)
%legend('MS','MSW','MinMax')
title(['N=' num2str(N)],'FontSize',FontSize,'FontName',FontName)
xlabel('\Omega (rad)','FontSize',FontSize,'FontName',FontName)
ylabel('|H(\Omega)| (dB)','FontSize',FontSize,'FontName',FontName)
axis tight, box on, grid on
v=axis; axis([v(1) 2 -2 2])
set(gca,'FontSize',FontSize,'FontName',FontName)
subplot(2,1,2), hold on
plot(w_1, 20*log10(abs(Hw_1)), kreska(1,:),'LineWidth',1)
plot(w_2, 20*log10(abs(Hw_2)), kreska(2,:),'LineWidth',1)
plot(w_3, 20*log10(abs(Hw_3)), kreska(3,:),'LineWidth',1)
legend('MS','MSW','MM','Location','SouthWest')
title(['N=' num2str(N)],'FontSize',FontSize,'FontName',FontName)
xlabel('\Omega (rad)','FontSize',FontSize,'FontName',FontName)
ylabel('|H(\Omega)| (dB)','FontSize',FontSize,'FontName',FontName)
axis tight, box on, grid on
v=axis; axis([v(1) v(2) -60 5])
set(gca,'FontSize',FontSize,'FontName',FontName)
% KONIEC PRZYKŁADU 4.4.A;
Przykład 4.5
Przykład (plik przyklad_4_5.m) przedstawia obliczenia wstępne przygotowujące do badania filtra interpolującego według (4.71) oraz – porównawczo – tego samego filtra z zastosowanym oknem Hamminga. Warto także porównać wyniki otrzymane dla innych okien dostępnych w pakiecie MATLAB. Przykłady wykresów można zobaczyć na rysunku 4.43. Należy zwrócić uwagę na szczególnie wnikliwe studium charakterystyki fazowej.% Przykład 4.5:
% m-plik: przyklad_4_5.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
M=11; n=0:M; d=0.25; % określamy rząd filtra, indeksy dla h[n] oraz opóźnienie ułamkowe;
h1=sinc(n-(M+1)/2+d); % wyliczamy h1[n] dla filtra według (4.71);
w=hamming(M+1); % wyliczamy wartości okna Hamminga;
h2=h1.*w'; % mnożymy h1[n] przez współczynniki okna;
[H1,W]=freqz(h1,1,10000); F=W/(2*pi); % wyznaczamy DTFT dla h1[n];
dF=F(2)-F(1); P1=unwrap(angle(H1)); % wyliczamy przyrost na osi F oraz ch-kę fazową;
D10=-P1(2:end)./F(2:end); % wyliczamy opóźnienie fazowe;
D1=D10/(2*pi); % ... i przeliczamy na opóźnienie w indeksach;
G10=-(P1(2:end)-P1(1:end-1))/dF; % wyliczamy opóźnienie grupowe;
G1=G10/(2*pi); % ... i przeliczamy na opóźnienie w indeksach;
figure(1);
plot(0:M,h2,'bo-'); grid on; axis tight
xlabel('n'); legend('elementy odpowiedzi impulsowej h1[n]');
figure(2);
subplot(4,1,1); plot(F,abs(H1),'b.-'); grid on;
legend('DTFT: modul(H1)',0);
subplot(4,1,2); plot(F,angle(H1)/pi,'b.-'); grid on;
legend('DTFT: faza(H1)/pi',0);
subplot(4,1,3); plot(F(1:end-1)+dF/2,D1,'b.-'); grid on;
legend('opoznienie fazowe dla h1',0);
subplot(4,1,4); plot(F(1:end-1)+dF/2,G1,'b.-'); grid on;
legend('opoznienie grupowe dla h1',0); xlabel('F');
% analogicznie dla h2, a następnie wykresy - ale to już proszę wykonać samodzielnie...
% KONIEC PRZYKŁADU 4.5.
Przykład 4.6
Przykład (plik przyklad_4_6.m) pokazuje zrealizowane w pakiecie MATLAB porównanie dwóch metod interpolacji wynikającej z K-krotnego zwiększenia tempa próbkowania sygnału wejściowego: a) za pomocą banku filtrów polifazowych typu sinc, b) za pomocą pojedynczej filtracji wykonanej po nadpróbkowaniu.% Przykład 4.6:
% m-plik: przyklad_4_6.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
M=10; n=0:M; K=32; % rząd filtra , indeksy h[n] oraz krotność nadpróbkowania;
dd=1/K; % ułamkowy krok opóźnienia w dziedzinie indeksów;
k=0; % uruchamiamy licznik pętli;
for d=0:dd:1-dd; k=k+1; % początek pętli oraz zwiększenie licznika "k" w każdym obiegu;
hh(k,:)=sinc(n-M/2+d); % w wierszach macierzy hh umieszczamy kolejne odp. impulsowe;
h0(k:K:K*(M+1))=hh(k,:); % te same odp. impulsowe układamy naprzemiennie w h0[n];
end;
N=2^8; x=rand(1,N); % to jest dyskretny "dowolny sygnał";
for k=1:K,
y1(k:K:K*N)=filter(hh(k,:),1,x); % wyliczamy poprzesuwane wyniki filtracji
% (patrz rys.4.42),
% układając próbki naprzemiennie w ciągu wynikowym;
end;
% sposób alternatywny:
x2=zeros(1,K*N); x2(1:K:K*N)=x; % ciąg x2, to nadpróbkowany K-krotnie sygnał x;
y2=filter(h0,1,x2); % wyznaczamy wynik filtracji za pomocą h0;
err=max(abs(y1(:)-y2(:))), % porównanie, czy to na pewno to samo;
figure(1);
subplot(2,1,1);
plot(y1,'r.-'); grid on; hold on;
plot(y2,'bo-');
title('Warto ten wykres rozciagnac');
legend('sygnal y1','sygnal y2'); xlabel('n');
subplot(2,1,2);
plot(h0,'ko-'); grid on;
legend('odp.imp. h0[n]'); xlabel('n');
figure(2);
mesh(hh); axis tight;
title('wszystkie odpowiedzi impulsowe razem'); xlabel('indeks n'); ylabel('numer odpowiedzi imp.');
% KONIEC PRZYKŁADU 4.6.
Przykład 4.7
Przykład (plik przyklad_4_7.m) pokazuje wzajemne powiązania między najważniejszymi wielkościami występującymi w opisie zagadnienia oraz konfrontuje wynik zastosowania wzoru (4.74) oraz dwóch sposobów wykorzystania funkcji firrcos.% Przykład 4.7:
% m-plik: przyklad_4_7.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all; % porządki wstępne;
Tb=0.01; % czas trwania pojedynczego symbolu;
M=200; K=20; % założony rząd filtra oraz liczba próbek na symbol;
a=0.5; % parametr pasma przejsciowego filtra;
fpr=K/Tb; % częstotliwość próbkowania;
Df=a/Tb; % szerokość pasma przejsciowego filtra;
h0 = firrcos( M, 1/(2*Tb), Df, fpr); % MATLAB: 1. sposób;
h1 = firrcos( M, 1/(2*Tb), a, fpr,'rollof'); % MATLAB: 2. sposób;
%% wzór (4.74):
n=0:M; Nd=M/2; % indeksy oraz wartość opóźnienia;
h2 = sinc((n-Nd)/K).* cos(pi*a*(n-Nd)/K)./(1-4*a^2*(n-Nd).^2/K^2);
n0=find((1-4*a^2*(n-Nd).^2/K^2)==0); % wzór może dać dzielenie 0/0;
h2(n0)=0; % wartość graniczna, gdy ze wzoru wynika 0/0;
h2=h2/K; % normalizacja;
err_0=max(abs(h0-h1)), % weryfikacja nr 1;
err_1=max(abs(h1-h2)), % weryfikacja nr 2;
[H,W]=freqz(h2,1,1000); % wyliczenie DTFT, czyli charakterystyk filtra;
F=W/(2*pi); f=F*fpr; % przeliczenie W na F (częstotl. unormowaną) oraz na f [Hz];
figure(1);
plot(0:M,h2,'bo-'); grid on;
legend('odp. impulsowa h2'); xlabel('n');
figure(2);
subplot(211);
plot(f,abs(H),'b.-'); grid on; legend('DTFT: modul(H)',0); title('Warto powiekszyc wybrane fragmenty');
subplot(212);
plot(f,angle(H)/pi,'r.-'); grid on; legend('DTFT: faza(H)/pi',0); xlabel('f [Hz]');
% KONIEC PRZYKŁADU 4.7.
Przykład 4.8
W przykładzie (plik przyklad_4_8.m) przedstawiono opisaną w rozdziale procedurę projektowania filtrów interpolowanych z wykorzystaniem filtrów FIR. Wyliczone odpowiedzi impulsowe i widma warto obejrzeć w postaci wykresów.% Przykład 4.8:
% m-plik: przyklad_4_8.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clear; close all; clc; % przygotowanie "biurka" do pracy;
Fg=0.01; Fz=0.012; % częstotliwości graniczne filtra;
A_dB=60; % tłumienie w paśmie zaporowym w dB;
NF=100000; % liczba punktów na osi F do wylicz. ch-k;
M=round(1/(Fg+Fz+sqrt(Fz-Fg))), % wyliczamy krotność nadpróbkowania;
Nref=round(A_dB/(22*(Fz-Fg))), % wyliczamy rząd filtra "ref";
Wref=[1,10]; % wagi dla procedury "remez";
Aref=[1,1,0,0]; % amplitudy docelowe dla funkcji "remez";
Fref=[0,Fg,Fz,0.5]*2; % wektor częstotl. dla funkcji "remez";
href=remez(Nref,Fref,Aref,Wref); % wyznaczamy odpowiedź imp. filtra "ref";
[Href,W]=freqz(href,1,NF); F=W/(2*pi); % wyliczamy ch-ki filtra "ref";
Fprot=[0,Fg*M,Fz*M,0.5]*2; % wektor częstotl. dla funkcji "remez";
Nprot=floor(Nref/M), % wyliczamy rząd filtra "prot";
hprot=remez(Nprot,Fprot,Aref,Wref); % wyznaczamy odp. imp. filtra "prot";
hnp=zeros(1,M*Nprot+1); hnp(1:M:end)=hprot; % nadpróbkowanie;
Hnp=freqz(hnp,1,NF); % ch-ki filtra po nadpróbkowaniu (F=Fref);
Fintrp_z=1/M-Fz; Fintrp_g=Fg; % częstotliwości graniczne dla filtra "intrp";
Fintrp=[0,Fintrp_g,Fintrp_z,0.5]*2; % wektor częst. dla funkcji "remez";
Nintrp=round(A_dB/(22*(Fintrp_z-Fintrp_g))), % wyliczamy rząd filtra "intrp";
hintrp=remez(Nintrp,Fintrp,Aref,Wref); % wyznaczamy odp. imp. filtra "intrp";
Hintrp=freqz(hintrp,1,NF); % ch-ki filtra "intrp" (F=Fref);
H_IFIR=Hnp.*Hintrp; % efekt zastępczy filtracji dwuetapowej;
Oszczednosc_w_procentach=100*(Nref-Nintrp-Nprot)/Nref,
figure(1);
subplot(2,1,1);
plot(F,abs(H_IFIR),'b.-'); grid on; hold on;
plot(F,abs(Href),'r.-');
xlim([0,0.012]);
title('DTFT: modul(H)'); legend('H_I_F_I_R','H_r_e_f',0);
subplot(2,1,2);
plot(F,angle(H_IFIR)/pi,'b.-'); grid on; hold on;
plot(F,angle(Href)/pi,'r.-');
xlim([0,0.012]);
title('DTFT: faza(H)/pi'); legend('H_I_F_I_R','H_r_e_f',0);
xlabel('F');
% KONIEC PRZYKŁADU 4.8.
Przykład 4.9
W przykładzie (plik przyklad_4_9.m) przedstawiono opisaną w rozdziale procedurę z filtrami IIR. Wyliczone widma warto obejrzeć w postaci wykresów.% Przykład 4.9:
% m-plik: przyklad_4_9.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 4. Projektowanie filtrów.
% Autorzy rozdziału: Krzysztof Duda, Przemysław Korohoda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all; % przygotowanie "biurka" do pracy;
Fg=0.01; Fz=0.012; % częstotliwości graniczne filtra;
A_dB=60; % tłumienie w paśmie zaporowym w dB;
NF=100000; % liczba punktów na osi F do wylicz. ch-k;
M=20; p=2*pi/(3*M); r=127/128; % parametry które warto zmieniać;
Nref=round(A_dB/(22*(Fz-Fg))), % wyliczamy rząd filtra "ref";
Wref=[1,10]; % wagi dla procedury "remez";
Aref=[1,1,0,0]; % amplitudy docelowe dla funkcji "remez";
Fref=[0,Fg,Fz,0.5]*2; % wektor częstotl. dla funkcji "remez";
href=remez(Nref,Fref,Aref,Wref);% wyznaczamy odpowiedź imp. filtra "ref";
[Href,W]=freqz(href,1,NF); F=W/(2*pi); % wyliczamy ch-ki filtra "ref";
Fprot=[0,Fg*M,Fz*M,0.5]*2; % wektor częstotl. dla funkcji "remez";
Nprot=floor(Nref/M), % wyliczamy rząd filtra "prot";
hprot=remez(Nprot,Fprot,Aref,Wref); % wyznaczamy odp. imp. filtra "prot";
hnp=zeros(1,M*Nprot+1); hnp(1:M:end)=hprot; % nadpróbkowanie;
Hnp=freqz(hnp,1,NF); % ch-ki filtra po nadpróbkowaniu (F=Fref);
A=1+2*cos(M*p); B=1+2*cos(p); % a1 i b1 opisują Hintrp1(z);
b1=[1,zeros(1,M-1),-r^M*A,zeros(1,M-1),r^(2*M)*A,zeros(1,M-1),-r^(3*M)*1];
a1=[1,-r*B,r^2*B,-r^3*1];
Hintrp1=freqz(b1,a1,NF);
C=2*cos(M*p/2); D=2*cos(p/2); % a2 i b2 opisują Hintrp2(z);
b2=[1,zeros(1,M-1),-r^M*C,zeros(1,M-1),r^(2*M)*1];
a2=[1,-r*D,r^2*1];
Hintrp2=freqz(b2,a2,NF);
E=5.2; % tę wartość dobieramy ręcznie;
hkor=[1,zeros(1,M-1),-E,zeros(1,M-1),1]; % odp. impulsowa filtra FIR korygującego;
Hkor=freqz(hkor,1,NF);
H_IFIR=Hnp.*Hintrp1.*Hintrp2.*Hkor; % efekt zastępczy filtracji;
figure(1);
subplot(4,1,1);
plot(F,abs(Hnp),'b.-'); grid on;; hold on;
xlim([0,0.012]); legend('DTFT: modul(H_n_p)',0);
subplot(4,1,2);
plot(F,abs(Hintrp1),'b.-'); grid on;; hold on;
xlim([0,0.012]); legend('DTFT: modul(H_i_n_t_r_p_1)',0);
subplot(4,1,3);
plot(F,abs(Hintrp2),'b.-'); grid on;; hold on;
xlim([0,0.012]); legend('DTFT: modul(H_i_n_t_r_p_2)',0);
subplot(4,1,4);
plot(F,abs(Hkor),'b.-'); grid on;; hold on;
xlim([0,0.012]); legend('DTFT: modul(H_k_o_r)',0);
xlabel('F');
figure(2);
subplot(2,1,1);
plot(F,abs(H_IFIR)/max(abs(H_IFIR)),'b.-'); grid on; hold on;
plot(F,abs(Href),'r.-');
xlim([0,0.012]);
title('DTFT: modul(H)'); legend('H_I_F_I_R','H_r_e_f',0);
subplot(2,1,2);
plot(F,angle(H_IFIR)/pi,'b.-'); grid on; hold on;
plot(F,angle(Href)/pi,'r.-');
xlim([0,0.012]);
title('DTFT: faza(H)/pi'); legend('H_I_F_I_R','H_r_e_f',0);
xlabel('F');
% KONIEC PRZYKŁADU 4.9.