1.  
2. % Wprowadzenie do filtracji cyfrowej - filtry rekurencyjne
3. % filtry o nieskończonej odpowiedzi impulsowej
4.  
5. clear all; close all;
6.  
7. % Sygnał #################################################################
8. if(0)
9.    Nx = 1000; fpr1=8000;
10.    x = zeros(1,Nx);
11.    x(1) = 1;
12. end    
13. if(0)
14.    fpr1 = 8000;
15.    Nx = 3*fpr1;
16.    dt = 1/fpr1; t=dt*(0:Nx-1);
17.    x = sin( 2*pi*200*t) + sin( 2*pi*2000*t);
18. end  
19. if(1)
20.    [ x, fpr1 ] = audioread('mowa.wav');
21.  % [ x, fpr1 ] = audioread('track1.wav');
22.    disp('Rozmiar audio'); size(x),
23.    x = x(:,1);
24.    Nx = length(x);
25.    pause
26. end
27. %figure; plot(x); pause
28. figure; stem(x); pause
29. figure;
30. f = fpr1/Nx * [ -Nx/2 : Nx/2-1 ];
31. %plot( f,20*log10( fftshift(abs(fft(x)/Nx)) )); pause
32. soundsc( x, fpr1 );
33.  
34. % Wagi filtra #############################################################
35.  
36. %b = [ 2 3 4 ];
37. %a = [ 1 2 ];
38.  
39. if(0)
40.   z = [ 1 1 ] .* exp( j*2*pi*[ 2000 3000 ]/fpr1 );
41.   p = [ 0.99 0.99 ] .* exp( j*2*pi*[ 100 200 ]/fpr1 );
42.   z = [ z conj(z) ];
43.   p = [ p conj(p) ];
44.  
45.   b = poly( z );
46.   a = poly( p );
47. end
48. if(1)
49.    N = 8;
50.    f0 = 1000;
51.    
52.    pasmo = [1000 2000]
53.    
54.   [b,a] = butter( N, pasmo/(fpr1/2) );
55.   [b1,a1] = cheby1( N, 3, pasmo/(fpr1/2) );
56.   [b2,a2] = cheby2( N, 100, pasmo/(fpr1/2) );
57.   [b3,a3] = ellip( N, 3, 100, pasmo/(fpr1/2) );
58.  
59.   z = roots( b );
60.   p = roots( a );
61.   z1 = roots( b1 );
62.   p1 = roots( a1 );
63.   z2 = roots( b2 );
64.   p2 = roots( a2 );
65.   z3 = roots( b3 );
66.   p3 = roots( a3 );
67.  
68. end
69.  
70. figure;
71. alfa = 0 : 0.01 : 2*pi; c=cos(alfa); s=sin(alfa);
72. %plot(c,s,'k-',real(z),imag(z),'bo',real(p),imag(p),'r*'); grid; title('Z & P'); pause
73.  
74.  
75.  
76. f = 0 : 1 : fpr1/2;
77. zz = exp( -j*2*pi*f/fpr1);
78. H = polyval( b, zz ) ./polyval( a, zz );
79. figure;
80. plot( f, 20*log10( abs( H ) ),'r' ); hold on;
81. H1 = polyval( b1, zz ) ./polyval( a1, zz );
82. plot( f, 20*log10( abs( H1 ) ),'g' ); hold on;
83. H2 = polyval( b2, zz ) ./polyval( a2, zz );
84. plot( f, 20*log10( abs( H2 ) ),'b' ); hold on;
85. H3 = polyval( b3, zz ) ./polyval( a3, zz );
86. plot( f, 20*log10( abs( H3 ) ),'y' ); xlabel('f [Hz]'); title('|H(f)| [dB]'); grid;
87. legend('butter','cheby1','cheby2','ellip'); hold on; pause
88.  
89.  
90.  
91.  
92. figure;
93. bode(a,b,'r'); hold on;
94. bode(a1,b1,'g'); hold on;
95. bode(a2,b2,'b'); hold on;
96. bode(a3,b3,'y');
97. grid; legend('butter','cheby1','cheby2','ellip'); hold on; pause
98.  
99.  
100. figure;
101. plot(f,unwrap(angle(H)),'r'); hold on;
102. plot(f,unwrap(angle(H1)),'g'); hold on;
103. plot(f,unwrap(angle(H2)),'b'); hold on;
104. plot(f,unwrap(angle(H3)),'y'); hold on;
105. xlabel('[Hz]'),ylabel('[rad/s]')
106. grid; legend('butter','cheby1','cheby2','ellip'); hold on; pause
107.  
108.  
109.  
110. % Filtracja ###############################################################
111. if(1)
112.   M = length(b);
113.   bx = zeros(1,M);
114.   N = length(a);
115.   a = a(2:N);
116.   N = N-1;
117.   by = zeros(1,N);
118.  
119.   M1 = length(b1);
120.   bx1 = zeros(1,M1);
121.   N1 = length(a1);
122.   a1 = a1(2:N1);
123.   N1 = N1-1;
124.   by1 = zeros(1,N1);
125.  
126.   M2 = length(b2);
127.   bx2 = zeros(1,M2);
128.   N2 = length(a2);
129.   a2 = a2(2:N2);
130.   N2 = N2-1;
131.   by2 = zeros(1,N2);
132.  
133.   M3 = length(b3);
134.   bx3 = zeros(1,M3);
135.   N3 = length(a3);
136.   a3 = a3(2:N3);
137.   N3 = N3-1;
138.   by3 = zeros(1,N3);
139.  
140.  
141.   for n = 1 : Nx
142.       bx = [ x(n) bx(1:M-1) ];
143.       y(n) = sum( bx .* b ) - sum( by .* a );
144.       by = [ y(n) by(1:N-1) ];
145.      
146.       bx1 = [ x(n) bx1(1:M1-1) ];
147.       y1(n) = sum( bx1 .* b1 ) - sum( by1 .* a1 );
148.       by1 = [ y1(n) by1(1:N1-1) ];
149.      
150.       bx2 = [ x(n) bx2(1:M2-1) ];
151.       y2(n) = sum( bx2 .* b2 ) - sum( by2 .* a2 );
152.       by2 = [ y2(n) by2(1:N2-1) ];
153.      
154.       bx3 = [ x(n) bx3(1:M3-1) ];
155.       y3(n) = sum( bx3 .* b3 ) - sum( by3 .* a3 );
156.       by3 = [ y3(n) by3(1:N3-1) ];
157.   end
158. else
159.     y = filter(b,a,x);
160.     y1 = filter(b1,a1,x);
161.     y2 = filter(b2,a2,x);
162.     y3 = filter(b3,a3,x);
163.  
164. end
165.  
166. figure
167. subplot(211); plot(x(1:400)); title('WE');
168. subplot(212);
169. plot(y(1:400),'r'); hold on;
170. plot(y1(1:400),'g'); hold on;
171. plot(y2(1:400),'b'); hold on;
172. plot(y3(1:400),'y'); hold on;
173. title('WY');legend('butter','cheby1','cheby2','ellip');
174. pause
175.  
176.  
177. figure
178. subplot(211); plot(f,20*log10(fftshift(abs(fft(x)/length(x))))); title('FFT WE');
179. subplot(212); plot(f,20*log10(fftshift(abs(fft(y)/length(y))))); title('FFT WY'); pause
180. soundsc( y, fpr1 );
181.