1 //
2 // matched_filter_example.c
3 //
4
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <string.h>
8 #include <getopt.h>
9 #include <math.h>
10
11 #include "liquid.h"
12
13 #define OUTPUT_FILENAME "matched_filter_example.m"
14
15 // print usage/help message
usage()16 void usage()
17 {
18 printf("matched_filter_example options:\n");
19 printf(" u/h : print usage/help\n");
20 printf(" t : filter type: [rrcos], rkaiser, arkaiser, hM3, gmsk, fexp, fsech, farcsech\n");
21 printf(" k : filter samples/symbol, k >= 2, default: 2\n");
22 printf(" m : filter delay (symbols), m >= 1, default: 3\n");
23 printf(" b : filter excess bandwidth factor, 0 < b < 1, default: 0.5\n");
24 printf(" n : number of symbols, default: 16\n");
25 }
26
27
main(int argc,char * argv[])28 int main(int argc, char*argv[]) {
29 // options
30 unsigned int k=2; // samples/symbol
31 unsigned int m=3; // symbol delay
32 float beta=0.7f; // excess bandwidth factor
33 unsigned int num_symbols=16;
34 int ftype_tx = LIQUID_FIRFILT_RRC;
35 int ftype_rx = LIQUID_FIRFILT_RRC;
36
37 int dopt;
38 while ((dopt = getopt(argc,argv,"uht:k:m:b:n:")) != EOF) {
39 switch (dopt) {
40 case 'u':
41 case 'h': usage(); return 0;
42 case 't':
43 if (strcmp(optarg,"gmsk")==0) {
44 ftype_tx = LIQUID_FIRFILT_GMSKTX;
45 ftype_rx = LIQUID_FIRFILT_GMSKRX;
46 } else {
47 ftype_tx = liquid_getopt_str2firfilt(optarg);
48 ftype_rx = liquid_getopt_str2firfilt(optarg);
49 }
50 if (ftype_tx == LIQUID_FIRFILT_UNKNOWN) {
51 fprintf(stderr,"error: %s, unknown filter type '%s'\n", argv[0], optarg);
52 exit(1);
53 }
54 break;
55 case 'k': k = atoi(optarg); break;
56 case 'm': m = atoi(optarg); break;
57 case 'b': beta = atof(optarg); break;
58 case 'n': num_symbols = atoi(optarg); break;
59 default:
60 exit(1);
61 }
62 }
63
64 if (k < 2) {
65 fprintf(stderr,"error: %s, k must be at least 2\n", argv[0]);
66 exit(1);
67 } else if (m < 1) {
68 fprintf(stderr,"error: %s, m must be at least 1\n", argv[0]);
69 exit(1);
70 } else if (beta <= 0.0f || beta >= 1.0f) {
71 fprintf(stderr,"error: %s, beta must be in (0,1)\n", argv[0]);
72 exit(1);
73 }
74
75 unsigned int i;
76
77 // derived values
78 unsigned int num_samples = num_symbols*k;
79 unsigned int h_len = 2*k*m+1; // transmit/receive filter length
80 unsigned int hc_len = 4*k*m+1; // composite filter length
81
82 // arrays
83 float ht[h_len]; // transmit filter
84 float hr[h_len]; // receive filter
85 float hc[hc_len]; // composite filter
86
87 // design the filter(s)
88 liquid_firdes_prototype(ftype_tx, k, m, beta, 0, ht);
89 liquid_firdes_prototype(ftype_rx, k, m, beta, 0, hr);
90
91 for (i=0; i<h_len; i++) printf("ht(%3u) = %12.8f;\n", i+1, ht[i]);
92 for (i=0; i<h_len; i++) printf("hr(%3u) = %12.8f;\n", i+1, hr[i]);
93
94 #if 0
95 // generate receive filter coefficients (reverse of transmit)
96 float hr[h_len];
97 for (i=0; i<h_len; i++)
98 hr[i] = ht[h_len-i-1];
99 #endif
100
101 // compute composite filter response
102 for (i=0; i<4*k*m+1; i++) {
103 int lag = (int)i - (int)(2*k*m);
104 hc[i] = liquid_filter_crosscorr(ht,h_len, hr,h_len, lag);
105 }
106
107 // compute filter inter-symbol interference
108 float rxy0 = liquid_filter_crosscorr(ht,h_len, hr,h_len, 0);
109 float isi_rms=0;
110 for (i=1; i<2*m; i++) {
111 float e = liquid_filter_crosscorr(ht,h_len, hr,h_len, i*k) / rxy0;
112 isi_rms += e*e;
113 }
114 isi_rms = sqrtf( isi_rms / (float)(2*m-1) );
115 printf(" isi (rms) : %12.8f dB\n", 20*log10f(isi_rms));
116
117 // compute relative stop-band energy
118 unsigned int nfft = 2048;
119 float As = liquid_filter_energy(ht, h_len, 0.5f*(1.0f + beta)/(float)k, nfft);
120 printf(" As : %12.8f dB\n", 20*log10f(As));
121
122 // generate signal
123 float sym_in[num_symbols];
124 float y[num_samples];
125 float sym_out[num_symbols];
126
127 // create interpolator and decimator
128 firinterp_rrrf interp = firinterp_rrrf_create(k, ht, h_len);
129 firdecim_rrrf decim = firdecim_rrrf_create( k, hr, h_len);
130
131 for (i=0; i<num_symbols; i++) {
132 // generate random symbol
133 sym_in[i] = (rand() % 2) ? 1.0f : -1.0f;
134
135 // interpolate
136 firinterp_rrrf_execute(interp, sym_in[i], &y[i*k]);
137
138 // decimate
139 firdecim_rrrf_execute(decim, &y[i*k], &sym_out[i]);
140
141 // normalize output
142 sym_out[i] /= k;
143
144 printf(" %3u : %8.5f", i, sym_out[i]);
145 if (i>=2*m) printf(" *\n");
146 else printf("\n");
147 }
148
149 // clean up objects
150 firinterp_rrrf_destroy(interp);
151 firdecim_rrrf_destroy(decim);
152
153 //
154 // export results
155 //
156 FILE * fid = fopen(OUTPUT_FILENAME,"w");
157 fprintf(fid,"%% %s : auto-generated file\n\n", OUTPUT_FILENAME);
158 fprintf(fid,"clear all;\n");
159 fprintf(fid,"close all;\n");
160 fprintf(fid,"k = %u;\n", k);
161 fprintf(fid,"m = %u;\n", m);
162 fprintf(fid,"beta = %12.8f;\n", beta);
163 fprintf(fid,"num_symbols = %u;\n", num_symbols);
164 fprintf(fid,"num_samples = k*num_symbols;\n");
165
166 fprintf(fid,"y = zeros(1,num_samples);\n");
167 for (i=0; i<num_samples; i++)
168 fprintf(fid," y(%3u) = %12.8f;\n", i+1, y[i]);
169
170 for (i=0; i<h_len; i++) fprintf(fid,"ht(%3u) = %20.8e;\n", i+1, ht[i]);
171 for (i=0; i<h_len; i++) fprintf(fid,"hr(%3u) = %20.8e;\n", i+1, hr[i]);
172 for (i=0; i<hc_len; i++) fprintf(fid,"hc(%3u) = %20.8e;\n", i+1, hc[i]);
173
174 fprintf(fid,"nfft=1024;\n");
175 fprintf(fid,"f = [0:(nfft-1)]/nfft - 0.5;\n");
176 fprintf(fid,"Ht = 20*log10(abs(fftshift(fft(ht/k, nfft))));\n");
177 fprintf(fid,"Hr = 20*log10(abs(fftshift(fft(hr/k, nfft))));\n");
178 fprintf(fid,"Hc = 20*log10(abs(fftshift(fft(hc/k^2, nfft))));\n");
179 fprintf(fid,"figure;\n");
180 fprintf(fid,"plot(f,Hc,'-','LineWidth',2,...\n");
181 fprintf(fid," [0.5/k],[-6],'or',...\n");
182 fprintf(fid," [0.5/k*(1-beta) 0.5/k*(1-beta)],[-100 10],'-r',...\n");
183 fprintf(fid," [0.5/k*(1+beta) 0.5/k*(1+beta)],[-100 10],'-r');\n");
184 fprintf(fid,"xlabel('normalized frequency');\n");
185 fprintf(fid,"ylabel('PSD');\n");
186 fprintf(fid,"axis([-0.5 0.5 -100 10]);\n");
187 fprintf(fid,"grid on;\n");
188
189 // compute composite filter
190 fprintf(fid,"figure;\n");
191 fprintf(fid,"hc = conv(ht,hr)/k;\n");
192 fprintf(fid,"t = [(-2*k*m):(2*k*m)]/k;\n");
193 fprintf(fid,"i0 = [0:k:4*k*m]+1;\n");
194 fprintf(fid,"plot(t, hc, '-s',...\n");
195 fprintf(fid," t(i0),hc(i0),'or');\n");
196 fprintf(fid,"xlabel('symbol index');\n");
197 fprintf(fid,"ylabel('matched filter response');\n");
198 fprintf(fid,"grid on;\n");
199
200 fclose(fid);
201 printf("results written to %s.\n", OUTPUT_FILENAME);
202
203 printf("done.\n");
204 return 0;
205 }
206
207