1 //
2 // firdes_gmskrx_test.c : test synthesis of GMSK receive filter
3 //
4
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <string.h>
8 #include <math.h>
9 #include <getopt.h>
10
11 #include "liquid.h"
12
13 #define OUTPUT_FILENAME "firdes_gmskrx_test.m"
14
15 // print usage/help message
usage()16 void usage()
17 {
18 printf("Usage: sandbox/firdes_gmskrx_test [OPTION]\n");
19 printf("Run example GMSK receive filter design\n");
20 printf("\n");
21 printf(" u/h : print usage/help\n");
22 printf(" k : filter samples/symbol, k >= 2, default: 4\n");
23 printf(" m : filter delay (symbols), m >= 1, default: 3\n");
24 printf(" b : filter excess bandwidth factor, 0 < b < 1, default: 0.3\n");
25 }
26
main(int argc,char * argv[])27 int main(int argc, char*argv[]) {
28 // options
29 unsigned int k=4; // samples/symbol
30 unsigned int m=3; // filter delay [symbols]
31 float BT = 0.3f; // bandwidth-time product
32
33 // read properties from command line
34 int dopt;
35 while ((dopt = getopt(argc,argv,"uhk:m:b:")) != EOF) {
36 switch (dopt) {
37 case 'u':
38 case 'h': usage(); return 0;
39 case 'k': k = atoi(optarg); break;
40 case 'm': m = atoi(optarg); break;
41 case 'b': BT = atof(optarg); break;
42 default:
43 exit(1);
44 }
45 }
46
47
48 // validate input
49 if (k < 2) {
50 fprintf(stderr,"error: %s, k must be at least 2\n", argv[0]);
51 exit(1);
52 } else if (m < 1) {
53 fprintf(stderr,"error: %s, m must be at least 1\n", argv[0]);
54 exit(1);
55 } else if (BT <= 0.0f || BT >= 1.0f) {
56 fprintf(stderr,"error: %s, BT must be in (0,1)\n", argv[0]);
57 exit(1);
58 }
59
60
61 unsigned int i;
62
63 // derived values
64 unsigned int h_len = 2*k*m+1; // filter length
65
66 // arrays
67 float ht[h_len]; // transmit filter coefficients
68 float hr[h_len]; // recieve filter coefficients
69
70 // design transmit filter
71 liquid_firdes_gmsktx(k,m,BT,0.0f,ht);
72
73 // print results to screen
74 printf("gmsk transmit filter:\n");
75 for (i=0; i<h_len; i++)
76 printf(" ht(%3u) = %12.8f;\n", i+1, ht[i]);
77
78 //
79 // start of filter design procedure
80 //
81
82 float beta = BT; // prototype filter cut-off
83 float delta = 1e-2f; // filter design correction factor
84
85 // temporary arrays
86 float h_primef[h_len]; // temporary buffer for real coefficients
87 float g_primef[h_len]; //
88
89 float complex h_tx[h_len]; // impulse response of transmit filter
90 float complex h_prime[h_len]; // impulse response of 'prototype' filter
91 float complex g_prime[h_len]; // impulse response of 'gain' filter
92 float complex h_hat[h_len]; // impulse response of receive filter
93
94 float complex H_tx[h_len]; // frequency response of transmit filter
95 float complex H_prime[h_len]; // frequency response of 'prototype' filter
96 float complex G_prime[h_len]; // frequency response of 'gain' filter
97 float complex H_hat[h_len]; // frequency response of receive filter
98
99 // create 'prototype' matched filter
100 // for now use raised-cosine
101 liquid_firdes_prototype(LIQUID_FIRFILT_RCOS,k,m,beta,0.0f,h_primef);
102
103 // create 'gain' filter to improve stop-band rejection
104 float fc = (0.7f + 0.1*beta) / (float)k;
105 float As = 60.0f;
106 liquid_firdes_kaiser(h_len, fc, As, 0.0f, g_primef);
107
108 // copy to fft input buffer, shifting appropriately
109 for (i=0; i<h_len; i++) {
110 h_prime[i] = h_primef[ (i+k*m)%h_len ];
111 g_prime[i] = g_primef[ (i+k*m)%h_len ];
112 h_tx[i] = ht[ (i+k*m)%h_len ];
113 }
114
115 // run ffts
116 fft_run(h_len, h_prime, H_prime, LIQUID_FFT_FORWARD, 0);
117 fft_run(h_len, g_prime, G_prime, LIQUID_FFT_FORWARD, 0);
118 fft_run(h_len, h_tx, H_tx, LIQUID_FFT_FORWARD, 0);
119
120 #if 0
121 // print results
122 for (i=0; i<h_len; i++)
123 printf("Ht(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(H_tx[i]), cimagf(H_tx[i]));
124 for (i=0; i<h_len; i++)
125 printf("H_prime(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(H_prime[i]), cimagf(H_prime[i]));
126 for (i=0; i<h_len; i++)
127 printf("G_prime(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(G_prime[i]), cimagf(G_prime[i]));
128 #endif
129
130 // find minimum of reponses
131 float H_tx_min = 0.0f;
132 float H_prime_min = 0.0f;
133 float G_prime_min = 0.0f;
134 for (i=0; i<h_len; i++) {
135 if (i==0 || crealf(H_tx[i]) < H_tx_min) H_tx_min = crealf(H_tx[i]);
136 if (i==0 || crealf(H_prime[i]) < H_prime_min) H_prime_min = crealf(H_prime[i]);
137 if (i==0 || crealf(G_prime[i]) < G_prime_min) G_prime_min = crealf(G_prime[i]);
138 }
139
140 // compute 'prototype' response, removing minima, and add correction factor
141 for (i=0; i<h_len; i++) {
142 // compute response necessary to yeild prototype response (not exact, but close)
143 H_hat[i] = crealf(H_prime[i] - H_prime_min + delta) / crealf(H_tx[i] - H_tx_min + delta);
144
145 // include additional term to add stop-band suppression
146 H_hat[i] *= crealf(G_prime[i] - G_prime_min) / crealf(G_prime[0]);
147 }
148
149 // compute ifft and copy response
150 fft_run(h_len, H_hat, h_hat, LIQUID_FFT_BACKWARD, 0);
151 for (i=0; i<h_len; i++)
152 hr[i] = crealf( h_hat[(i+k*m+1)%h_len] ) / (float)(k*h_len);
153
154 //
155 // end of filter design procedure
156 //
157
158 // print results to screen
159 printf("gmsk receive filter:\n");
160 for (i=0; i<h_len; i++)
161 printf(" hr(%3u) = %12.8f;\n", i+1, hr[i]);
162
163 // compute isi
164 float rxy0 = liquid_filter_crosscorr(ht,h_len, hr,h_len, 0);
165 float isi_rms = 0.0f;
166 for (i=1; i<2*m; i++) {
167 float e = liquid_filter_crosscorr(ht,h_len, hr,h_len, i*k) / rxy0;
168 isi_rms += e*e;
169 }
170 isi_rms = sqrtf(isi_rms / (float)(2*m-1));
171 printf("\n");
172 printf("ISI (RMS) = %12.8f dB\n", 20*log10f(isi_rms));
173
174
175 //
176 // export output file
177 //
178 FILE*fid = fopen(OUTPUT_FILENAME,"w");
179 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
180 fprintf(fid,"clear all;\n");
181 fprintf(fid,"close all;\n");
182 fprintf(fid,"\n\n");
183 fprintf(fid,"k = %u;\n", k);
184 fprintf(fid,"m = %u;\n", m);
185 fprintf(fid,"beta = %f;\n", BT);
186 fprintf(fid,"h_len = 2*k*m+1;\n");
187 fprintf(fid,"nfft = 1024;\n");
188 fprintf(fid,"ht = zeros(1,h_len);\n");
189 fprintf(fid,"hp = zeros(1,h_len);\n");
190 fprintf(fid,"hr = zeros(1,h_len);\n");
191
192 // print results
193 for (i=0; i<h_len; i++) fprintf(fid,"ht(%3u) = %12.4e;\n", i+1, ht[i] / k);
194 for (i=0; i<h_len; i++) fprintf(fid,"hr(%3u) = %12.4e;\n", i+1, hr[i] * k);
195 for (i=0; i<h_len; i++) fprintf(fid,"hp(%3u) = %12.4e;\n", i+1, h_primef[i]);
196
197 fprintf(fid,"hc = k*conv(ht,hr);\n");
198
199 // plot results
200 fprintf(fid,"f = [0:(nfft-1)]/nfft - 0.5;\n");
201 fprintf(fid,"Ht = 20*log10(abs(fftshift(fft(ht, nfft))));\n");
202 fprintf(fid,"Hp = 20*log10(abs(fftshift(fft(hp/k,nfft))));\n");
203 fprintf(fid,"Hr = 20*log10(abs(fftshift(fft(hr, nfft))));\n");
204 fprintf(fid,"Hc = 20*log10(abs(fftshift(fft(hc/k,nfft))));\n");
205 fprintf(fid,"\n");
206 fprintf(fid,"figure;\n");
207 fprintf(fid,"plot(f,Ht,'LineWidth',1,'Color',[0.00 0.25 0.50],...\n");
208 fprintf(fid," f,Hp,'LineWidth',1,'Color',[0.80 0.80 0.80],...\n");
209 fprintf(fid," f,Hr,'LineWidth',1,'Color',[0.00 0.50 0.25],...\n");
210 fprintf(fid," f,Hc,'LineWidth',2,'Color',[0.50 0.00 0.00],...\n");
211 fprintf(fid," [-0.5/k 0.5/k], [1 1]*20*log10(0.5),'or');\n");
212 fprintf(fid,"legend('transmit','prototype','receive','composite','alias points',1);\n");
213 fprintf(fid,"xlabel('Normalized Frequency');\n");
214 fprintf(fid,"ylabel('PSD');\n");
215 fprintf(fid,"grid on;\n");
216 fprintf(fid,"axis([-0.5 0.5 -100 20]);\n");
217
218 fprintf(fid,"\n");
219 fprintf(fid,"figure;\n");
220 fprintf(fid,"tr = [ -k*m:k*m]/k;\n");
221 fprintf(fid,"tc = [-2*k*m:2*k*m]/k;\n");
222 fprintf(fid,"ic = [0:k:(4*k*m)]+1;\n");
223 fprintf(fid,"subplot(2,1,1);\n");
224 fprintf(fid," plot(tr,ht,'-x', tr,hr,'-x');\n");
225 fprintf(fid," legend('transmit','receive',1);\n");
226 fprintf(fid," xlabel('Time');\n");
227 fprintf(fid," ylabel('GMSK Tx/Rx Filters');\n");
228 fprintf(fid," grid on;\n");
229 fprintf(fid," axis([-2*m 2*m floor(5*min([hr ht]))/5 ceil(5*max([hr ht]))/5]);\n");
230 fprintf(fid,"subplot(2,1,2);\n");
231 fprintf(fid," plot(tc,hc,'-x', tc(ic),hc(ic),'or');\n");
232 fprintf(fid," xlabel('Time');\n");
233 fprintf(fid," ylabel('GMSK Composite Response');\n");
234 fprintf(fid," grid on;\n");
235 fprintf(fid," axis([-2*m 2*m -0.2 1.2]);\n");
236 fprintf(fid," axis([-2*m 2*m floor(5*min(hc))/5 ceil(5*max(hc))/5]);\n");
237
238 fclose(fid);
239 printf("results written to %s.\n", OUTPUT_FILENAME);
240
241 return 0;
242 }
243
244