1 /*
2 ACfax - Fax reception with X11-interface for amateur radio
3 Copyright (C) 1995-1998 Andreas Czechanowski, DL4SDC
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software Foundation,
17 Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18
19 andreas.czechanowski@ins.uni-stuttgart.de
20 */
21
22 /*
23 * mod_demod.c - FM and AM modulator/demodulator for ACfax
24 */
25
26 #include <string.h>
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <math.h>
30 #include <unistd.h>
31 #include "mod_demod.h"
32
33 #define PI M_PI
34
35 SHORT int firwide[] = { 6, 20, 7, -42, -74, -12, 159, 353, 440 };
36 SHORT int firmiddle[] = { 0, -18, -38, -39, 0, 83, 191, 284, 320 };
37 SHORT int firnarrow[] = { -7, -18, -15, 11, 56, 116, 177, 223, 240 };
38
39 /* these some more coefficient tables were for testing only. */
40 /* int firwide[] = { 12, 24, 7, -42, -74, -12, 159, 353, 440 }; */
41 /* int firnarrow[] = { -16, -21, -15, 11, 56, 116, 177, 223, 240 }; */
42 /* int fircheap[] = { 0, 0, 0, 0, 0, 0, 36, 99, 128 }; */
43 /* int firtoobig[] = { -14, -12, 18, 40, -1, -80, -71, 114, 381, 510 } */
44 /* int firnotbad[] = { 11, 16, 1, -29, -44, -3, 98, 208, 256 }; */
45
46 static SHORT int *fcoeff = firmiddle;
47 static int *sqrt_lo; /* square-root lookup-table for AM demodulator */
48 static int *sqrt_hi; /* (splitted in 2 parts) */
49 static int *sintab; /* sine wave table for FM modulator */
50 static int *asntab; /* arcsin-table for FM-demodulator */
51 static int *fmphsinc; /* phase-increment-table for FM modulator */
52 static char *amreal; /* amplitude tables for AM modulator */
53 static int compval; /* comparator value for APT detecion */
54
55 /*
56 void main(int argc, char **argv)
57 {
58 char indata[2048];
59 int outdata[1024];
60 int q;
61
62 modem_init();
63 set_modem_param (FIL_MIDL, 255, 400);
64 while ((q = fread(indata, 1, 2048, stdin)) > 0) {
65 q &= 0xfffe;
66 fm_demod(indata, q, outdata);
67 }
68 }
69 */
70
71 /* Initialize the modulator/demodulator functions. This function allocates
72 space for the arrays and lookup-tables
73 */
74
modem_init(void)75 void modem_init(void)
76 {
77 int i;
78 static int inited = 0;
79
80 if (inited) return;
81
82 sqrt_lo = (int *)malloc(2048*sizeof(int));
83 sqrt_hi = (int *)malloc(2048*sizeof(int));
84 if (!(sqrt_lo) || !(sqrt_hi)) {
85 perror("modem_init:sqrt_*");
86 exit(1);
87 }
88 amreal = (char *)malloc(258);
89 if (!(amreal)) {
90 perror("modem_init:am{real|imag}");
91 exit(1);
92 }
93 sintab = (int *)malloc(1024*sizeof(int));
94 if (!(sintab)) {
95 perror("modem_init:sintab");
96 exit(1);
97 }
98 asntab = (int *)malloc(512*sizeof(int));
99 if (!(asntab)) {
100 perror("modem_init:asntab");
101 exit(1);
102 }
103 fmphsinc = (int *)malloc(258*sizeof(int));
104 for (i=0; i<1024; i++)
105 sintab[i] = 127.5 + 127*sin(i*PI/512.0);
106
107 inited = -1;
108 }
109
110 /* set up a lookup-table for the given deviation which gives values from
111 0 to maxval with clipping beyond the edges for reception, initialize
112 a lookup-table for FM transmission.
113 */
set_modem_param(int filter,int maxval,int devi)114 void set_modem_param(int filter, int maxval, int devi)
115 {
116 int phmax;
117 int i;
118
119 /* initialize decoder first, if not already done */
120 modem_init();
121 /* do some range-checking */
122 if (devi < 100) devi = 100;
123 if (devi > 1200) devi = 1200;
124 if (maxval < 1) maxval = 1;
125 if (maxval > 255) maxval = 255;
126 phmax = 255.9 * sin(devi * PI / 4000.0) + 0.5;
127 for (i=0; i<512; i++) {
128 if (i <= 256-phmax)
129 asntab[i] = 0;
130 else if (i >= 256+phmax)
131 asntab[i] = maxval;
132 else
133 asntab[i] = maxval * 2000 / devi * asin((i-256.0)/256.5) / PI
134 + maxval/ 2;
135 }
136 for (i=0; i<2048; i++) {
137 sqrt_lo[i] = (sqrt(i) * maxval / 256.0) + 0.5;
138 sqrt_hi[i] = (sqrt(32*i) * maxval / 256.0) + 0.5;
139 }
140 for (i=0; i<=maxval; i++)
141 fmphsinc[i] = devi / 2000.0 * 65536;
142 for (i=0; i<= maxval; i++)
143 amreal[i] = 127.5 + 127.5*i/maxval;
144 compval = maxval >> 1;
145 switch (filter) {
146 case FIL_NARR:
147 fcoeff = firnarrow;
148 fprintf(stderr, "selecting narrow filter\n");
149 break;
150 case FIL_MIDL:
151 fcoeff = firmiddle;
152 fprintf(stderr, "selecting middle filter\n");
153 break;
154 case FIL_WIDE:
155 fcoeff = firwide;
156 fprintf(stderr, "selecting wide filter\n");
157 break;
158 }
159 }
160
161 /* Demodulate the samples taken with 8kHz to a frequency sampled with 4kHz.
162 This is done with a delay-demodulator :
163
164 +----------+ +-------+ +----+ +---+
165 +->|*sin(2kHz)|->|FIR-LPF|--*->|z^-1|->|mul|--+
166 | +----------+ +-------+ | +----+ +---+ |
167 | *-->-- / | AM^2
168 | | \ / |+ |
169 | +---------+ \ / +---+ +---+
170 <input>--+ AM^2<--|amp^2 det| X |add|->|div|--<FM-out>
171 | +---------+ / \ +---+ +---+
172 | | / \ |-
173 | *-->-- \ |
174 | +----------+ +-------+ | +----+ +---+ |
175 +->|*cos(2kHz)|->|FIR-LPF|--*->|z^-1|->|mul|--+
176 +----------+ +-------+ +----+ +---+
177
178 The cosine-signal is simply a sequence of 1,0,-1,0 , and the sine-
179 signal is a sequence of 0,1,0,-1 values because the frequency of
180 these is exactly a fouth of the sample-rate. The values multiplied by
181 zero need not be evaluated by the FIR-filter, what reduces the
182 calculation-effort to one half. Taps of the FIR-chain which are to
183 be multiplied with the same coefficient (possible due to the
184 symmetry of the pulse-response) are first added or subtracted,
185 so again reducing the number of multiplications to one half.
186
187 NOTE: incnt must be a multiple of 2 to operate properly !
188 output range is determined by asntab, which is set up by
189 set_modem_param. This table also performs the arcsin-correction
190 and clips values that are beyond the given deviation.
191 */
fm_demod(char * smplin,int incnt,char * decout)192 void fm_demod(char *smplin, int incnt, char *decout)
193 {
194 static SHORT int firbuf[32]; /* buffer holding the input-values */
195 /* static int inptr = 0; / pointer for data-entry into firbuf */
196 /* static int rdptr = 0; / pointer for data-output from firbuf */
197 static int neg = 0; /* flag changing every 2 sample-points */
198 static int px, py, qx, qy; /* the filtered output-results */
199 static int pamp, qamp, pfrq, qfrq;
200 /*
201 static int n = 0;
202 int smplptr = 0;
203 */
204
205 while (incnt >= 2) {
206
207 /* shift buffer "left" by 2 : copy firbuf[2..17] to firbuf[0..15] */
208 memcpy(firbuf, firbuf+2, (16*sizeof(int)));
209 /* enter 2 new samle-points into buffer */
210 firbuf[16] = *(unsigned char *)smplin++ - 128;
211 firbuf[17] = *(unsigned char *)smplin++ - 128;
212 incnt -= 2;
213
214 /* do the first quarter : multiply with +0-0+0-0+0-0+0-0+ for px
215 and 0-0+0-0+0-0+0-0+0 for py.
216 2 sample-periodes later, multiply with 0-0+0 for px and -0+0- for py,
217 what is done using the neg variable as flag. In this case, the result
218 for frequency is negated (amplitude is always positive) */
219 px = (firbuf[0] + firbuf[16]) * fcoeff[0];
220 px -= (firbuf[2] + firbuf[14]) * fcoeff[2];
221 px += (firbuf[4] + firbuf[12]) * fcoeff[4];
222 px -= (firbuf[6] + firbuf[10]) * fcoeff[6];
223 px += (firbuf[8]) * fcoeff[8];
224 px >>= 2;
225 py = (firbuf[15] - firbuf[1]) * fcoeff[1];
226 py += (firbuf[3] - firbuf[13]) * fcoeff[3];
227 py += (firbuf[11] - firbuf[5]) * fcoeff[5];
228 py += (firbuf[7] - firbuf[9]) * fcoeff[7];
229 py >>= 2;
230 /*
231 dfft->src[smplptr].r = (neg) ? -px : px;
232 dfft->src[smplptr].i = 0;
233 smplptr++;
234 */
235 /* rdptr = (rdptr+1) & 31; */
236 pamp = ((px*px + py*py) >> 8) + 1;
237 pfrq = (px*qy - py*qx) / (pamp + qamp);
238
239 /* do the second quarter : multiply with 0-0+0-0+0-0+0-0+0 for px,
240 -0+0-0+0-0+0-0+0- for py */
241 qx = (firbuf[15+1] - firbuf[1+1]) * fcoeff[1];
242 qx += (firbuf[3+1] - firbuf[13+1]) * fcoeff[3];
243 qx += (firbuf[11+1] - firbuf[5+1]) * fcoeff[5];
244 qx += (firbuf[7+1] - firbuf[9+1]) * fcoeff[7];
245 qx >>= 2;
246 qy = -(firbuf[0+1] + firbuf[16+1]) * fcoeff[0];
247 qy += (firbuf[2+1] + firbuf[14+1]) * fcoeff[2];
248 qy -= (firbuf[4+1] + firbuf[12+1]) * fcoeff[4];
249 qy += (firbuf[6+1] + firbuf[10+1]) * fcoeff[6];
250 qy -= (firbuf[8+1]) * fcoeff[8];
251 qy >>= 2;
252 /*
253 dfft->src[smplptr].r = (neg) ? -qx : qx;
254 dfft->src[smplptr].i = 0;
255 smplptr++;
256 */
257 /* rdptr = (rdptr+1) & 31; */
258 qamp = ((qx*qx + qy*qy) >> 8) + 1;
259 qfrq = (px*qy - py*qx) / (qamp + pamp);
260
261 /* asntab gives values from minval to maxval for given deviation */
262 *decout = (unsigned char)asntab[(pfrq+qfrq+256) & 511];
263 /*
264 printf("%3d %5d %5d\n", n, (px*qy)/1024, (py*qx)/1024);
265 printf("%3d %5d %5d\n", n, *decout, pamp+qamp);
266 n++;
267 */
268 decout++;
269 neg = ~neg;
270 }
271 /*
272 do_dfft(dfft);
273 for (smplptr=0; smplptr<1024; smplptr++) {
274 dx = dfft->dest[smplptr].r;
275 dy = dfft->dest[smplptr].i;
276 printf("%4.1f %3.1f\n", smplptr * 3.906, 20*log10(sqrt(dx*dx + dy*dy)));
277 }
278 */
279 }
280
281
282 /* Demodulate the samples taken with 9.6kHz to an amplitude sampled with 4.8kHz.
283 This done with a synchronous demodulator :
284
285 +------------+ +-------+ +------+
286 +->|*sin(2.4kHz)|->|FIR-LPF|->|square|-+
287 | +------------+ +-------+ +------+ |
288 | |
289 | |
290 | +-----+ +-------+
291 <input>--+ |adder|->|sq.root|--<AM-out>
292 | +-----+ +-------+
293 | |
294 | |
295 | +------------+ +-------+ +------+ |
296 +->|*cos(2.4kHz)|->|FIR-LPF|->|square|-+
297 +------------+ +-------+ +------+
298
299 The cosine-signal is simply a sequence of 1,0,-1,0 , and the sine-
300 signal is a sequence of 0,1,0,-1 values because the frequency of
301 these is exactly a fouth of the sample-rate. The values multiplied by
302 zero need not be evaluated by the FIR-filter, what reduces the
303 calculation-effort to one half.
304
305 NOTE: incnt must be a multiple of 2 to operate properly !
306 output range is determined by sqrt_lo and sqrt_hi, which
307 are set up by set_modem_param.
308 */
am_demod(char * smplin,int incnt,char * decout)309 void am_demod(char *smplin, int incnt, char *decout)
310 {
311 static SHORT int firbuf[32]; /* buffer holding the input-values */
312 /* static int inptr = 0; / pointer for data-entry into firbuf */
313 /* static int rdptr = 0; / pointer for data-output from firbuf */
314 static int neg = 0; /* flag changing every 2 sample-points */
315 static int px, py, qx, qy; /* the filtered output-results */
316 static int pamp, qamp;
317 /*
318 static int n = 0;
319 int smplptr = 0;
320 */
321
322 while (incnt >= 2) {
323
324 /* shift buffer "left" by 2 : copy firbuf[2..17] to firbuf[0..15] */
325 memcpy(firbuf, firbuf+2, (16*sizeof(int)));
326 /* enter 2 new samle-points into buffer */
327 firbuf[16] = *(unsigned char *)smplin++ - 128;
328 firbuf[17] = *(unsigned char *)smplin++ - 128;
329 incnt -= 2;
330
331 /* do the first quarter : multiply with +0-0+0-0+0-0+0-0+ for px
332 and 0-0+0-0+0-0+0-0+0 for py.
333 2 sample-periodes later, multiply with 0-0+0 for px and -0+0- for py,
334 what is done using the neg variable as flag. In this case, the result
335 for frequency is negated (amplitude is always positive) */
336 px = (firbuf[0] + firbuf[16]) * fcoeff[0];
337 px -= (firbuf[2] + firbuf[14]) * fcoeff[2];
338 px += (firbuf[4] + firbuf[12]) * fcoeff[4];
339 px -= (firbuf[6] + firbuf[10]) * fcoeff[6];
340 px += firbuf[8] * fcoeff[8];
341 px >>= 2;
342 py = (firbuf[15] - firbuf[1]) * fcoeff[1];
343 py += (firbuf[3] - firbuf[13]) * fcoeff[3];
344 py += (firbuf[11] - firbuf[5]) * fcoeff[5];
345 py += (firbuf[7] - firbuf[9]) * fcoeff[7];
346 py >>= 2;
347 /*
348 dfft->src[smplptr].r = (neg) ? -px : px;
349 dfft->src[smplptr].i = 0;
350 smplptr++;
351 */
352 /* rdptr = (rdptr+1) & 31; */
353 pamp = (px*px + py*py);
354
355 /* do the second quarter : multiply with 0-0+0-0+0-0+0-0+0 for px,
356 -0+0-0+0-0+0-0+0- for py */
357 qx = (firbuf[15+1] - firbuf[1+1]) * fcoeff[1];
358 qx += (firbuf[3+1] - firbuf[13+1]) * fcoeff[3];
359 qx += (firbuf[11+1] - firbuf[5+1]) * fcoeff[5];
360 qx += (firbuf[7+1] - firbuf[9+1]) * fcoeff[7];
361 qx >>= 2;
362 qy = -(firbuf[0+1] + firbuf[16+1]) * fcoeff[0];
363 qy += (firbuf[2+1] + firbuf[14+1]) * fcoeff[2];
364 qy -= (firbuf[4+1] + firbuf[12+1]) * fcoeff[4];
365 qy += (firbuf[6+1] + firbuf[10+1]) * fcoeff[6];
366 qy -= firbuf[8+1] * fcoeff[8];
367 qy >>= 2;
368 /*
369 dfft->src[smplptr].r = (neg) ? -qx : qx;
370 dfft->src[smplptr].i = 0;
371 smplptr++;
372 */
373 /* rdptr = (rdptr+1) & 31; */
374 qamp = (qx*qx + qy*qy + pamp) / 6400;
375
376 if (qamp >= 65535)
377 *decout = sqrt_hi[2047];
378 else if (qamp >= 2047)
379 *decout = (unsigned char)sqrt_hi[qamp >> 5];
380 else
381 *decout = (unsigned char)sqrt_lo[qamp];
382
383 /*
384 printf("%3d %5d %5d\n", n, (px*px)/1024, (py*py)/1024);
385 printf("%3d %5d %5d\n", n, *decout, qamp);
386 n++;
387 */
388 decout++;
389 neg = ~neg;
390 }
391 }
392
393 /* Encode a frequency sampled at 4kHz to an AF-signal sampled at 8kHz.
394 This is done by shifting the phase by the appropriate value for each
395 sample.
396 */
fm_modulate(char * codin,int incnt,char * smplout)397 void fm_modulate(char *codin, int incnt, char *smplout)
398 {
399 static int phs = 0;
400
401 while (incnt-- > 0) {
402 phs += 0x10000 + fmphsinc[*(unsigned char *)codin++];
403 phs &= 0x3ffff;
404 *smplout++ = (char)sintab[phs >> 8];
405 phs += 0x10000 + fmphsinc[*(unsigned char *)codin++];
406 phs &= 0x3ffff;
407 *smplout++ = (char)sintab[phs >> 8];
408 }
409 }
410
411 /* Encode an amplitude sampled at 4.8kHz to an AF-signal sampled at 9.6kHz.
412 This is done by multiplying the amplitude-values with a sampled sine-
413 function represented by a 0,1,0,-1 sequence.
414 */
am_modulate(char * codin,int incnt,char * smplout)415 void am_modulate(char *codin, int incnt, char *smplout)
416 {
417 static int neg = 0;
418
419 while (incnt-- > 0) {
420 if (neg) {
421 *smplout++ = amreal[*(unsigned char *)codin++];
422 *smplout++ = 128; /*amimag[*codin++]*/
423 } else {
424 *smplout++ = amreal[*(unsigned char *)codin++] ^ 0xff;
425 *smplout++ = 128; /*amimag[*codin++] ^ 0xff*/
426 }
427 neg = ~neg;
428 }
429 }
430
431