1 //
2 // iirdes_analog_example.c
3 //
4 // Tests infinite impulse reponse (IIR) analog filter design.
5 // While this example seems purely academic as IIR filters
6 // used in liquid are all digital, it is important to realize
7 // that they are all derived from their analog counterparts.
8 // This example serves to check the response of the analog
9 // filters to ensure they are correct. The results of design
10 // are written to a file.
11 // SEE ALSO: iirdes_example.c
12 // iir_filter_crcf_example.c
13 //
14
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <string.h>
18 #include <getopt.h>
19 #include <math.h>
20 #include <complex.h>
21
22 #include "liquid.h"
23
24 #define OUTPUT_FILENAME "iirdes_analog_example.m"
25
26 // print usage/help message
usage()27 void usage()
28 {
29 printf("iirdes_analog_example -- infinite impulse response filter design\n");
30 printf("options (default values in []):\n");
31 printf(" u/h : print usage/help\n");
32 printf(" t : filter type: [butter], cheby1, cheby2, ellip, bessel\n");
33 // printf(" b : filter transformation: [LP], HP, BP, BS\n");
34 printf(" n : filter order, n > 0 [5]\n");
35 printf(" r : passband ripple in dB (cheby1, ellip), r > 0 [3.0]\n");
36 printf(" s : stopband attenuation in dB (cheby2, ellip), s > 0 [60.0]\n");
37 printf(" f : angular passband cut-off frequency, f > 0 [1.0]\n");
38 // printf(" c : center frequency (BP, BS cases), 0 < c < 0.5 [0.25]\n");
39 // printf(" o : format [sos], tf\n");
40 // printf(" sos : second-order sections form\n");
41 // printf(" tf : regular transfer function form (potentially\n");
42 // printf(" unstable for large orders\n");
43 }
44
main(int argc,char * argv[])45 int main(int argc, char*argv[]) {
46 // options
47 unsigned int order=3; // filter order
48 float wc=1.0f; // angular cutoff frequency
49 float Ap=3.0f; // pass-band Ap
50 float As=60.0f; // stop-band attenuation
51
52 // filter type
53 liquid_iirdes_filtertype ftype = LIQUID_IIRDES_BUTTER;
54
55 int dopt;
56 while ((dopt = getopt(argc,argv,"uht:n:r:s:f:")) != EOF) {
57 switch (dopt) {
58 case 'u':
59 case 'h':
60 usage();
61 return 0;
62 case 't':
63 if (strcmp(optarg,"butter")==0) {
64 ftype = LIQUID_IIRDES_BUTTER;
65 } else if (strcmp(optarg,"cheby1")==0) {
66 ftype = LIQUID_IIRDES_CHEBY1;
67 } else if (strcmp(optarg,"cheby2")==0) {
68 ftype = LIQUID_IIRDES_CHEBY2;
69 } else if (strcmp(optarg,"ellip")==0) {
70 ftype = LIQUID_IIRDES_ELLIP;
71 } else if (strcmp(optarg,"bessel")==0) {
72 ftype = LIQUID_IIRDES_BESSEL;
73 } else {
74 fprintf(stderr,"error: %s, unknown filter type \"%s\"\n", argv[0], optarg);
75 usage();
76 exit(1);
77 }
78 break;
79 case 'n': order = atoi(optarg); break;
80 case 'r': Ap = atof(optarg); break;
81 case 's': As = atof(optarg); break;
82 case 'f': wc = atof(optarg); break;
83 default:
84 exit(1);
85 }
86 }
87
88 // validate input
89 if (wc <= 0) {
90 fprintf(stderr,"error: %s, cutoff frequency out of range\n", argv[0]);
91 usage();
92 exit(1);
93 } else if (Ap <= 0) {
94 fprintf(stderr,"error: %s, pass-band ripple out of range\n", argv[0]);
95 usage();
96 exit(1);
97 } else if (As <= 0) {
98 fprintf(stderr,"error: %s, stop-band ripple out of range\n", argv[0]);
99 usage();
100 exit(1);
101 }
102
103 // number of analaog poles/zeros
104 unsigned int npa = order;
105 unsigned int nza = 0;
106
107 // complex analog zeros, poles, gain
108 float complex za[order];
109 float complex pa[order];
110 float complex ka;
111
112 unsigned int i;
113
114 unsigned int r = order % 2;
115 unsigned int L = (order-r)/2;
116
117 float Gp, Gs;
118 float ep = sqrtf( powf(10.0f, Ap / 10.0f) - 1.0f );
119 float es = powf(10.0f, -As / 20.0f);
120
121 switch (ftype) {
122 case LIQUID_IIRDES_BUTTER:
123 printf("Butterworth filter design:\n");
124 nza = 0;
125 butter_azpkf(order,za,pa,&ka);
126 break;
127 case LIQUID_IIRDES_CHEBY1:
128 printf("Cheby-I filter design:\n");
129 nza = 0;
130 cheby1_azpkf(order,ep,za,pa,&ka);
131 break;
132 case LIQUID_IIRDES_CHEBY2:
133 printf("Cheby-II filter design:\n");
134 nza = 2*L;
135 float epsilon = powf(10.0f, -As/20.0f);
136 cheby2_azpkf(order,epsilon,za,pa,&ka);
137 break;
138 case LIQUID_IIRDES_ELLIP:
139 printf("elliptic filter design:\n");
140 nza = 2*L;
141 Gp = powf(10.0f, -Ap / 20.0f);
142 Gs = powf(10.0f, -As / 20.0f);
143 printf(" Gp = %12.8f\n", Gp);
144 printf(" Gs = %12.8f\n", Gs);
145
146 // epsilon values
147 ep = sqrtf(1.0f/(Gp*Gp) - 1.0f);
148 es = sqrtf(1.0f/(Gs*Gs) - 1.0f);
149
150 ellip_azpkf(order,ep,es,za,pa,&ka);
151 break;
152 case LIQUID_IIRDES_BESSEL:
153 printf("Bessel filter design:\n");
154 bessel_azpkf(order,za,pa,&ka);
155 nza = 0;
156 break;
157 default:
158 fprintf(stderr,"error: %s: unknown filter type\n", argv[0]);
159 exit(1);
160 }
161
162 // transform zeros, poles, gain
163 for (i=0; i<npa; i++) {
164 pa[i] *= wc;
165 ka *= wc;
166 }
167 for (i=0; i<nza; i++) {
168 za[i] *= wc;
169 ka /= wc;
170 }
171
172 // print zeros, poles, gain to screen
173 for (i=0; i<nza; i++)
174 printf("z(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(za[i]), cimagf(za[i]));
175 printf("\n");
176 for (i=0; i<npa; i++)
177 printf("p(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(pa[i]), cimagf(pa[i]));
178 printf("\n");
179 printf("ka = %12.4e + j*%12.4e;\n", crealf(ka), cimagf(ka));
180
181 // open output file
182 FILE*fid = fopen(OUTPUT_FILENAME,"w");
183 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
184 fprintf(fid,"clear all;\n");
185 fprintf(fid,"close all;\n");
186 fprintf(fid,"\n");
187 fprintf(fid,"order=%u;\n", order);
188 fprintf(fid,"wc = %12.8f;\n", wc);
189 fprintf(fid,"npa = %u;\n", npa);
190 fprintf(fid,"nza = %u;\n", nza);
191 fprintf(fid,"pa = zeros(1,npa);\n");
192 fprintf(fid,"za = zeros(1,nza);\n");
193 for (i=0; i<npa; i++)
194 fprintf(fid,"pa(%3u) = %16.8e + j*%16.8e;\n", i+1, crealf(pa[i]), cimagf(pa[i]));
195 for (i=0; i<nza; i++)
196 fprintf(fid,"za(%3u) = %16.8e + j*%16.8e;\n", i+1, crealf(za[i]), cimagf(za[i]));
197
198 fprintf(fid,"k = %16.8e;\n", crealf(ka));
199
200 fprintf(fid,"b = 1;\n");
201 fprintf(fid,"for i=1:nza,\n");
202 fprintf(fid," b = conv(b,[1 za(i)]);\n");
203 fprintf(fid,"end;\n");
204
205 fprintf(fid,"a = 1;\n");
206 fprintf(fid,"for i=1:npa,\n");
207 fprintf(fid," a = conv(a,[1 pa(i)]);\n");
208 fprintf(fid,"end;\n");
209
210 fprintf(fid,"b = real(b)*k;\n");
211 fprintf(fid,"a = real(a);\n");
212
213 fprintf(fid,"w = 10.^[-2:0.002:2]*wc;\n");
214 fprintf(fid,"s = j*w;\n");
215 fprintf(fid,"h = polyval(b,s) ./ polyval(a,s);\n");
216 fprintf(fid,"H = 20*log10(abs(h));\n");
217 fprintf(fid,"figure;\n");
218 fprintf(fid,"t = [0:0.02:1]*2*pi;\n");
219 fprintf(fid,"plot(wc*cos(t),wc*sin(t),'-','Color',[1 1 1]*0.7,...\n");
220 fprintf(fid," real(pa),imag(pa),'x');\n");
221 fprintf(fid,"if nza > 0,\n");
222 fprintf(fid," hold on; plot(real(za),imag(za),'ob'); hold off;\n");
223 fprintf(fid," legend('(\\omega_c)','poles','zeros',0);\n");
224 fprintf(fid,"else,\n");
225 fprintf(fid," legend('(\\omega_c)','poles',0);\n");
226 fprintf(fid,"end;\n");
227 fprintf(fid," axis([-1 1 -1 1]*1.2*max([wc abs(pa) abs(za)]));\n");
228 //fprintf(fid," real(za),imag(za),'x');\n");
229 fprintf(fid,"axis square;\n");
230 fprintf(fid,"grid on;\n");
231 fprintf(fid,"xlabel('real');\n");
232 fprintf(fid,"ylabel('imag');\n");
233 fprintf(fid,"\n");
234
235 // plot group delay
236 fprintf(fid,"figure;\n");
237 fprintf(fid,"gd = gradient(unwrap(arg(h)))./gradient(w);\n");
238 fprintf(fid,"semilogx(w,gd);\n");
239 fprintf(fid,"xlabel('Angular frequency, \\omega [rad/s]');\n");
240 fprintf(fid,"ylabel('group delay [s]');\n");
241 fprintf(fid,"gd_min = min(gd); if gd_min < 0, gd_min=0; end;\n");
242 fprintf(fid,"axis([wc/100 wc*100 gd_min 1.1*max(gd)]);\n");
243 fprintf(fid,"grid on;\n");
244
245 // plot magnitude response
246 fprintf(fid,"figure;\n");
247 fprintf(fid,"subplot(2,1,1);\n");
248 fprintf(fid," semilogx(w,H,'-'); grid on;\n");
249 fprintf(fid," axis([wc/100 wc*100 -5 1]);\n");
250 fprintf(fid," xlabel('Angular frequency, \\omega [rad/s]');\n");
251 fprintf(fid," ylabel('PSD [dB]');\n");
252 fprintf(fid,"subplot(2,1,2);\n");
253 fprintf(fid," semilogx(w,H,'-'); grid on;\n");
254 fprintf(fid," axis([wc/100 wc*100 -100 10]);\n");
255 fprintf(fid," xlabel('Angular frequency, \\omega [rad/s]');\n");
256 fprintf(fid," ylabel('PSD [dB]');\n");
257
258 fclose(fid);
259 printf("results written to %s.\n", OUTPUT_FILENAME);
260
261 printf("done.\n");
262 return 0;
263 }
264
265