57409dbc4b9677388c289bab48f72172054a5409
[scilab.git] / scilab / modules / signal_processing / help / en_US / howto / DesignEllipticFilter.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <refentry version="5.0-subset Scilab" xml:id="DesignEllipticFilter"
3           xmlns="http://docbook.org/ns/docbook"
4           xmlns:xlink="http://www.w3.org/1999/xlink"
5           xmlns:xi="http://www.w3.org/2001/XInclude"
6           xmlns:svg="http://www.w3.org/2000/svg"
7           xmlns:mml="http://www.w3.org/1998/Math/MathML"
8           xmlns:html="http://www.w3.org/1999/xhtml"
9           xmlns:db="http://docbook.org/ns/docbook">
10   <refnamediv>
11     <refname>How to design an elliptic filter</refname>
12
13     <refpurpose>How to design an elliptic filter (analog and
14     digital)</refpurpose>
15   </refnamediv>
16
17   <refsection>
18     <title>Description</title>
19
20     <para>The goal is to design a simple analog and digital elliptic
21     filter.</para>
22   </refsection>
23
24   <refsection>
25     <title>Designing an analog elliptic filter</title>
26
27     <para>There are several possibilities to design an elliptic lowpass
28     filter. We can use <link linkend="analpf">analpf</link> or <link
29     linkend="zpell">zpell</link>. We will use zpell to produce the poles and
30     zeros of the filter. Once we have got these poles and zeros, we will have
31     to translate this representation into a <link
32     linkend="syslin">syslin</link> one.</para>
33
34     <para>And then, the filter can be represented in bode plot.</para>
35
36     <programlisting> 
37 // analog elliptic (Bessel), order 2, cutoff 1 Hz
38 Epsilon = 3;  // ripple of filter in pass band (0&lt;epsilon&lt;1)
39 A       = 60; // attenuation of filter in stop band (A&gt;1)
40 OmegaC  = 10; // pass band cut-off frequency in Hertz
41 OmegaR  = 50; // stop band cut-off frequency in Hertz
42
43 // Generate the filter
44 [_zeros,pols,gain] = zpell(3,60,10,50);
45
46 // Generate the equivalent linear system of the filter
47 num   = gain * real(poly(_zeros,'s'));;
48 den   = real(poly(pols,'s'));
49 elatf = syslin('c',num,den);
50
51 // Plot the resulting filter
52 bode(elatf,0.01,100);
53 title('Analog Elliptic filter');
54  </programlisting>
55
56     <para>Bode plot is only suited for analog filters.</para>
57
58     <mediaobject>
59       <imageobject>
60         <imagedata align="center"
61                    fileref="../../images/analog_elliptic_filter.png" />
62       </imageobject>
63     </mediaobject>
64
65     <para>If you want to design a highpass, bandpass or bandstop filter, you
66     can first design a lowpass and then transfrom this lowpass filter using
67     the <link linkend="trans">trans</link> function.</para>
68   </refsection>
69
70   <refsection>
71     <title>Designing a digital elliptic filter</title>
72
73     <para>Now, let's focus on how to produce a digital lowpass elliptic
74     filter.</para>
75
76     <para>We can produce two kinds of digital filters:</para>
77
78     <itemizedlist>
79       <listitem>
80         <para>an IIR (Infinite Impulse Response).</para>
81
82         <para>To compute such a filter, we can use the following
83         functions:</para>
84
85         <itemizedlist>
86           <listitem>
87             <para><link linkend="iir">iir</link></para>
88           </listitem>
89
90           <listitem>
91             <para><link linkend="eqiir">eqiir</link></para>
92           </listitem>
93         </itemizedlist>
94       </listitem>
95
96       <listitem>
97         <para>a FIR (Finite Impulse Response).</para>
98
99         <para>To compute such a filter, we can use the following
100         functions:</para>
101
102         <itemizedlist>
103           <listitem>
104             <para><link linkend="eqfir">eqfir</link></para>
105           </listitem>
106
107           <listitem>
108             <para><link linkend="ffilt">ffilt</link></para>
109           </listitem>
110
111           <listitem>
112             <para><link linkend="wfir">wfir</link></para>
113           </listitem>
114
115           <listitem>
116             <para><link linkend="fsfirlin">fsfirlin</link></para>
117           </listitem>
118         </itemizedlist>
119       </listitem>
120     </itemizedlist>
121
122     <para>For our demonstration, we will use the <link
123     linkend="iir">iir</link> function.</para>
124
125     <programlisting role="example"> 
126 Order   = 2; // The order of the filter
127 Fs      = 1000; // The sampling frequency
128 Fcutoff = 40;   // The cutoff frequency
129
130 // We design a low pass elliptic filter
131 hz = iir(Order,'lp','ellip',[Fcutoff/Fs/2 0],[0.1 0.1]);
132
133 // We compute the frequency response of the filter
134 [frq,repf]=repfreq(hz,0:0.001:0.5);
135 [db_repf, phi_repf] = dbphi(repf);
136
137 // And plot the bode like representation of the digital filter
138 subplot(2,1,1);
139 plot2d(Fs*frq,db_repf);
140 xtitle('Obtained Frequency Response (Magnitude)');
141 subplot(2,1,2);
142 plot2d(Fs*frq,phi_repf);
143 xtitle('Obtained Frequency Response (Phase in degree)');
144  </programlisting>
145
146     <para>Here is the representation of the digital elliptic filter.</para>
147
148     <mediaobject>
149       <imageobject>
150         <imagedata align="center"
151                    fileref="../../images/digital_elliptic_filter.png" />
152       </imageobject>
153     </mediaobject>
154
155     <para>To represent the filter in phase and magnitude, we need first to
156     convert the discrete impulse response into magnitude and phase using the
157     <link linkend="dbphi">dbphi</link> function. This convertion is done using
158     a set of normalized frequencies.</para>
159   </refsection>
160
161   <refsection>
162     <title>Filtering a signal using the digital filter</title>
163
164     <para>Designing a filter is a first step. Once done, this filter will be
165     used to transform a signal. To get rid of some noise for example.</para>
166
167     <para>In the following examples, we will filter a gaussian noise.</para>
168
169     <programlisting>
170 rand('normal');
171 Input = rand(1,1000); // Produce a random gaussian noise
172 t     = 1:1000;
173
174 sl= tf2ss(hz); // From transfert function to syslin representation
175 y = flts(Input,sl); // Filter the signal
176
177 subplot(2,1,1);
178 plot(t,Input);
179 xtitle('The gaussian noise','t','y');
180 subplot(2,1,2);
181 plot(t,y);
182 xtitle('The filtered gaussian noise','t','y');
183  </programlisting>
184
185     <para>Here is the representation of the signal before and after
186     filtering.</para>
187
188     <mediaobject>
189       <imageobject>
190         <imagedata fileref="../../images/digital_filtered_noise.png" />
191       </imageobject>
192     </mediaobject>
193
194     <para>As we can see in the result, the high frequencies of the noise have
195     been removed and it remains only the low frequencies. The signal is still
196     noisy, but it contains mainly low frequencies.</para>
197   </refsection>
198
199   <refsection>
200     <title>Filtering a signal using the analog filter</title>
201
202     <para>To filter a signal using an analog filter, we have two
203     strategies:</para>
204
205     <itemizedlist>
206       <listitem>
207         <para>transform the analog filter into a discrete one using the <link
208         linkend="dscr">dscr</link> function</para>
209       </listitem>
210
211       <listitem>
212         <para>apply the <link linkend="csim">csim</link> function to filter
213         the signal</para>
214       </listitem>
215     </itemizedlist>
216
217     <para>First, we try using the <link linkend="dscr">dscr</link> + <link
218     linkend="flts">flts</link> functions.</para>
219
220     <programlisting> 
221 rand('normal');
222 Input = rand(1,1000); // Produce a random gaussian noise
223 n     = 1:1000; // The sample index
224
225 eldtf = dscr(elatf,1/100); // Discretization of the linear filter
226 y = flts(Input,eldtf); // Filter the signal
227
228 subplot(2,1,1);
229 plot(n,Input);
230 xtitle('The gaussian noise','n','y');
231 subplot(2,1,2);
232 plot(n,y);
233 xtitle('The filtered gaussian noise','n','y');
234  </programlisting>
235
236     <para>Here is the representation of the signal before and after filtering
237     using the <link linkend="dscr">dscr</link> + <link
238     linkend="flts">flts</link> approach.</para>
239
240     <mediaobject>
241       <imageobject>
242         <imagedata fileref="../../images/analog_filtered_noise.png" />
243       </imageobject>
244     </mediaobject>
245
246     <para>Next, we use the <link linkend="csim">csim</link> function.</para>
247
248     <programlisting> 
249 rand('normal');
250 Input = rand(1,1000); // Produce a random gaussian noise
251 t     = 1:1000;
252 t     = t*0.01; // Convert sample index into time steps
253
254 y = csim(Input,t,elatf); // Filter the signal
255
256 subplot(2,1,1);
257 plot(t,Input);
258 xtitle('The gaussian noise','t','y');
259 subplot(2,1,2);
260 plot(t,y);
261 xtitle('The filtered gaussian noise','t','y');
262  </programlisting>
263
264     <para>Here is the representation of the signal before and after filtering
265     using the <link linkend="csim">csim</link> approach.</para>
266
267     <mediaobject>
268       <imageobject>
269         <imagedata fileref="../../images/analog_filtered_noise_csim.png" />
270       </imageobject>
271     </mediaobject>
272
273     <para>The main difference between the <link linkend="dscr">dscr</link> +
274     <link linkend="flts">flts</link> approach and the <link
275     linkend="csim">csim</link> approach: the <link linkend="dscr">dscr</link>
276     + <link linkend="flts">flts</link> uses samples whereas the <link
277     linkend="csim">csim</link> functions uses time steps.</para>
278   </refsection>
279
280   <refsection>
281     <title>See Also</title>
282
283     <simplelist type="inline">
284       <member><link linkend="bode">bode</link></member>
285
286       <member><link linkend="iir">iir</link></member>
287
288       <member><link linkend="poly">poly</link></member>
289
290       <member><link linkend="syslin">syslin</link></member>
291
292       <member><link linkend="zpell">zpell</link></member>
293
294       <member><link linkend="flts">flts</link></member>
295
296       <member><link linkend="tf2ss">tf2ss</link></member>
297
298       <member><link linkend="dscr">dscr</link></member>
299
300       <member><link linkend="csim">csim</link></member>
301
302       <member><link linkend="trans">trans</link></member>
303
304       <member><link linkend="analpf">analpf</link></member>
305     </simplelist>
306   </refsection>
307 </refentry>