% Demonstracia overenia navrhu IIR filtra (vystup z programu Filter
% Express, ktory je pouzivany v predmete SPvT).
% Vysledkom navrhu su koeficienty bikvadov. Kaskadnym zapojenim
% bikvadov dostavame IIR filter s pozadovanou charakteristikou.
% Realizacia IIR pomocou bikvadov je vyhodna jednak z pohladu
% stability realizacia v prostredi s konecnou presnostou, ako aj
% vzhladom na optimalizaciu architektur DSP prave pre realizaciu IIR
% filtrov pomocou bikvadov.
%
% SPvT, M.D. 26-02-2007, v.1.01

% Parametre navrhu (pre rpogram Filter Express v.5.1):
% Fvz = 48 kHz  (Input Clock = vzorkovacia frekvencia)
% Quantize Coefficients
% Fixed Point Quantization (Integer Bits = 1, Fraction Bits = 15)
% IIR filter
% Bi-Linear Transform
% Elliptic
% Low Pass
% F1 = 5 kHz (pasmo priepustnosti)
% F2 = 7 kHz (pasmo tlmenia)
% Pass Band 1 dB
% Stop Band 60 dB

fvz = 48000;       % vzorkovacia frekvencia

% Vysledkom navrhu su tieto koeficenty (3 bikvady)
% Program Filter Express realizuje navrh pre prenosovu
% funkciu bikvadu v tvare
% H(z) = (B0 + B1*z^-1 + B2*z^-2)/(1 - A1*z^-1 - A2*z^2)

% bikvad 1
B01 = 0.02215576171875;
B11 = 0.020538330078125;
B21 = 0.02215576171875;
A11 = 1.6365966796875;
A21 = -0.70147705078125;

% bikvad 2
B02 = 0.213592529296875;
B12 = -0.174652099609375;
B22 = 0.213592529296875;
A12 = 1.57125854492188;
A22 = -0.823760986328125;

% bikvad3
B03 = 0.49176025390625;
B13 = -0.58294677734375;
B23 = 0.49176025390625;
A13 = 1.54629516601563;
A23 = -0.946868896484375;

% Pre simulaciu uvazujeme signal zlozeny z 2 harmonickych zloziek
N = 1024;               % pocet simulovanych vzoriek
f1 = 2000;              % frekvencia 1. harmonickej zlozky (z pasma priep.)
f2 = 8000;              % frekvencia 2. harmonckej zlozky (z pasma tlmenia)
a1 = 0.4;               % amplituda 1. zlozky
a2 = 0.4;               % amplituda 2. zlozky
n = (0:N-1);            % diskretny cas
% sucet zloziek je z int <-1,1)
x = a1*sin(2*pi*n*f1/fvz) + a2*sin(2*pi*n*f2/fvz);

% Na filtrovanie pouzijeme funkciu Filter
% !!! Matlabovska funkcia Filter implementuje v pripade bikvadu prenosovu
% funkciu
% H(z) = (B0 + B1*z^-1 + B2*z^-2)/(1 + A1*z^-1 + A2*z^2)
% t.j. ma v menovateli opacne znamienka a preto je potrebne u koeficientov
% A** zmenit znamienko. Pri pouziti koeficientov bikvadu v konkretnych
% funkciach (napr. vo funkciach VisualDSP) je potrebne vzdy preverit, pre
% aku prenosovu funkciu su cielove funkcie implementovane. Nespravne
% pouzitie znamienok vedie zvycajne k nestabilite celeho IIR filtra.

% Vystup 1. bikavdu
y1 = filter( [B01 B11 B21], [1 -A11 -A21], x);

% Vystup 2. bikavdu
y2 = filter( [B02 B12 B22], [1 -A12 -A22], y1);

% Vystup 3. bikavdu
y3 = filter( [B03 B13 B23], [1 -A13 -A23], y2);

figure(1);
plot(x);
ylabel('Amplituda');
xlabel('Diskretny cas');
title('Vstupny signal filtra');

figure(2);
plot(y3);
ylabel('Amplituda');
xlabel('Diskretny cas');
title('Vystupny signal filtra (s prech. dejom na zaciatku)');

% vypocet amplitudovej frekvencnej charakteristiky v M boddoch
M = 1000;

% 1. bikvad
h1 = freqz([B01 B11 B21], [1 -A11 -A21], M);

% 2. bikvad
h2 = freqz([B02 B12 B22], [1 -A12 -A22], M);

% 3. bikvad
h3 = freqz([B03 B13 B23], [1 -A13 -A23], M);

% Hodnoty h1,h2, h3 reprezentuju komplexne hodnoty frekvencnych
% charakteristik jednotlivych bikvadov. Kedze bikvady su zapojene v kaskade,
% celkova frekvencna charakteristika je dana sucinom (PO ZLOZKACH) hodnot
% h1, h2, h3.
h = h1.*h2.*h3;

% Vysledna amplitudova frkvencna charakteristika v dB
figure(3);
frek = (0:length(h)-1)*fvz/(2*length(h));       % pre zobrazenie do fvz/2
plot(frek,20*log10(abs(h)));
xlabel('Frekvencia');
ylabel('Amplituda [dB]');
title('Amplitudova frekvencna charakteristika IIR filtra');
grid;





