add pictures & MATLAB

infotro/Eln/labs/06-filters/README.md

@@ -11,21 +11,117 @@ We will use MATLAB to simulate and implement these filters on a digital signal p
 
 We must design Butterworth low-pass filter that cuts at 1kHz as shown with the following equation :
 
-$$
-F(s) = \frac{1}{\frac{s^{2}}{\omega_c^{2}} + s\cdot\frac{1.414}{\omega_c} + 1}
-$$
+$$F(s) = \frac{1}{\frac{s^{2}}{\omega_c^{2}} + s\cdot\frac{1.414}{\omega_c} + 1} $$
 
 In according to the sample theory, we have to choose at least 2 times the maximum frequency of the signal.
-In this case, the maximum frequency is 1 kHz, so we have to choose a sample frequency of 2 kHz.
+In this case, the maximum frequency is 1 kHz, so we have to choose a sample frequency at least at 2 kHz.
+We take 10 time $f_c$ for the sample frequency, so we choose 20 kHz.
 
 ## Transfer function
 For calculate the transfer function, we have to calculate the poles and zeros of the filter.
 
-$$
-p_{di} = e^{p_{ai} \cdot T} \qquad z_{di} = e^{z_{ai} \cdot T}
-$$
+$$p_{di} = e^{p_{ai} \cdot T} \qquad z_{di} = e^{z_{ai} \cdot T} $$
 
-- $k_a$, $p_{ai}$ and $z_{ai}$ are the gain, poles and zeros of the analog filter that the representation `zpk`.
+- $k_a$, $p_{ai}$ and $z_{ai}$ are the gain, poles and zeros of the analog filter that the representation *zpk*.
 Be careful, $k_a$ is not the DC gain.
 - $k_d$, $p_{di}$ and $z_{di}$ are the gain, poles and zeros of the digital filter.
 
+For a low pass filter without zeros, the "digital" gain is :
+$$k_{d} = k_{a}\cdot \prod_{i} \frac{1-e^{P_{ai}\cdot T}}{P_{ai}}$$
+
+----------
+MatLab:
+% tfchk         % Check for proper transfer function
+% zp2tf         % Convert zero-pole-gain filter parameters to transfer function form
+% filter        % 1-D digital filter
+Signal Processing Toolbox :
+% eqtflength    % Equalize lengths of transfer function numerator and denominator
+% residuez      % Z-transform partial-fraction expansion
+% sos2tf        % Convert digital filter second-order section data to transfer function form
+% sos2zp        % Convert digital filter second-order section parameters to zero-pole-gain form
+% tf2sos        % Convert digital filter transfer function data to second-order sections form
+% tf2zpk        % Convert transfer function filter parameters to zero-pole-gain form
+% zp2sos        % Convert zero-pole-gain filter parameters to second-order sections form
+% freqz         % Frequency response of digital filter
+% freqs         % Frequency response of analog filters
+% impz          % Impulse response of digital filter
+% stepz         % Step response of digital filter
+% zplane        % Zero-pole plot for discrete-time systems
+% bilinear      % Bilinear transformation method for analog-to-digital filter conversion
+% ellipord      % Minimum order for elliptic filters
+% ellip         % Elliptic filter design
+% chirp         % Swept-frequency cosine
+Control Toolbox :
+% bode          % Bode plot of frequency response, or magnitude and phase data
+% zpk           % Zero-pole-gain model
+% zpkdata       % Access zero-pole-gain data
+% tf            % Transfer function model
+% tfdata        % Access transfer function data
+% lsim          % Plot simulated time response of dynamic system to arbitrary inputs; simulated response data
+
+
+fc = 1000;
+wc = 2*pi*fc;
+fs = 10*fc;
+ts = 1/fs;
+
+num = [0 0 1];
+den = [1/(wc*wc) 1.414/wc 1];
+sysa = tf(num, den)
+[za,pa,ka] = tf2zp(num, den);
+%fvtool(num,den,'polezero')
+
+
+kd = ka .* ((1-exp(pa(1).*ts))./pa(1)) .* ((1-exp(pa(2).*ts))./pa(2));
+kd = real(kd);
+pd = (1+pa.*(ts./2)) ./ (1-pa.*(ts./2));
+zd = (1+za.*(ts/2)) ./ (1-za.*(ts/2));
+sysd = zpk(zd, pd, kd, ts)
+[zdb, pdb, kdb] = bilinear(za, pa, ka, fs);
+sysdb = zpk(zdb, pdb, kdb, ts)
+
+bode(sysa, sysd, sysdb);
+title('Bode Diagram', 'Color', '#808080');
+xlabel('Frequency (rad/s)', 'Color', '#808080');
+ylabel('Magnitude (dB)', 'Color', '#808080');
+ylabel('Phase (deg)', 'Color', '#808080');
+legend('analog','digital','bilinear', 'Color', 'white')
+
+
+![bode_design_filter.svg](bode_design_filter.svg)
+
+-----------
+clear;
+close all;
+clc;
+
+H0 = 0;      % dB
+H1 = -3;     % dB
+H2 = -40;    % dB
+fn.p = 7000;   % Hz
+fn.s = 9000   % Hz
+fs = 48000;  % Hz
+
+fa.p = (fs/pi) * tan((pi*fn.p)/fs);
+fa.s = (fs/pi) * tan((pi*fn.s)/fs)
+
+wp = 2*pi*fa.p;
+ws = 2*pi*fa.s;
+ts = 1/fs;
+
+[n, wp] = ellipord(wp, ws, 3, 40, 's')
+[z, p, k] = ellip(n, 3, 40, wp, 's');
+sys = zpk(z, p, k);
+[zd, pd, kd] = bilinear(z, p, k, fs);
+sysd = zpk(zd, pd, kd, ts)
+bode(sys, sysd)
+title('Bode Diagram', 'Color', '#808080');
+xlabel('Frequency (rad/s)', 'Color', '#808080');
+ylabel('Magnitude (dB)', 'Color', '#808080');
+ylabel('Phase (deg)', 'Color', '#808080');
+legend('digital','bilinear', 'Color', 'white')
+sos = zp2sos(z,p,k);
+
+![iir_filter.svg](iir_filter.svg)
+
+