%%
% ZZS - 10.cvičení: Návrhy FIR a IIR filtrů, 
% R.Čmejla, 1.prosince 2014
%                                eeg, data.asc, spektrum, spektrum2
% příklad 2 a 7

%% Příklad 1:
clear all; close all
% a) Navrhněme dolní propust pro N = 50 , = 0,4 n w a pravoúhlé okno. Porovnejme
% lineární a logaritmickou amplitudovou frekvenční charakteristiku.
% b) Sledujme vliv délky okna na šířku přechodného pásma.
% c) Mějme dolní propust s parametry z bodu a).
% Porovnejme frekvenční charakteristiky pro různá okna.
% d) Navrhněme horní propust pro N = 33 a normovanou zlomovou frekvenci = 0,4 n w
% pomocí funkce fir1(). Proč MATLAB generuje varování?
% e) S použitím funkce grpdelay() vypočítejme skupinové zpoždění dolní propusti
% navržené v bodě a).
% a)
fn=0.4;				% mezni frekvence (normovana na Fs/2)
N= 50;				% delka filtru (pocet zpozdeni)
b=fir1(N,fn,boxcar(N+1));	% vypocet koeficientu oriznute idealni 
% impulzni odezvy
[H,f] = freqz(b,1);		% vypocet frekvencni charakteristiky
amp = 20*log10(abs(H));		% filtru a jeji zobrazeni
subplot(211)
plot(f/pi,amp,'k','Linewidth',1.5), grid on
xlabel('normovana frekvence [Hz]')
ylabel('amplituda [dB]')
subplot(212)
plot(f/pi,abs(H),'k','Linewidth',1.5), grid on
xlabel('normovana frekvence [Hz]')
ylabel('[-]')

% b)
figure
fn=0.5; N= 100;
b=fir1(N,fn,boxcar(N+1));
[H,f] = freqz(b,1);
plot(f/pi,abs(H),'k','Linewidth',1.5), grid on
xlabel('normovana frekvence [Hz]'), ylabel('[-]')
hold on
fn=0.5;					      
N= 10;						     
b=fir1(N,fn,boxcar(N+1));	
[H,f] = freqz(b,1);		     
plot(f/pi,abs(H),':k','Linewidth',1.5),
legend('N=100','N=10')

figure
clear b H
b(1,:) = fir1(50,0.4,boxcar(51)); 
[A, w] = freqz(b(1,:),1,512); 
H(1,:)=A';
b(2,:) = fir1(50,0.4,hanning(51)); 
H(2,:) = freqz(b(2,:),1,512)';
b(3,:) = fir1(50,0.4,hamming(51)); 
H(3,:) = freqz(b(3,:),1,512)';
b(4,:) = fir1(50,0.4,blackman(51)); 
H(4,:) = freqz(b(4,:),512)';
nadpis(1,:)=' Pravouhle ';
nadpis(2,:)='  Hannovo  ';
nadpis(3,:)='Hammingovo ';
nadpis(4,:)='Blackmanovo';
for k=1:4
    subplot(2,2,k)
    plot(w/pi,20*log10(abs(H(k,:)/max(H(k,:)))),'k','Linewidth',1.5);
    xlabel('\Omega/\pi'); ylabel('|H| [dB]'); 
    grid; axis([0 1 -110 0]); title(nadpis(k,:))
end;

% d)
b = fir1(33,0.4,'high'); [H w] = freqz(b,1,512);
figure;
plot(w/pi,20*log10(abs(H/max(H))),'k','Linewidth',2); grid;
xlabel('\Omega/\pi'); ylabel('|H| [dB]'); title('Horni propust');
figure
%e)
b = fir1(50,0.4); 
[H,w]=grpdelay(b,1,256);
plot(w,H,'k','Linewidth',2), grid
title('Skupinove zpozdeni'); ylabel(''); xlabel('\Omega/\pi');
axis([0 1 24 26]);

