1 //
2 // iirfilt_cccf_example.c
3 //
4 // Complex infinite impulse response filter example. Demonstrates the
5 // functionality of iirfilt with complex coefficients by designing a
6 // filter with specified parameters and then filters noise.
7 //
8
9 #include <stdlib.h>
10 #include <stdio.h>
11 #include <string.h>
12 #include <math.h>
13 #include <complex.h>
14 #include <getopt.h>
15
16 #include "liquid.h"
17
18 #define OUTPUT_FILENAME "iirfilt_cccf_example.m"
19
20 // print usage/help message
usage()21 void usage()
22 {
23 printf("iirfilt_cccf_example -- infinite impulse response filter example\n");
24 printf("options (default values in []):\n");
25 printf(" h : print help\n");
26 printf(" t : filter type: [butter], cheby1, cheby2, ellip, bessel\n");
27 printf(" b : filter transformation: [LP], HP, BP, BS\n");
28 printf(" n : filter order, n > 0 [5]\n");
29 printf(" r : passband ripple in dB (cheby1, ellip), r > 0 [1.0]\n");
30 printf(" s : stopband attenuation in dB (cheby2, ellip), s > 0 [40.0]\n");
31 printf(" f : passband cut-off, 0 < f < 0.5 [0.2]\n");
32 printf(" c : center frequency (BP, BS cases), 0 < c < 0.5 [0.25]\n");
33 printf(" o : format [sos], tf\n");
34 printf(" sos : second-order sections form\n");
35 printf(" tf : regular transfer function form (potentially\n");
36 printf(" unstable for large orders\n");
37 }
38
main(int argc,char * argv[])39 int main(int argc, char*argv[])
40 {
41 // options
42 unsigned int order=4; // filter order
43 float fc=0.1f; // cutoff frequency
44 float f0=0.0f; // center frequency
45 float Ap=1.0f; // pass-band ripple
46 float As=40.0f; // stop-band attenuation
47 unsigned int n=128; // number of samples
48 liquid_iirdes_filtertype ftype = LIQUID_IIRDES_ELLIP;
49 liquid_iirdes_bandtype btype = LIQUID_IIRDES_LOWPASS;
50 liquid_iirdes_format format = LIQUID_IIRDES_SOS;
51
52 int dopt;
53 while ((dopt = getopt(argc,argv,"ht:b:n:r:s:f:c:o:")) != EOF) {
54 switch (dopt) {
55 case 'h': usage(); return 0;
56 case 't':
57 if (strcmp(optarg,"butter")==0) ftype = LIQUID_IIRDES_BUTTER;
58 else if (strcmp(optarg,"cheby1")==0) ftype = LIQUID_IIRDES_CHEBY1;
59 else if (strcmp(optarg,"cheby2")==0) ftype = LIQUID_IIRDES_CHEBY2;
60 else if (strcmp(optarg,"ellip") ==0) ftype = LIQUID_IIRDES_ELLIP;
61 else if (strcmp(optarg,"bessel")==0) ftype = LIQUID_IIRDES_BESSEL;
62 else {
63 fprintf(stderr,"error: iirdes_example, unknown filter type '%s'\n", optarg);
64 exit(1);
65 }
66 break;
67 case 'b':
68 if (strcmp(optarg,"LP")==0) btype = LIQUID_IIRDES_LOWPASS;
69 else if (strcmp(optarg,"HP")==0) btype = LIQUID_IIRDES_HIGHPASS;
70 else if (strcmp(optarg,"BP")==0) btype = LIQUID_IIRDES_BANDPASS;
71 else if (strcmp(optarg,"BS")==0) btype = LIQUID_IIRDES_BANDSTOP;
72 else {
73 fprintf(stderr,"error: iirdes_example, unknown band type '%s'\n", optarg);
74 exit(1);
75 }
76 break;
77 case 'n': order = atoi(optarg); break;
78 case 'r': Ap = atof(optarg); break;
79 case 's': As = atof(optarg); break;
80 case 'f': fc = atof(optarg); break;
81 case 'c': f0 = atof(optarg); break;
82 case 'o':
83 if (strcmp(optarg,"sos")==0) format = LIQUID_IIRDES_SOS;
84 else if (strcmp(optarg,"tf") ==0) format = LIQUID_IIRDES_TF;
85 else {
86 fprintf(stderr,"error: iirdes_example, unknown output format '%s'\n", optarg);
87 exit(1);
88 }
89 break;
90 default:
91 exit(1);
92 }
93 }
94
95 // design filter from prototype
96 iirfilt_cccf q = iirfilt_cccf_create_prototype(ftype, btype, format, order, fc, f0, Ap, As);
97 iirfilt_cccf_print(q);
98
99 unsigned int i;
100
101 // allocate memory for data arrays
102 float complex x[n];
103 float complex y[n];
104
105 // generate input signal (noisy sine wave with decaying amplitude)
106 unsigned int wlen = (3*n)/4;
107 for (i=0; i<n; i++) {
108 // input signal (windowed noise)
109 x[i] = randnf() + _Complex_I*randnf();
110 x[i] *= i < wlen ? hamming(i,wlen) : 0.0f;
111
112 // run filter
113 iirfilt_cccf_execute(q, x[i], &y[i]);
114 }
115
116 // compute two-sided frequency response
117 unsigned int nfft=512;
118 float complex H[nfft];
119 for (i=0; i<nfft; i++) {
120 float freq = (float)i / (float)nfft - 0.5f;
121 iirfilt_cccf_freqresponse(q, freq, &H[i]);
122 }
123
124 // destroy filter object
125 iirfilt_cccf_destroy(q);
126
127 //
128 // plot results to output file
129 //
130 FILE * fid = fopen(OUTPUT_FILENAME,"w");
131 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
132 fprintf(fid,"clear all;\n");
133 fprintf(fid,"close all;\n");
134 fprintf(fid,"\n");
135 fprintf(fid,"order=%u;\n", order);
136
137 // save input, output arrays
138 fprintf(fid,"n=%u;\n",n);
139 fprintf(fid,"x=zeros(1,n);\n");
140 fprintf(fid,"y=zeros(1,n);\n");
141 for (i=0; i<n; i++) {
142 fprintf(fid,"x(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(x[i]), cimagf(x[i]));
143 fprintf(fid,"y(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
144 }
145
146 // save frequency response
147 fprintf(fid,"nfft=%u;\n",nfft);
148 fprintf(fid,"H=zeros(1,nfft);\n");
149 for (i=0; i<nfft; i++)
150 fprintf(fid,"H(%4u) = %12.8f + j*%12.8f;\n", i+1, crealf(H[i]), cimagf(H[i]));
151
152 // plot time-domain output
153 fprintf(fid,"t=0:(n-1);\n");
154 fprintf(fid,"figure;\n");
155 fprintf(fid,"subplot(2,1,1);\n");
156 fprintf(fid," plot(t,real(x),'-','Color',[1 1 1]*0.5,'LineWidth',1,...\n");
157 fprintf(fid," t,real(y),'-','Color',[0 0.5 0.25],'LineWidth',2);\n");
158 fprintf(fid," xlabel('time');\n");
159 fprintf(fid," ylabel('real');\n");
160 fprintf(fid," legend('input','filtered output',1);\n");
161 fprintf(fid," grid on;\n");
162 fprintf(fid,"subplot(2,1,2);\n");
163 fprintf(fid," plot(t,imag(x),'-','Color',[1 1 1]*0.5,'LineWidth',1,...\n");
164 fprintf(fid," t,imag(y),'-','Color',[0 0.25 0.5],'LineWidth',2);\n");
165 fprintf(fid," xlabel('time');\n");
166 fprintf(fid," ylabel('imag');\n");
167 fprintf(fid," legend('input','filtered output',1);\n");
168 fprintf(fid," grid on;\n");
169
170 // plot spectral output
171 fprintf(fid,"X = 20*log10(abs(fftshift(fft(x))));\n");
172 fprintf(fid,"Y = 20*log10(abs(fftshift(fft(y))));\n");
173 fprintf(fid,"figure;\n");
174 fprintf(fid,"plot([0:(n-1)]/n-0.5, X, 'Color', [1 1 1]*0.5,\n");
175 fprintf(fid," [0:(n-1)]/n-0.5, Y, 'Color', [0 0.25 0.50]);\n");
176 fprintf(fid,"xlabel('Normalized Frequency');\n");
177 fprintf(fid,"ylabel('PSD [dB]');\n");
178 fprintf(fid,"legend('input','filtered output',1);\n");
179 fprintf(fid,"grid on;\n");
180
181 // plot ideal frequency response
182 fprintf(fid,"f=[0:(nfft-1)]/nfft - 0.5;\n");
183 fprintf(fid,"figure;\n");
184 fprintf(fid,"subplot(3,1,1);\n");
185 fprintf(fid," plot(f,20*log10(abs(H)));\n");
186 fprintf(fid," axis([-0.5 0.5 -3 0]);\n");
187 fprintf(fid," grid on;\n");
188 fprintf(fid," legend('Pass band (dB)',0);\n");
189 fprintf(fid,"subplot(3,1,2);\n");
190 fprintf(fid," plot(f,20*log10(abs(H)));\n");
191 fprintf(fid," axis([-0.5 0.5 -100 0]);\n");
192 fprintf(fid," grid on;\n");
193 fprintf(fid," legend('Stop band (dB)',0);\n");
194 fprintf(fid,"subplot(3,1,3);\n");
195 fprintf(fid," plot(f,180/pi*arg(H));\n");
196 fprintf(fid," axis([-0.5 0.5 -180 180]);\n");
197 fprintf(fid," grid on;\n");
198 fprintf(fid," legend('Phase (degrees)',0);\n");
199 fclose(fid);
200 printf("results written to %s.\n", OUTPUT_FILENAME);
201
202 printf("done.\n");
203 return 0;
204 }
205
206