1 //
2 // firpfbch_crcf_example.c
3 //
4 // Finite impulse response (FIR) polyphase filter bank (PFB)
5 // channelizer example. This example demonstrates the functionality
6 // of the polyphase filter bank channelizer and how its output
7 // is mathematically equivalent to a series of parallel down-
8 // converters (mixers/decimators). Both the synthesis and analysis
9 // filter banks are presented.
10 //
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <math.h>
15 #include <assert.h>
16
17 #include "liquid.h"
18
19 #define OUTPUT_FILENAME "firpfbch_crcf_example.m"
20
main()21 int main() {
22 // options
23 unsigned int num_channels=4; // number of channels
24 unsigned int p=3; // filter length (symbols)
25 unsigned int num_symbols=6; // number of symbols
26
27 // derived values
28 unsigned int num_samples = num_channels * num_symbols;
29
30 unsigned int i;
31 unsigned int j;
32
33 // generate synthesis filter
34 // NOTE : these coefficients can be random; the purpose of this
35 // exercise is to demonstrate mathematical equivalence
36 unsigned int h_len = p*num_channels;
37 float h[h_len];
38 for (i=0; i<h_len; i++) h[i] = randnf();
39
40 // generate analysis filter
41 unsigned int g_len = p*num_channels;
42 float g[g_len];
43 for (i=0; i<g_len; i++)
44 g[i] = h[g_len-i-1];
45
46 // create synthesis/analysis filter objects
47 firfilt_crcf fs = firfilt_crcf_create(h, h_len);
48 firfilt_crcf fa = firfilt_crcf_create(g, g_len);
49
50 // create synthesis/analysis filterbank channelizer objects
51 firpfbch_crcf qs = firpfbch_crcf_create(LIQUID_SYNTHESIZER, num_channels, p, h);
52 firpfbch_crcf qa = firpfbch_crcf_create(LIQUID_ANALYZER, num_channels, p, g);
53
54 float complex x[num_samples]; // random input (noise)
55 float complex Y0[num_symbols][num_channels]; // channelized output (filterbank)
56 float complex Y1[num_symbols][num_channels]; // channelized output
57 float complex z0[num_samples]; // time-domain output (filterbank)
58 float complex z1[num_samples]; // time-domain output
59
60 // generate input sequence (complex noise)
61 for (i=0; i<num_samples; i++)
62 x[i] = randnf() * cexpf(_Complex_I*randf()*2*M_PI);
63
64 //
65 // ANALYZERS
66 //
67
68 //
69 // run analysis filter bank
70 //
71 for (i=0; i<num_symbols; i++)
72 firpfbch_crcf_analyzer_execute(qa, &x[i*num_channels], &Y0[i][0]);
73
74
75 //
76 // run traditional down-converter (inefficient)
77 //
78 float dphi; // carrier frequency
79 unsigned int n=0;
80 for (i=0; i<num_channels; i++) {
81
82 // reset filter
83 firfilt_crcf_reset(fa);
84
85 // set center frequency
86 dphi = 2.0f * M_PI * (float)i / (float)num_channels;
87
88 // reset symbol counter
89 n=0;
90
91 for (j=0; j<num_samples; j++) {
92 // push down-converted sample into filter
93 firfilt_crcf_push(fa, x[j]*cexpf(-_Complex_I*j*dphi));
94
95 // compute output at the appropriate sample time
96 assert(n<num_symbols);
97 if ( ((j+1)%num_channels)==0 ) {
98 firfilt_crcf_execute(fa, &Y1[n][i]);
99 n++;
100 }
101 }
102 assert(n==num_symbols);
103
104 }
105
106
107 //
108 // SYNTHESIZERS
109 //
110
111 //
112 // run synthesis filter bank
113 //
114 for (i=0; i<num_symbols; i++)
115 firpfbch_crcf_synthesizer_execute(qs, &Y0[i][0], &z0[i*num_channels]);
116
117 //
118 // run traditional up-converter (inefficient)
119 //
120
121 // clear output array
122 for (i=0; i<num_samples; i++)
123 z1[i] = 0.0f;
124
125 float complex y_hat;
126 for (i=0; i<num_channels; i++) {
127 // reset filter
128 firfilt_crcf_reset(fs);
129
130 // set center frequency
131 dphi = 2.0f * M_PI * (float)i / (float)num_channels;
132
133 // reset input symbol counter
134 n=0;
135
136 for (j=0; j<num_samples; j++) {
137
138 // interpolate sequence
139 if ( (j%num_channels)==0 ) {
140 assert(n<num_symbols);
141 firfilt_crcf_push(fs, Y1[n][i]);
142 n++;
143 } else {
144 firfilt_crcf_push(fs, 0);
145 }
146 firfilt_crcf_execute(fs, &y_hat);
147
148 // accumulate up-converted sample
149 z1[j] += y_hat * cexpf(_Complex_I*j*dphi);
150 }
151 assert(n==num_symbols);
152 }
153
154 // destroy objects
155 firfilt_crcf_destroy(fs);
156 firfilt_crcf_destroy(fa);
157 firpfbch_crcf_destroy(qs);
158 firpfbch_crcf_destroy(qa);
159
160
161 //
162 // RESULTS
163 //
164
165 //
166 // analyzers
167 //
168
169 // print filterbank channelizer
170 printf("\n");
171 printf("filterbank channelizer:\n");
172 for (i=0; i<num_symbols; i++) {
173 printf("%3u: ", i);
174 for (j=0; j<num_channels; j++) {
175 printf(" %8.5f+j%8.5f, ", crealf(Y0[i][j]), cimagf(Y0[i][j]));
176 }
177 printf("\n");
178 }
179
180 // print traditional channelizer
181 printf("\n");
182 printf("traditional channelizer:\n");
183 for (i=0; i<num_symbols; i++) {
184 printf("%3u: ", i);
185 for (j=0; j<num_channels; j++) {
186 printf(" %8.5f+j%8.5f, ", crealf(Y1[i][j]), cimagf(Y1[i][j]));
187 }
188 printf("\n");
189 }
190
191
192 float mse_analyzer[num_channels];
193 float complex d;
194 for (i=0; i<num_channels; i++) {
195 mse_analyzer[i] = 0.0f;
196 for (j=0; j<num_symbols; j++) {
197 d = Y0[j][i] - Y1[j][i];
198 mse_analyzer[i] += crealf(d*conjf(d));
199 }
200
201 mse_analyzer[i] /= num_symbols;
202 }
203 printf("\n");
204 printf("rmse: ");
205 for (i=0; i<num_channels; i++)
206 printf("%12.4e ", sqrt(mse_analyzer[i]));
207 printf("\n");
208
209
210 //
211 // synthesizers
212 //
213
214 printf("\n");
215 printf("output: filterbank: traditional:\n");
216 for (i=0; i<num_samples; i++) {
217 printf("%3u: %10.5f+%10.5fj %10.5f+%10.5fj\n",
218 i,
219 crealf(z0[i]), cimagf(z0[i]),
220 crealf(z1[i]), cimagf(z1[i]));
221 }
222
223 float mse_synthesizer = 0.0f;
224 for (i=0; i<num_samples; i++) {
225 d = z0[i] - z1[i];
226 mse_synthesizer += crealf(d*conjf(d));
227 }
228 mse_synthesizer /= num_samples;
229 printf("\n");
230 printf("rmse: %12.4e\n", sqrtf(mse_synthesizer));
231
232 //
233 // EXPORT DATA TO FILE
234 //
235
236 FILE * fid = fopen(OUTPUT_FILENAME,"w");
237 fprintf(fid,"%% %s: auto-generated file\n\n", OUTPUT_FILENAME);
238 fprintf(fid,"clear all;\n");
239 fprintf(fid,"close all;\n");
240 fprintf(fid,"num_channels=%u;\n", num_channels);
241 fprintf(fid,"num_symbols=%u;\n", num_symbols);
242 fprintf(fid,"num_samples = num_channels*num_symbols;\n");
243
244 fprintf(fid,"x = zeros(1,num_samples);\n");
245 fprintf(fid,"y0 = zeros(num_symbols,num_channels);\n");
246 fprintf(fid,"y1 = zeros(num_symbols,num_channels);\n");
247 fprintf(fid,"z0 = zeros(1,num_samples);\n");
248 fprintf(fid,"z1 = zeros(1,num_samples);\n");
249
250 // input
251 for (i=0; i<num_samples; i++)
252 fprintf(fid,"x(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(x[i]), cimag(x[i]));
253
254 // analysis
255 for (i=0; i<num_symbols; i++) {
256 for (j=0; j<num_channels; j++) {
257 fprintf(fid,"y0(%4u,%4u) = %12.4e + j*%12.4e;\n", i+1, j+1, crealf(Y0[i][j]), cimag(Y0[i][j]));
258 fprintf(fid,"y1(%4u,%4u) = %12.4e + j*%12.4e;\n", i+1, j+1, crealf(Y1[i][j]), cimag(Y1[i][j]));
259 }
260 }
261
262 // synthesis
263 for (i=0; i<num_samples; i++) {
264 fprintf(fid,"z0(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(z0[i]), cimag(z0[i]));
265 fprintf(fid,"z1(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(z1[i]), cimag(z1[i]));
266 }
267 fprintf(fid,"z0 = z0 / num_channels;\n");
268 fprintf(fid,"z1 = z1 / num_channels;\n");
269
270 // plot results
271 fprintf(fid,"\n\n");
272 fprintf(fid,"ts = 0:(num_symbols-1);\n");
273 fprintf(fid,"for i=1:num_channels,\n");
274 fprintf(fid,"figure;\n");
275 fprintf(fid,"title(['channel ' num2str(i)]);\n");
276 fprintf(fid,"subplot(2,1,1);\n");
277 fprintf(fid," plot(ts,real(y0(:,i)),'-x', ts,real(y1(:,i)),'s');\n");
278 //fprintf(fid," axis([0 (num_symbols-1) -2 2]);\n");
279 fprintf(fid,"subplot(2,1,2);\n");
280 fprintf(fid," plot(ts,imag(y0(:,i)),'-x', ts,imag(y1(:,i)),'s');\n");
281 //fprintf(fid," axis([0 (num_symbols-1) -2 2]);\n");
282 fprintf(fid,"end;\n");
283
284 fprintf(fid,"t = 0:(num_samples-1);\n");
285 fprintf(fid,"figure;\n");
286 fprintf(fid,"title('composite');\n");
287 fprintf(fid,"subplot(2,1,1);\n");
288 fprintf(fid," plot(t,real(z0),'-x', t,real(z1),'s');\n");
289 fprintf(fid,"subplot(2,1,2);\n");
290 fprintf(fid," plot(t,imag(z0),'-x', t,imag(z1),'s');\n");
291
292
293 fclose(fid);
294 printf("results written to '%s'\n", OUTPUT_FILENAME);
295
296 printf("done.\n");
297 return 0;
298 }
299