%% Příklad 2:
% Při zpracování EEG signálů se používá rozklad signálů do 
% frekvenčních pásem.
%           Pásmo         Amplituda                Typická činnost
%         < 0,5 Hz                                 Pohybové a oční artefakty
% delta 0,5 - 4 Hz        Velká                    Vlny vznikají v hlubokém spánku, tranzu
% théta   4 - 8 Hz        Střední                  Souvisí se stavem během denního snění, jsou
%                                                  příznačné pro některé psychické poruchy
% alfa    8 - 13 Hz       Malá                     Souvisí s relaxací
% beta   13 - 30 Hz       Nejmenší                 Souvisí s iritací, zlostí, frustrací,
%                                                  starostmi, duševním napětím; vznikají rovněž %                                                  při usilovném přemýšlení
%         > 30 Hz                                  EMG aktivita
% REM      60-70 Hz                                Rychlé pohyby očí (Rapid eye movement) 
%                                                  během spánku
% a) Navrhněme FIR filtry pro filtraci EEG signálů do 
%    jednotlivých frekvenčních pásem.
% b) Proveďme rozklad EEG signálu.
% c) Porovnejme vliv řádu filtru
% d) Porovnejme součet filtrovaných signálů s původním signálem


%% Příklad 3: Návrh FIR filtrů Remezovým algoritmem
% Ve funkci remez() (Remezův výměnný algoritmus založený na Čebyševově 
% aproximaci a realizovaný Američany Parksem a-McClelannem) je 
% implementována návrhová metoda pro FIR filtrů plně využívajícího 
% informace z tolerančního schématu.
% Funkce remez(m,F,G) navrhuje filtr m-tého řádu, jehož specifikace 
% je dána vektory F a G. F zahrnuje frekvence (mezi 0 a 1), frekvence 
% počátku a konce propustného a nepropustného pásma a G obsahuje vzorky 
% amplitudové frekvenční charakteristiky. Počet prvků vektorů F a
% G musí být stejný.
% 
% Příklad: F = [0 0.3 0.5 1] G = [1 1 0 0]
% Navrhněme dolní propust podle níže uvedených specifikací:
% Propustné pásmo: 0,0 - 0,3 Nepropustné pásmo: 0,5 - 1
% Řád filtru N: 10

clear, close all
% Vypocet optimalnich koeficientu FIR filtru 
N=10;                       % rad
f = [0 0.3 0.5 1.0];        % frekvence
m = [1 1   0   0  ];        % amplitudy
b = remez(N,f,m); 
[H,w]=freqz(b,1);
string = sprintf('Rad filtru: %d', N);
figure;
plot(f,m,w/pi,abs(H),'k','Linewidth',1.5);      % vykresleni asymptot a 
                            % frekv.charakteristiky
DPtolerance(f,H);           % vykresleni tolerancnich mezi
xlabel('\Omega/\pi'); ylabel('|H|'); title(string);

%% Příklad 4: Potlačení 50 Hz v EEG signálu s použitím FIR filtrů
% a) Navrhněme FIR filtry (paralelní zapojení horní a dolní propusti) 
% pro potlačení rušení o frekvenci 50 Hz v EEG signálu. Paralelní spojení 
% dvou dílčích filtrů umožňuje návrh realizovat. Návrh pásmové zádrže jako 
% jednoho celku není možný.
% b) Signál filtrujme a zobrazme (vzorkovací frekvence je 128 Hz).

load data.asc;                      % zaruseny signal
x=data(:,9);                        % vyber elektody
x=x-mean(x);                        % odstraneni ss slozky
x=x./max(abs(x));                   % normovani signalu
[n1,fo,mo,w] = remezord( [47 50], [1 0], [0.001 0.05], 128);
b1 = remez(94,fo,mo,w);
[n2,fo,mo,w] = remezord( [50 53], [0 1], [0.05 0.001], 128);
b2 = remez(94,fo,mo,w);

% [H1,w]=freqz(b1,1),title('Remez DP + HP'), hold on, freqz(b2,1)
[H1,w]=freqz(b1,1);
[H2,w]=freqz(b2,1);

