1 // ----------------------------------------------------------------------------
2 //
3 // filters.cxx -- Several Digital Filter classes used in fldigi
4 //
5 // Copyright (C) 2006-2008 Dave Freese, W1HKJ
6 //
7 // These filters are based on the gmfsk design and the design notes given in
8 // "Digital Signal Processing, A Practical Guid for Engineers and Scientists"
9 // by Steven W. Smith.
10 //
11 // This file is part of fldigi.
12 //
13 // Fldigi is free software: you can redistribute it and/or modify
14 // it under the terms of the GNU General Public License as published by
15 // the Free Software Foundation, either version 3 of the License, or
16 // (at your option) any later version.
17 //
18 // Fldigi is distributed in the hope that it will be useful,
19 // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 // GNU General Public License for more details.
22 //
23 // You should have received a copy of the GNU General Public License
24 // along with fldigi. If not, see <http://www.gnu.org/licenses/>.
25 // ----------------------------------------------------------------------------
26
27 #include <config.h>
28
29 #include <stdlib.h>
30 #include <stdio.h>
31 #include <string.h>
32
33 #include "filters.h"
34
35 #include <iostream>
36
37
38 //=====================================================================
39 // C_FIR_filter
40 //
41 // a class of Finite Impulse Response (FIR) filters with
42 // decimate in time capability
43 //
44 //=====================================================================
45
C_FIR_filter()46 C_FIR_filter::C_FIR_filter () {
47 pointer = counter = length = 0;
48 decimateratio = 1;
49 ifilter = qfilter = (double *)0;
50 ffreq = 0.0;
51 }
52
~C_FIR_filter()53 C_FIR_filter::~C_FIR_filter() {
54 if (ifilter) delete [] ifilter;
55 if (qfilter) delete [] qfilter;
56 }
57
init(int len,int dec,double * itaps,double * qtaps)58 void C_FIR_filter::init(int len, int dec, double *itaps, double *qtaps) {
59 length = len;
60 decimateratio = dec;
61 if (ifilter) {
62 delete [] ifilter;
63 ifilter = (double *)0;
64 }
65 if (qfilter) {
66 delete [] qfilter;
67 qfilter = (double *)0;
68 }
69
70 for (int i = 0; i < FIRBufferLen; i++)
71 ibuffer[i] = qbuffer[i] = 0.0;
72
73 if (itaps) {
74 ifilter = new double[len];
75 for (int i = 0; i < len; i++) ifilter[i] = itaps[i];
76 }
77 if (qtaps) {
78 qfilter = new double[len];
79 for (int i = 0; i < len; i++) qfilter[i] = qtaps[i];
80 }
81
82 pointer = len;
83 counter = 0;
84 }
85
86
87 //=====================================================================
88 // Create a band pass FIR filter with 6 dB corner frequencies
89 // of 'f1' and 'f2'. (0 <= f1 < f2 <= 0.5)
90 //=====================================================================
91
bp_FIR(int len,int hilbert,double f1,double f2)92 double * C_FIR_filter::bp_FIR(int len, int hilbert, double f1, double f2)
93 {
94 double *fir;
95 double t, h, x;
96
97 fir = new double[len];
98
99 for (int i = 0; i < len; i++) {
100 t = i - (len - 1.0) / 2.0;
101 h = i * (1.0 / (len - 1.0));
102
103 if (!hilbert) {
104 x = (2 * f2 * sinc(2 * f2 * t) -
105 2 * f1 * sinc(2 * f1 * t)) * hamming(h);
106 } else {
107 x = (2 * f2 * cosc(2 * f2 * t) -
108 2 * f1 * cosc(2 * f1 * t)) * hamming(h);
109 // The actual filter code assumes the impulse response
110 // is in time reversed order. This will be anti-
111 // symmetric so the minus sign handles that for us.
112 x = -x;
113 }
114
115 fir[i] = x;
116 }
117
118 return fir;
119 }
120
121 //=====================================================================
122 // Filter will be a lowpass with
123 // length = len
124 // decimation = dec
125 // 0.5 frequency point = freq
126 //=====================================================================
127
init_lowpass(int len,int dec,double freq)128 void C_FIR_filter::init_lowpass (int len, int dec, double freq) {
129 double *fi = bp_FIR(len, 0, 0.0, freq);
130 ffreq = freq;
131 init (len, dec, fi, fi);
132 delete [] fi;
133 }
134
135 //=====================================================================
136 // Filter will be a bandpass with
137 // length = len
138 // decimation = dec
139 // 0.5 frequency points of f1 (low) and f2 (high)
140 //=====================================================================
141
init_bandpass(int len,int dec,double f1,double f2)142 void C_FIR_filter::init_bandpass (int len, int dec, double f1, double f2) {
143 double *fi = bp_FIR (len, 0, f1, f2);
144 init (len, dec, fi, fi);
145 delete [] fi;
146 }
147
148 //=====================================================================
149 // Filter will the Hilbert form
150 //=====================================================================
151
init_hilbert(int len,int dec)152 void C_FIR_filter::init_hilbert (int len, int dec) {
153 double *fi = bp_FIR(len, 0, 0.05, 0.45);
154 double *fq = bp_FIR(len, 1, 0.05, 0.45);
155 init (len, dec, fi, fq);
156 delete [] fi;
157 delete [] fq;
158 }
159
160
161 //=====================================================================
162 // Run
163 // passes a cmplx value (in) and receives the cmplx value (out)
164 // function returns 0 if the filter is not yet stable
165 // returns 1 when stable and decimated cmplx output value is valid
166 //=====================================================================
167
run(const cmplx & in,cmplx & out)168 int C_FIR_filter::run (const cmplx &in, cmplx &out) {
169 ibuffer[pointer] = in.real();
170 qbuffer[pointer] = in.imag();
171 counter++;
172 if (counter == decimateratio)
173 out = cmplx ( mac(&ibuffer[pointer - length], ifilter, length),
174 mac(&qbuffer[pointer - length], qfilter, length) );
175 pointer++;
176 if (pointer == FIRBufferLen) {
177 /// memmove is necessary if length >= FIRBufferLen/2 , theoretically possible.
178 memmove (ibuffer, ibuffer + FIRBufferLen - length, length * sizeof (double) );
179 memmove (qbuffer, qbuffer + FIRBufferLen - length, length * sizeof (double) );
180 pointer = length;
181 }
182 if (counter == decimateratio) {
183 counter = 0;
184 return 1;
185 }
186 return 0;
187 }
188
189 //=====================================================================
190 // Run the filter for the Real part of the cmplx variable
191 //=====================================================================
192
Irun(const double & in,double & out)193 int C_FIR_filter::Irun (const double &in, double &out) {
194 double *iptr = ibuffer + pointer;
195
196 pointer++;
197 counter++;
198
199 *iptr = in;
200
201 if (counter == decimateratio) {
202 out = mac(iptr - length, ifilter, length);
203 }
204
205 if (pointer == FIRBufferLen) {
206 iptr = ibuffer + FIRBufferLen - length;
207 memcpy(ibuffer, iptr, length * sizeof(double));
208 pointer = length;
209 }
210
211 if (counter == decimateratio) {
212 counter = 0;
213 return 1;
214 }
215
216 return 0;
217 }
218
219 //=====================================================================
220 // Run the filter for the Imaginary part of the cmplx variable
221 //=====================================================================
222
Qrun(const double & in,double & out)223 int C_FIR_filter::Qrun (const double &in, double &out) {
224 double *qptr = ibuffer + pointer;
225
226 pointer++;
227 counter++;
228
229 *qptr = in;
230
231 if (counter == decimateratio) {
232 out = mac(qptr - length, qfilter, length);
233 }
234
235 if (pointer == FIRBufferLen) {
236 qptr = qbuffer + FIRBufferLen - length;
237 memcpy(qbuffer, qptr, length * sizeof(double));
238 pointer = length;
239 }
240
241 if (counter == decimateratio) {
242 counter = 0;
243 return 1;
244 }
245
246 return 0;
247 }
248
249
250 //=====================================================================
251 // Moving average filter
252 //
253 // Simple in concept, sublime in implementation ... the fastest filter
254 // in the west. Also optimal for the processing of time domain signals
255 // characterized by a transition edge. The is the perfect signal filter
256 // for CW, RTTY and other signals of that type. For a given filter size
257 // it provides the greatest s/n improvement while retaining the sharpest
258 // leading edge on the filtered signal.
259 //=====================================================================
260
Cmovavg(int filtlen)261 Cmovavg::Cmovavg (int filtlen)
262 {
263 len = filtlen;
264 in = new double[len];
265 empty = true;
266 }
267
~Cmovavg()268 Cmovavg::~Cmovavg()
269 {
270 if (in) delete [] in;
271 }
272
run(double a)273 double Cmovavg::run(double a)
274 {
275 if (!in) {
276 return a;
277 }
278 if (empty) {
279 empty = false;
280 out = 0;
281 for (int i = 0; i < len; i++) {
282 in[i] = a;
283 out += a;
284 }
285 pint = 0;
286 return a;
287 }
288 out = out - in[pint] + a;
289 in[pint] = a;
290 if (++pint >= len) pint = 0;
291 return out / len;
292 }
293
setLength(int filtlen)294 void Cmovavg::setLength(int filtlen)
295 {
296 if (filtlen > len) {
297 if (in) delete [] in;
298 in = new double[filtlen];
299 }
300 len = filtlen;
301 empty = true;
302 }
303
reset()304 void Cmovavg::reset()
305 {
306 empty = true;
307 }
308
309 //=====================================================================
310 // Sliding FFT filter
311 // Sliding Fast Fourier Transform
312 //
313 // The sliding FFT ingeniously exploits the properties of a time-delayed
314 // input and the property of linearity for its derivation.
315 //
316 // First of all, the N-point transform of a sequence x(n) is equal to the
317 // summation of the transforms of N separate transforms where each transform
318 // has just one of the original samples at it's original sample time.
319 //
320 //i.e.
321 // transform of [x0, x1, x2, x3, x4, x5,...xN-1]
322 // is equal to
323 // transform of [x0, 0, 0, 0, 0, 0,...0]
324 // + transform of [0, x1, 0, 0, 0, 0,...0]
325 // + transform of [0, 0, x2, 0, 0, 0,...0]
326 // + transform of [0, 0, 0, x3, 0, 0,...0]
327 // .
328 // .
329 // .
330 // + transform of [0, 0, 0, 0, 0, 0,...xN-1]
331 //
332 // Secondly, the transform of a time-delayed sequence is a phase-rotated
333 // version of the transform of the original sequence. i.e.
334 //
335 // If x(n) transforms to X(k),
336 // Then x(n-m) transforms to X(k)(Wn)^(-mk)
337 //
338 // where N is the FFT size, m is the delay in sample periods, and WN is the
339 // familiar phase-rotating coefficient or twiddle factor e^(-j2p/N)
340 //
341 // Therefore, if the N-point transform X(k) of an individual sample is considered,
342 // and then the sample is moved back in time by one sample period, all frequency
343 // bins of X(k) are phase-rotated by 2pk/N radians.
344 //
345 // The important thing here is that the transform is not performed again because
346 // the previous frequency results can be used by simply application of the correct
347 // coefficients.
348 //
349 // This is the technique that is applied when the rectangular sampling window
350 // slides along by one sample. The contributions of all samples that are
351 // included in both the original and the new windows are simply phase rotated.
352 // The end effects are that the transform of the new sample must be added, and
353 // the transform of the oldest sample that disappeared off the end must be
354 // subtracted. These end-effects are easy to perform if we treat the new sample
355 // as occurring at time t = 0, because the transform of a single sample at t = 0,
356 // say (a + bj), simply has all frequency bins equal to (a + bj). Similarly, the
357 // oldest sample that has just disappeared off the end of the window is exactly N
358 // samples old. I.e. it occurred at t = -N. The transform of this sample,
359 // say (c + dj), is also straightforward since every frequency bin has now been
360 // phase-rotated an integer number of times from when the sample was at t = 0.
361 // (The kth frequency bin has been rotated by 2pk radians). The transform of the
362 // sample at t = -N is therefore the same as if it was still at t = 0. I.e. it
363 // has all frequency bins equal to (c + dj).
364 //
365 // All that is needed therefore is to
366 // phase rotate each frequency bin in F(k) by WN^(k) and then
367 // add [(a + bj) + (c + dj)] to each frequency bin.
368 //
369 // One cmplx multiplication and two cmplx additions per frequency bin are
370 // therefore required, per sample period, regardless of the size of the transform.
371 //
372 // For example, a traditional 1024-point FFT needs 5120 cmplx multiplies
373 // and 10240 cmplx additions to calculate all 1024 frequency bins. A 1024-point
374 // Sliding FFT however needs 1024 cmplx multiplies and 2048 cmplx additions
375 // for all 1024 frequency bins, and as each frequency bin is calculated separately,
376 // it is only necessary to calculate the ones that are of interest.
377 //
378 // One drawback of the Sliding FFT is that in using feedback from previous
379 // frequency bins, there is potential for instability if the coefficients are not
380 // infinitely precise. Without infinite precision, stability can be guaranteed by
381 // making each phase-rotation coefficient have a magnitude of slightly less than
382 // unity. E.g. 0.9999.
383 //
384 // This then has to taken into account when the Nth sample is subtracted, because
385 // the factor 0.9999 has been applied N times to the transform of this sample.
386 // The sample cannot therefore be directly subtracted, it must first be multiplied
387 // by the factor of 0.9999^N. This unfortunately means there is another multipli-
388 // cation to perform per frequency bin. Another drawback is that a circular buffer
389 // is needed in which to keep N samples, so that the oldest sample, (from t= -N),
390 // can be subtracted each time.
391 //
392 // This filter is ideal for extracting a finite number of frequency bins
393 // with a very long kernel length. The filter only needs to calculate the
394 // values for the bins of interest and not the entire spectrum. It does
395 // require the store of the history associated with those bins over the
396 // kernel length.
397 //
398 // Use in the MFSK / DOMINO modem for extraction of the frequency spectra
399 //
400 //=====================================================================
401
402 struct sfft::vrot_bins_pair {
403 cmplx vrot;
404 cmplx bins;
405 } ;
406
sfft(int len,int _first,int _last)407 sfft::sfft(int len, int _first, int _last)
408 {
409 vrot_bins = new vrot_bins_pair[len];
410 delay = new cmplx[len];
411 fftlen = len;
412 first = _first;
413 last = _last;
414 ptr = 0;
415 double phi = 0.0, tau = 2.0 * M_PI/ len;
416 k2 = 1.0;
417 for (int i = 0; i < fftlen; i++) {
418 vrot_bins[i].vrot = cmplx( K1 * cos (phi), K1 * sin (phi) );
419 phi += tau;
420 delay[i] = vrot_bins[i].bins = 0.0;
421 k2 *= K1;
422 }
423 count = 0;
424 }
425
~sfft()426 sfft::~sfft()
427 {
428 delete [] vrot_bins;
429 delete [] delay;
430 }
431
reset()432 void sfft::reset()
433 {
434 for (int i = 0; i < fftlen; i++) delay[i] = vrot_bins[i].bins = 0.0;
435 count = 0;
436 }
437
is_stable()438 bool sfft::is_stable()
439 {
440 return (count >= fftlen);
441 }
442
443 // Sliding FFT, cmplx input, cmplx output
444 // FFT is computed for each value from first to last
445 // Values are not stable until more than "len" samples have been processed.
446 // Copies the frequencies to a pointer with a given stride.
run(const cmplx & input,cmplx * __restrict__ result,int stride)447 void sfft::run(const cmplx& input, cmplx * __restrict__ result, int stride )
448 {
449 cmplx & de = delay[ptr];
450 const cmplx z( input.real() - k2 * de.real(), input.imag() - k2 * de.imag());
451 de = input;
452
453 ++ptr ;
454 if( ptr >= fftlen ) ptr = 0 ;
455
456 // It is more efficient to have vrot and bins very close to each other.
457 for( vrot_bins_pair
458 * __restrict__ itr = vrot_bins + first,
459 * __restrict__ end = vrot_bins + last ;
460 itr != end ;
461 ++itr, result += stride ) {
462 *result = itr->bins = itr->bins * itr->vrot + z * itr->vrot;
463 }
464 if (count < fftlen) count++;
465 }
466
467 // ============================================================================
468 // Goertzel filter
469 // Optimized implementation of a DFT for a single frequency of interest
470 // SR = sample rate
471 // N = Block size (does not need to be a factor of 2!)
472 // bin size = SR / N
473 // K = frequency bin of interest = (N * freq / SR)
474 // N should be selected to make K an integer if possible
475 //
476 // Q0 = current sample
477 // Q1 = previous sample (1 delay)
478 // Q2 = previous sample (2 delay)
479 // w = (2 * pi * K / N)
480 // k1 = cos(w)
481 // k2 = sin(w)
482 // k3 = 2.0 * k1
483
484 // Q0, Q1, Q2 are initialized to zero
485 // Iterate N times:
486 // Q0 = k3*Q1 - Q2 + sample
487 // Q2 = Q1
488 // Q1 = Q0
489 //
490 // After N interations:
491 // real = (Q1 - Q2 * k1)
492 // imag = Q2 * k2
493 // or
494 // mag = Q1*Q1 + Q2*Q2 - Q1*Q2*k1
495 // ============================================================================
496
goertzel(int n,double freq,double sr)497 goertzel::goertzel(int n, double freq, double sr)
498 {
499 double w;
500 w = 2 * M_PI * freq / sr;
501 k1 = cos(w);
502 k2 = sin(w);
503 k3 = 2.0 * k1;
504 Q0 = Q1 = Q2 = 0.0;
505 count = N = n;
506 }
507
~goertzel()508 goertzel::~goertzel()
509 {
510 }
511
reset()512 void goertzel::reset()
513 {
514 Q0 = Q1 = Q2 = 0.0;
515 count = N;
516 }
517
reset(int n,double freq,double sr)518 void goertzel::reset(int n, double freq, double sr)
519 {
520 double w;
521 w = 2 * M_PI * freq / sr;
522 k1 = cos(w);
523 k2 = sin(w);
524 k3 = 2.0 * k1;
525 Q0 = Q1 = Q2 = 0.0;
526 count = N = n;
527 }
528
run(double sample)529 bool goertzel::run(double sample)
530 {
531 Q0 = sample + k3*Q1 - Q2;
532 Q2 = Q1;
533 Q1 = Q0;
534 if (count) { --count; return false; }
535 return true;
536 }
537
real()538 double goertzel::real()
539 {
540 return ((0.5*k3*Q1 - Q2)/N);
541 }
542
imag()543 double goertzel::imag()
544 {
545 return ((k2*Q1)/N);
546 }
547
mag()548 double goertzel::mag()
549 {
550 return (Q2*Q2 + Q1*Q1 - k3*Q2*Q1);
551 }
552