1 //
2 // sandbox/ofdmoqam_firpfbch_cfo_test.c
3 //
4 // Tests the validity of OFDM/OQAM carrier frequency offset
5 // recovery using firpfbch_cccf channelizer objects.
6 //
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <string.h>
12 #include <complex.h>
13
14 #include "liquid.h"
15
16 #define OUTPUT_FILENAME "ofdmoqam_firpfbch_cfo_test.m"
17
18 #define DEBUG 1
19
20
main()21 int main() {
22 // options
23 unsigned int num_channels=64; // must be even number
24 unsigned int num_symbols=16; // number of symbols
25 unsigned int m=3; // filter delay (symbols)
26 float beta = 0.9f; // filter excess bandwidth factor
27 float phi = 0.0f; // carrier phase offset;
28 float dphi = 0.04f; // carrier frequency offset
29
30 // number of frames (compensate for filter delay)
31 unsigned int num_frames = num_symbols + 2*m;
32 unsigned int num_samples = num_channels * num_frames;
33 unsigned int i;
34 unsigned int j;
35
36 // create filter prototype
37 unsigned int h_len = 2*num_channels*m + 1;
38 float h[h_len];
39 float complex hc[h_len];
40 float complex gc[h_len];
41 liquid_firdes_rkaiser(num_channels, m, beta, 0.0f, h);
42 unsigned int g_len = 2*num_channels*m;
43 for (i=0; i<g_len; i++) {
44 hc[i] = h[i];
45 gc[i] = h[g_len-i-1] * cexpf(_Complex_I*dphi*i);
46 }
47
48 // data arrays
49 float complex s[num_channels]; // input symbols
50 float complex y[num_samples]; // time-domain samples
51 float complex Y0[num_frames][num_channels]; // channelized output
52 float complex Y1[num_frames][num_channels]; // channelized output
53
54 // create ofdm/oqam generator object and generate data
55 ofdmoqam qs = ofdmoqam_create(num_channels, m, beta, 0.0f, LIQUID_SYNTHESIZER, 0);
56 for (i=0; i<num_frames; i++) {
57 for (j=0; j<num_channels; j++) {
58 if (i<num_symbols) {
59 #if 0
60 // QPSK on all subcarriers
61 s[j] = (rand() % 2 ? 1.0f : -1.0f) +
62 (rand() % 2 ? 1.0f : -1.0f) * _Complex_I;
63 s[j] *= 1.0f / sqrtf(2.0f);
64 #else
65 // BPSK on even subcarriers
66 s[j] = rand() % 2 ? 1.0f : -1.0f;
67 s[j] *= (j%2)==0 ? 1.0f : 0.0f;
68 #endif
69 } else {
70 s[j] = 0.0f;
71 }
72 }
73
74 // run synthesizer
75 ofdmoqam_execute(qs, s, &y[i*num_channels]);
76 }
77 ofdmoqam_destroy(qs);
78
79 // channel
80 for (i=0; i<num_samples; i++)
81 y[i] *= cexpf(_Complex_I*(phi + dphi*i));
82
83
84 //
85 // analysis filterbank (receiver)
86 //
87
88 // create filterbank manually
89 dotprod_cccf dp[num_channels]; // vector dot products
90 windowcf w[num_channels]; // window buffers
91
92 #if DEBUG
93 // print coefficients
94 printf("h_prototype:\n");
95 for (i=0; i<h_len; i++)
96 printf(" h[%3u] = %12.8f\n", i, h[i]);
97 #endif
98
99 // create objects
100 unsigned int gc_sub_len = 2*m;
101 float complex gc_sub[gc_sub_len];
102 for (i=0; i<num_channels; i++) {
103 // sub-sample prototype filter, loading coefficients in
104 // reverse order
105 #if 0
106 for (j=0; j<gc_sub_len; j++)
107 gc_sub[j] = h[j*num_channels+i];
108 #else
109 for (j=0; j<gc_sub_len; j++)
110 gc_sub[gc_sub_len-j-1] = gc[j*num_channels+i];
111 #endif
112
113 // create window buffer and dotprod objects
114 dp[i] = dotprod_cccf_create(gc_sub, gc_sub_len);
115 w[i] = windowcf_create(gc_sub_len);
116
117 #if DEBUG
118 printf("gc_sub[%u] : \n", i);
119 for (j=0; j<gc_sub_len; j++)
120 printf(" g[%3u] = %12.8f + %12.8f\n", j, crealf(gc_sub[j]), cimagf(gc_sub[j]));
121 #endif
122 }
123
124 // generate DFT object
125 float complex x[num_channels]; // time-domain buffer
126 float complex X[num_channels]; // freq-domain buffer
127 #if 0
128 fftplan fft = fft_create_plan(num_channels, X, x, LIQUID_FFT_BACKWARD, 0);
129 #else
130 fftplan fft = fft_create_plan(num_channels, X, x, LIQUID_FFT_FORWARD, 0);
131 #endif
132
133 //
134 // run analysis filter bank
135 //
136 #if 0
137 unsigned int filter_index = 0;
138 #else
139 unsigned int filter_index = num_channels-1;
140 #endif
141 float complex y_hat; // input sample
142 float complex * r; // read pointer
143 for (i=0; i<num_frames; i++) {
144
145 // load buffers
146 for (j=0; j<num_channels; j++) {
147 // grab sample
148 y_hat = y[i*num_channels + j];
149
150 // push sample into buffer at filter index
151 windowcf_push(w[filter_index], y_hat);
152
153 // decrement filter index
154 filter_index = (filter_index + num_channels - 1) % num_channels;
155 //filter_index = (filter_index + 1) % num_channels;
156 }
157
158 // execute filter outputs, reversing order of output (not
159 // sure why this is necessary)
160 for (j=0; j<num_channels; j++) {
161 windowcf_read(w[j], &r);
162 dotprod_cccf_execute(dp[j], r, &X[num_channels-j-1]);
163 }
164
165 #if 1
166 // compensate for carrier frequency offset (before transform)
167 for (j=0; j<num_channels; j++) {
168 X[j] *= cexpf(-_Complex_I*(dphi*i*num_channels));
169 }
170 #endif
171
172 // execute DFT, store result in buffer 'x'
173 fft_execute(fft);
174
175 #if 0
176 // compensate for carrier frequency offset (after transform)
177 for (j=0; j<num_channels; j++) {
178 x[j] *= cexpf(-_Complex_I*(dphi*i*num_channels));
179 }
180 #endif
181
182 // move to output array
183 for (j=0; j<num_channels; j++)
184 Y0[i][j] = x[j];
185 }
186
187
188 // destroy objects
189 for (i=0; i<num_channels; i++) {
190 dotprod_cccf_destroy(dp[i]);
191 windowcf_destroy(w[i]);
192 }
193 fft_destroy_plan(fft);
194
195 #if 0
196 // print filterbank channelizer
197 printf("\n");
198 printf("filterbank channelizer:\n");
199 for (i=0; i<num_symbols; i++) {
200 printf("%3u: ", i);
201 for (j=0; j<num_channels; j++) {
202 printf(" %8.5f+j%8.5f, ", crealf(Y0[i][j]), cimagf(Y0[i][j]));
203 }
204 printf("\n");
205 }
206 #endif
207
208 //
209 // export data
210 //
211 FILE*fid = fopen(OUTPUT_FILENAME,"w");
212 fprintf(fid,"%% %s: auto-generated file\n\n", OUTPUT_FILENAME);
213 fprintf(fid,"clear all;\nclose all;\n\n");
214 fprintf(fid,"num_channels=%u;\n", num_channels);
215 fprintf(fid,"num_symbols=%u;\n", num_symbols);
216 fprintf(fid,"num_frames = %u;\n", num_frames);
217 fprintf(fid,"num_samples = num_frames*num_channels;\n");
218
219 fprintf(fid,"y = zeros(1,%u);\n", num_samples);
220 fprintf(fid,"Y0 = zeros(%u,%u);\n", num_frames, num_channels);
221 fprintf(fid,"Y1 = zeros(%u,%u);\n", num_frames, num_channels);
222
223 for (i=0; i<num_frames; i++) {
224 for (j=0; j<num_channels; j++) {
225 fprintf(fid,"Y0(%4u,%4u) = %12.4e + j*%12.4e;\n", i+1, j+1, crealf(Y0[i][j]), cimagf(Y0[i][j]));
226 fprintf(fid,"Y1(%4u,%4u) = %12.4e + j*%12.4e;\n", i+1, j+1, crealf(Y1[i][j]), cimagf(Y1[i][j]));
227 }
228 }
229
230 // plot BPSK results
231 fprintf(fid,"figure;\n");
232 fprintf(fid,"plot(Y0(:,1:2:end),'x');\n");
233 fprintf(fid,"axis([-1 1 -1 1]*1.2*sqrt(num_channels));\n");
234 fprintf(fid,"axis square;\n");
235 fprintf(fid,"grid on;\n");
236
237 fclose(fid);
238 printf("results written to '%s'\n", OUTPUT_FILENAME);
239
240 printf("done.\n");
241 return 0;
242 }
243
244