subplot(311), plot(w/pi,20*log10(H1),'k','Linewidth',1.5), grid, axis([0 1 -60 5])
legend('Dolni propust',2), ylabel('---> dB')
subplot(312), plot(w/pi,20*log10(H2),'k','Linewidth',1.5), grid, axis([0 1 -60 5])
legend('Horni propust',2), ylabel('---> dB')
subplot(313), plot(w/pi,20*log10(H1),'k','Linewidth',1.5), grid
hold on, plot(w/pi,20*log10(H2),'k','Linewidth',1.5), hold off , axis([0 1 -60 5])
legend('Remez DP + HP',2)
xlabel('---> Normovana frekvence')
ylabel('---> dB')

figure
subplot(211), plot(x,'k'), axis tight, title('EEG signal v sumu (50 Hz)')
y1=filter(b1,1,x); y2=filter(b2,1,x);
subplot(212), plot(y1+y2,'k'), title('EEG signal po filtraci DP + HP'), axis tight
xlabel('---> n')

%% Příklad 5: Návrh FIR filtrů metodou frekvenčního vzorkování
clear
close all
% Navrhněme kombinovanou dolní propust popsanou v tabulce.
% Frekvence 0   0,15 0,25 0,45 0,5  0,75 0,85 1
% Amplituda 1   1    0,3  0,3  0,1  0,1  0    0
% vypocet FIR filtru s libovolnou amplitudovou charakteristikou
Fs=2000;
N=110;
fd=[0 0.15 0.25 0.45 0.5 0.75 0.85 1];	
Hd=[1 1 0.3 0.3 0.1 0.1 0 0];			
hn=fir2(N-1, fd, Hd);			% vypocet impulzni charakteristiky
[H, f] = freqz(hn, 1, 512, Fs);	% zobrazeni amplitudove char.
plot(f, abs(H),'k','Linewidth',1.5), grid on
xlabel('frekvence (Hz)')
ylabel('amplituda')


%% Příklad 6: Uživatelský návrh IIR filtrů
% Určeme řád filtrů splňující následující toleranční schéma pro
% =1 p f kHz, =1,4 s f kHz, = 0,5 p R dB, = 30 s R dB, = 8 vz f kHz.
% Porovnejme amplitudové frekvenční charakteristiky a z-roviny nul a pólů 
% jednotlivých typů filtrů:
% a) Butterworthův filtr
% b) Čebyševův filtr I.typu
% c) Čebyševův filtr II.typu
% d) Cauerův filtr
close all, clear
fp  = 1000;     % mezni frekvence propustneho pasma   [Hz]
fs  = 1400;     % mezni frekvence nepropustneho pasma [Hz]
Rp  = 0.5;      % zvlneni v propustnem pasmu          [dB]
Rs  = 30;       % odstup nepropustneho pasmu          [dB]
fvz = 8000;     % vzorkovaci frekvence                [Hz]
Wp  = 2*fp/fvz; % normovane frekvence
Ws  = 2*fs/fvz;

for typ=1:4 
  if typ == 1;
    [N,Wn] = buttord(Wp,Ws,Rp,Rs);    % Butterworthuv filtr
    [b,a]  = butter(N,Wn);
    FTyp = 'Butterworthuv';
  elseif typ == 2 ;  
    [N,Wn] = cheb1ord(Wp,Ws,Rp,Rs);   % Cebysevuv filtr I.typu
    [b,a]  = cheby1(N,Rp,Wn);
    FTyp = 'Cebysevuv I.typu';
  elseif typ == 3 ;  
    [N,Wn] = cheb2ord(Wp,Ws,Rp,Rs);   % Cebysevuv filtr II.typu
    [b,a]  = cheby2(N,Rs,Wn);
    FTyp = 'Cebysevuv II.typu';
  elseif typ == 4;
    [N,Wn] = ellipord(Wp,Ws,Rp,Rs);   % Caueruv filtr
    [b,a]  = ellip(N,Rp,Rs,Wn);
    FTyp = 'Caueruv';
  end;
 
  [H,W]=freqz(b,a,512);
  nadpis = sprintf('%s filtr %d. radu', FTyp, N);
  figure; 
  subplot(211), plot(W/pi,abs(H),'k','Linewidth',2); legend(nadpis); grid
  xlabel('---> \Omega/\pi'); ylabel('---> |H| [-]'); axis([0 1 0 1.2]);
  subplot(212), plot(W/pi,20*log10(abs(H)),'k','Linewidth',2); legend(nadpis); grid
  xlabel('---> \Omega/\pi'); ylabel('---> |H| [dB]'); 
  %%%
  figure; zplane2(b,a,{[0 0 0],[0 0 0]}); xlabel('---> Re'); ylabel('---> Im'); 
  title(nadpis); legend('nuly','poly')
