CZS cvičení 6(Matlab)

%% cvičení 6
% FIR filtr(half band)
%
%|
%|-############
%|-                        #
%|-                        #
%|-                        ##############
%|------------------------------------
%0                  1/2                     Fs/2
%                   Fs/4

%FIR1 navrhuje metodou okna

Fs = 8000; % vzorkovací frekvence
N = 30; % řád filtru
frekvence = 0:Fs/(2*N):Fs/2;



%% první část (Half band různých řádů)
figure(1);
f1 = 0.5*1/(Fs/2); % Half band - tedy polovina pásma
b = fir1(N,0.5); % vygeneruje koeficienty, N je řád filtru
freqz(b,1,1000,Fs); % a zobrazíme frekvenční charakteristiku
title('Dolní propust');
%freqz(koeficienty_b-tedy_čitatel,koeficienty_a-tedy_jmenovatel,počet_bodů_grafu,vzorkovací_frekvence)
figure(2);
%vykreslíme si charakteristiku filtru ve spektrální oblasti
plot(frekvence,b);

figure(3); % porovnáme filtry různých řádů
b = fir1(10,0.5);
% freqz vrací fázovou a amplitudovou charakteristiku [H fáz]
[frk10,ff1] = freqz(b,1,1000,Fs);
b = fir1(100,0.5);
[frk100 ff2] = freqz(b,1,1000,Fs);
b = fir1(200,0.5);
[frk200 ff3] = freqz(b,1,1000,Fs);
b = fir1(500,0.5);
[frk500 ff1] = freqz(b,1,1000,Fs);
hold on;
frek = 1:4:Fs/2;
% semilogy - y-lonová osa bude logaritmická
% zobrazujeme absolutní hodnotu amplitudy
semilogy(frek,abs(frk10));
semilogy(frek,abs(frk100),'r');
semilogy(frek,abs(frk200),'g');
semilogy(frek,abs(frk500),'y');
% legenda pro graf
legend('N = 10','N = 100','N = 200', 'N = 500');

figure(4);
% tady budeme dělat horní propust, jedná se o převrácený jev proti DP
b = fir1(N,0.5,'high');
freqz(b,1,1000,Fs); % a zobrazíme frekvenční charakteristiku
title('Horní propust');
%% druhá část(filtrace pomocí Half band různých řádů)
figure(5);
data = loadbin('s0001.bin');
subplot(311);
plot(data);
%soundsc(data);
b10 = fir1(10,0.5);
b500 = fir1(500,0.5);
[frk500 ff1] = freqz(b,1,1000,Fs);
subplot(312);
title('Filtrováno N = 10');
dataout = filter(b10,1,data);
%soundsc(dataout);
plot(dataout);
subplot(313);
dataout = filter(b10,1,data);
plot(dataout,'r');
%soundsc(dataout);
title('Filtrováno N = 500');
%% Návrh pásmové propusti
Fs = 16e3;
Fd = 300;
Fh = 3400;
N = 500;

Fd_norm = Fd/(Fs/2);
Fh_norm = Fh/(Fs/2);
b = fir1(N,[Fd_norm Fh_norm],'DC-0'); % tohle bude filtr typu pásmová propust
figure(6);
freqz(b,1,1000,Fs);
data = loadbin('SA106S06.CS0');%fc44be305265.ils_a');%SA106S06.CS0');%s0001.bin');
soundsc(data,Fs);
soundsc(filter(b,1,data),Fs);
% Pozor zde se nezabýváme stabilitou filtru!!!
http://noel.feld.cvut.cz/vyu/a2m99czs/cv_Filtry.html