vzoreèek pro výpoèet kv. rovnice poskytl 
Tomá¹ Studený xstude02@stud.fit.vutbr.cz

% aux pro diskretni systemy 

%%%%%%%%%%%%%%%%%%%%%%%% H pro filtr 3 2 1
n = 0:2; Fs = 8000; gj0 = find(n==0); T = 1/Fs; 
x = [3 2 1]; 
subplot(311); stem(n,x); grid; title('h')
f = -2*Fs:100:2*Fs; 
fn = f / Fs; 
Om = 2*pi*f;
om = Om / Fs; 
% dtft
X = zeros(size(om)); 
for ii=1:length(om)
  ejon = exp(j*om(ii).*n); 
  X(ii) = sum(x .* ejon); 
end
subplot(312); plot (om,abs(X)); grid; axis tight; title('module');
subplot(313); plot (om,angle(X)); grid; axis tight; title('phase');
print -depsc FIG/h.eps

subplot(221); plot (om,abs(X)); grid; axis tight; title('normalized omega');
subplot(222); plot (Om,abs(X)); grid; axis tight; title('omega');
subplot(223); plot (fn,abs(X)); grid; axis tight; title('normalized f');
subplot(224); plot (f,abs(X)); grid; axis tight;  title('f'); 
print -depsc FIG/hfs.eps


%%%%%%%%%%%%%%%%% priklady filtry %%%%%%%%%%%%%%%%%%%%%
H=freqz([1 0.5],[1],256); om=(0:255)/256 * pi;
subplot(211); plot(om,abs(H)); grid
subplot(212); plot(om,angle(H)); grid
print -depsc FIG/h_1+0.5.eps

zplane([1 0.5],[1]); 
print -depsc FIG/np_1+0.5.eps
%%% iir
n = 0:30;
y= filter([1],[1 0.5],[1 zeros(1,30)]); stem(n,y); 
print -depsc FIG/imp_rek_1+0.5.eps
H=freqz([1],[1 0.5],256); om=(0:255)/256 * pi;
subplot(211); plot(om,abs(H)); grid
subplot(212); plot(om,angle(H)); grid
print -depsc FIG/h_rek_1+0.5.eps

zplane([1],[1 0.5]); 
print -depsc FIG/np_rek_1+0.5.eps

%%%%%%%%%%%%%%%%%%%%% priklady filtru 2014-15 %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nejaky signal - -123 min ? ne, to je v labinach ... liska. 
vezmu si Smyfoniettu ze zacatku ... 
    
% sinf_1ch.wav - 1 channel, 48k. this is OK. 

Fs = 48000; Fs2 = Fs / 2; 
s = wavread ('sinf_1ch.wav'); 
sum = 0.1*randn(1,2*Fs); wavwrite(sum,Fs,16,'sum.wav'); 
f = (0:511) / 512 * Fs2; 
%%%%%%%%%%%%%%%%%% low pass, till 3400 ... 
% vyrobime filtr
Wp = 3400 / Fs2;  % pass-band
Ws = 4000 / Fs2;  % stop-band
Rp = 3;               % zvlneni pass
Rs = 40;              % potlaceni stop

[N, Wn] = ellipord(Wp, Ws, Rp, Rs)
[B,A] = ellip(N,Rp,Rs,Wp); 

figure(1); 
H = freqz (B,A,512,Fs); 
subplot (211); plot (f, abs(H)); xlabel ('f [Hz]'); ylabel ('|H|'); grid
subplot (212); plot (f, angle(H)); xlabel ('f [Hz]'); ylabel ('arg (H)'); grid
  print -depsc lp.eps
figure(2); zplane (B,A); grid; 
y = filter(B,A,s);   print -depsc lp_z.eps
wavwrite(y,Fs,16,'lp.wav'); 
ys = filter(B,A,sum); wavwrite(ys,Fs,16,'lp_sum.wav'); 

%%%%%%%%%%%%%%%%%% high  pass, from  5000 ... 
% vyrobime filtr
Wp = 5000 / Fs2;  % pass-band
Ws = 4500 / Fs2;  % stop-band
Rp = 3;               % zvlneni pass
Rs = 40;              % potlaceni stop

[N, Wn] = ellipord(Wp, Ws, Rp, Rs)
[B,A] = ellip(N,Rp,Rs,Wp,'high'); 

figure(1); 
H = freqz (B,A,512,Fs); 
subplot (211); plot (f, abs(H)); xlabel ('f [Hz]'); ylabel ('|H|'); grid
subplot (212); plot (f, angle(H)); xlabel ('f [Hz]'); ylabel ('arg (H)'); grid
  print -depsc hp.eps
figure(2); zplane (B,A); grid; 
y = filter(B,A,s);   print -depsc hp_z.eps
wavwrite(y,Fs,16,'hp.wav'); 
ys = filter(B,A,sum); wavwrite(ys,Fs,16,'hp_sum.wav'); 

%%%%%%%%%%%%%%%%%% band  pass, 1000 az 2000
% vyrobime filtr
Wp = [1000 2000] / Fs2;  % pass-band
Ws = [500 2500] / Fs2;  % stop-band
Rp = 3;               % zvlneni pass
Rs = 40;              % potlaceni stop

[N, Wn] = ellipord(Wp, Ws, Rp, Rs)
[B,A] = ellip(N,Rp,Rs,Wn); 

figure(1); 
H = freqz (B,A,512,Fs); 
subplot (211); plot (f, abs(H)); xlabel ('f [Hz]'); ylabel ('|H|'); grid
subplot (212); plot (f, angle(H)); xlabel ('f [Hz]'); ylabel ('arg (H)'); grid
  print -depsc bp.eps
figure(2); zplane (B,A); grid; 
y = filter(B,A,s);   print -depsc bp_z.eps
wavwrite(y,Fs,16,'bp.wav'); 
ys = filter(B,A,sum); wavwrite(ys,Fs,16,'bp_sum.wav'); 

%%%%%%%%%%%%%%%%%% band  stop  440 - 880 Hz 
Wp = [400 950] / Fs2;  % pass-band
Ws = [440 880] / Fs2;  % stop-band
Rp = 3;               % zvlneni pass
Rs = 40;              % potlaceni stop

[N, Wn] = ellipord(Wp, Ws, Rp, Rs)
[B,A] = ellip(N,Rp,Rs,Wn,'stop'); 

figure(1); 
H = freqz (B,A,512,Fs); 
subplot (211); plot (f, abs(H)); xlabel ('f [Hz]'); ylabel ('|H|'); grid
subplot (212); plot (f, angle(H)); xlabel ('f [Hz]'); ylabel ('arg (H)'); grid
  print -depsc bs.eps
figure(2); zplane (B,A); grid; 
y = filter(B,A,s);   print -depsc bs_z.eps
wavwrite(y,Fs,16,'bs.wav'); 
ys = filter(B,A,sum); wavwrite(ys,Fs,16,'bs_sum.wav'); 