end;

%% Příklad 7: Návrh pásmových propustí IIR

% |--------|---------------------------------------|---------|----------|
% | Typ    |           Parametry filtru            | Vstupní | Výstupní |
% |--------|---------------------------------------|---------|----------|
% | filtru | fp [Hz]  | fs [Hz] | rp [dB] | rp [dB]| signál  | signál   |
% |--------|---------------------------------------|---------|----------|
% |  DP    |   600    |   1600  |     1   |    40  |Obdélník | Flétna   |
% |  HP    |   3400   |    340  |     1   |    40  |Pila     | Trubka   |
% |  PP    |[700 2400]|[70 3400]|     1   |    40  |Obdélník | Hoboj    |
% |--------|---------------------------------------|---------|----------|


%% Příklad 8: Návrh vícepásmových IIR filtrů metodou nejmenších čtverců
% Pomocí funkce yulewalk navrhněme následující pásmové propusti:
%  Typ filtru   fp [Hz]    fs [Hz]
%  DP           600        1600
%  HP           3400       340
%  PP           [700 2400] [70 3400]
clear, close all
f=[0 600/4000 1600/4000 4000/4000]; % 'DP'
m=[1 1 0 0 ];
[dpb,dpa]=yulewalk(5,f,m);
subplot(1,3,1), plot(abs(freqz(dpb,dpa,4000)))
axis([0 4000 -0.1 1.1]), xlabel('---> f [Hz]'), legend('DP')
f=[0 340/4000 3400/4000 4000/4000]; % 'HP'
m=[0 0 1 1 ];
[hpb,hpa]=yulewalk(5,f,m);
subplot(1,3,2), plot(abs(freqz(hpb,hpa,4000)))
axis([0 4000 -0.1 1.1]), xlabel('---> f [Hz]'), legend('HP')
f=[0 70/4000 700/4000 2400/4000 3400/4000 4000/4000]; % 'PP'
m=[0 0 1 1 0 0 ];
[ppb,ppa]=yulewalk(5,f,m);
subplot(1,3,3), plot(abs(freqz(ppb,ppa,4000)))
axis([0 4000 -0.1 1.1]), xlabel('---> f [Hz]'), legend('PP')

%% Příklad 9: Barevné šumy
% Barevné šumy lze generovat filtrací bílého šumu. Je to stejné, jako když při průchodu bílého
% světla hranolem dochází k rozkladu bílého světla na barevné složky. Rovněž pojmenování
% akustických barevných šumů vychází z analogie s viditelným spektrem. Šumy s vyšším
% obsahem nízkofrekvenčních složek nesou jména barev s delší vlnovou délkou (červený,
% růžový), kdežto barvám s kratší délkou (modrý, fialový) odpovídají šumy s obsahem vyšších
% frekvencí ve spektru.
% Červený šum získáme, když bílý šum filtrujeme integrátorem, který má pokles –6dB/oktávu =
% –20 dB/dekádu). Pro tento šum se často také používá pojmenování hnědý šum, což není
% pojmenování po hnědé složce spektra, ale název vznikl zkomolením jména Brown. Viz
% Brownův pohyb, který lze modelovat autoregresním procesem 1. řádu a má výše uvedené
% frekvenční vlastnosti.
% Podobně jako růžová barva leží ve světelném spektru mezi bílou a červenou, rovněž i
% frekvenční charakteristka filtru generujícího růžový šum je méně strmá a klesá –3 dB/oktávu
% (–10 dB/dekádu). Návrh těchto filtrů bývá obtížnější; často se používají filtry se 3 póly a 3
% nulami.
% Modrý šum je generován pomocí horní propusti, jejíž frekvenční charakteristika v daném
% frekvenčním pásmu roste +3 dB/oktávu (+10 dB/dekádu).
% Fialový šum je diferencovaným bílým šumem. Frekvenční spektrum tohoto šumu se s rostoucí
% frekvencí zvětšuje o +6 dB/oktávu (+20 dB/dekádu).

% Navrhněme příslušné filtry a zobrazme jejich frekvenční charakteristiky, generujme
% barevné šumy, poslechněme si je a zobrazme jejich spektra.
% a) červený šum,
% b) růžový šum,
% c) modrý šum,
% d) fialový šum.

cas=1;                                   % doba generovani [s]
fs=44100;                                % vzorkovaci frekvence [Hz]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%bily_sum  = rand(1,fs*cas)-.5;          % bily sum = nahodne vzorky
bily_sum  = randn(1,fs*cas);             % gaussovsky

f=[0.01 0.1 1]; % cerveny sum
%f=[0.0009 0.009 0.09 0.9];
fi=pi*(1-f)/2;
H=10.^([20 0 -20]/20);
%[b1,a1]=invfreqz(H.*exp(-j*fi),f*pi,1,1);
[b1,a1]=invfreqz(H.*exp(-j*fi),f*pi,1,1);
[HH,om]=freqz(b1,a1);
figure(1), 
subplot(211), semilogx(om,20*log10(abs(HH)),'k','Linewidth',2), axis([5e-3 4 -40 40]) 
ylabel('dB'), xlabel('---> \omega'), legend('cerveny sum -20dB/dek (-6dB/okt)'), grid
subplot(223),zplane(b1,a1),
zplane2(b1,a1,{[0 0 0],[0 0 0]})
subplot(224), stem2(filter(b1,a1,[1 zeros(1,999)]),'.k')
cerveny_sum   = filter(b1,a1,bily_sum);
soundsc(cerveny_sum,fs), pause(1)

%soundsc(bily_sum,fs), pause(1)

f=[0.001 0.01 0.1 1]; % ruzovy sum
fi=pi*(1-f)/4;
H=10.^([20 10 0 -10]/20);
[b2,a2]=invfreqz(H.*exp(-j*fi),f*pi,3,3);
[HH,om]=freqz(b2,a2);
figure(2), 
subplot(211), semilogx(om,20*log10(abs(HH)),'k','Linewidth',2), axis([5e-3 4 -40 40])
ylabel('dB'), xlabel('---> \omega'), legend('ruzovy sum -10dB/dek (-3dB/okt)'), grid
subplot(223),zplane(b2,a2),zplane2(b2,a2,{[0 0 0],[0 0 0]})
subplot(224),
stem2(filter(b2,a2,[1 zeros(1,59)]),'.k')
ruzovy_sum   = filter(b2,a2,bily_sum);
soundsc(ruzovy_sum,fs), pause(1)

figure(3), 
[HH,om]=freqz(a2,b2);
subplot(211), semilogx(om,20*log10(abs(HH)),'k','Linewidth',2), axis([5e-3 4 -40 40])
ylabel('dB'), xlabel('---> \omega'), legend('modry sum +10dB/dek (+3dB/okt)'), grid
subplot(223),zplane(a2,b2),zplane2(a2,b2,{[0 0 0],[0 0 0]})
subplot(224),stem2(filter(a2,b2,[1 zeros(1,59)]),'.k')
modry_sum   = filter(a2,b2,bily_sum);
soundsc(modry_sum,fs), pause(1)

figure(4), 
[HH,om]=freqz(a1,b1);
subplot(211), semilogx(om,20*log10(abs(HH)),'k','Linewidth',2), axis([5e-3 4 -40 40])
ylabel('dB'), xlabel('---> \omega'), legend('fialovy sum +20dB/dek (+6dB/okt)'), grid
subplot(223),zplane(a1,b1),zplane2(a1,b1,{[0 0 0],[0 0 0]})
subplot(224),stem2(filter(a1,b1,[1 zeros(1,59)]),'.k')
fialovy_sum   = filter(a1,b1,bily_sum);
soundsc(fialovy_sum,fs), pause(1)