SDRAngel  4.11.5
Developer docs for <a href="https://github.com/f4exb/sdrangel">SDRangel<\a>, an Open Source Qt5 / OpenGL 3.0+ SDR and signal analyzer frontend to various hardware.
fftfilt.cpp
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // fftfilt.cxx -- Fast convolution Overlap-Add filter
3 //
4 // Filter implemented using overlap-add FFT convolution method
5 // h(t) characterized by Windowed-Sinc impulse response
6 //
7 // Reference:
8 // "The Scientist and Engineer's Guide to Digital Signal Processing"
9 // by Dr. Steven W. Smith, http://www.dspguide.com
10 // Chapters 16, 18 and 21
11 //
12 // Copyright (C) 2006-2008 Dave Freese, W1HKJ
13 //
14 // This file is part of fldigi.
15 //
16 // Fldigi is free software: you can redistribute it and/or modify
17 // it under the terms of the GNU General Public License as published by
18 // the Free Software Foundation, either version 3 of the License, or
19 // (at your option) any later version.
20 //
21 // Fldigi is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with fldigi. If not, see <http://www.gnu.org/licenses/>.
28 // ----------------------------------------------------------------------------
29 
30 #include <memory.h>
31 #include <algorithm>
32 #include <iostream>
33 #include <fstream>
34 #include <cstdlib>
35 #include <cmath>
36 #include <typeinfo>
37 
38 #include <stdio.h>
39 #include <sys/types.h>
40 #include <memory.h>
41 
42 #include <dsp/misc.h>
43 #include <dsp/fftfilt.h>
44 
45 //------------------------------------------------------------------------------
46 // initialize the filter
47 // create forward and reverse FFTs
48 //------------------------------------------------------------------------------
49 
50 // Only need a single instance of g_fft, used for both forward and reverse
52 {
53  flen2 = flen >> 1;
54  fft = new g_fft<float>(flen);
55 
56  filter = new cmplx[flen];
57  filterOpp = new cmplx[flen];
58  data = new cmplx[flen];
59  output = new cmplx[flen2];
60  ovlbuf = new cmplx[flen2];
61 
62  memset(filter, 0, flen * sizeof(cmplx));
63  memset(filterOpp, 0, flen * sizeof(cmplx));
64  memset(data, 0, flen * sizeof(cmplx));
65  memset(output, 0, flen2 * sizeof(cmplx));
66  memset(ovlbuf, 0, flen2 * sizeof(cmplx));
67 
68  inptr = 0;
69 }
70 
71 //------------------------------------------------------------------------------
72 // fft filter
73 // f1 < f2 ==> band pass filter
74 // f1 > f2 ==> band reject filter
75 // f1 == 0 ==> low pass filter
76 // f2 == 0 ==> high pass filter
77 //------------------------------------------------------------------------------
78 fftfilt::fftfilt(float f1, float f2, int len)
79 {
80  flen = len;
81  pass = 0;
82  window = 0;
83  init_filter();
84  create_filter(f1, f2);
85 }
86 
87 fftfilt::fftfilt(float f2, int len)
88 {
89  flen = len;
90  pass = 0;
91  window = 0;
92  init_filter();
94 }
95 
97 {
98  if (fft) delete fft;
99 
100  if (filter) delete [] filter;
101  if (filterOpp) delete [] filterOpp;
102  if (data) delete [] data;
103  if (output) delete [] output;
104  if (ovlbuf) delete [] ovlbuf;
105 }
106 
107 void fftfilt::create_filter(float f1, float f2)
108 {
109  // initialize the filter to zero
110  memset(filter, 0, flen * sizeof(cmplx));
111 
112  // create the filter shape coefficients by fft
113  bool b_lowpass, b_highpass;
114  b_lowpass = (f2 != 0);
115  b_highpass = (f1 != 0);
116 
117  for (int i = 0; i < flen2; i++) {
118  filter[i] = 0;
119  // lowpass @ f2
120  if (b_lowpass)
121  filter[i] += fsinc(f2, i, flen2);
122  // highighpass @ f1
123  if (b_highpass)
124  filter[i] -= fsinc(f1, i, flen2);
125  }
126  // highpass is delta[flen2/2] - h(t)
127  if (b_highpass && f2 < f1)
128  filter[flen2 / 2] += 1;
129 
130  for (int i = 0; i < flen2; i++)
131  filter[i] *= _blackman(i, flen2);
132 
133  fft->ComplexFFT(filter); // filter was expressed in the time domain (impulse response)
134 
135  // normalize the output filter for unity gain
136  float scale = 0, mag;
137  for (int i = 0; i < flen2; i++) {
138  mag = abs(filter[i]);
139  if (mag > scale) scale = mag;
140  }
141  if (scale != 0) {
142  for (int i = 0; i < flen; i++)
143  filter[i] /= scale;
144  }
145 }
146 
147 // Double the size of FFT used for equivalent SSB filter or assume FFT is half the size of the one used for SSB
149 {
150  // initialize the filter to zero
151  memset(filter, 0, flen * sizeof(cmplx));
152 
153  for (int i = 0; i < flen2; i++) {
154  filter[i] = fsinc(f2, i, flen2);
155  filter[i] *= _blackman(i, flen2);
156  }
157 
158  fft->ComplexFFT(filter); // filter was expressed in the time domain (impulse response)
159 
160  // normalize the output filter for unity gain
161  float scale = 0, mag;
162  for (int i = 0; i < flen2; i++) {
163  mag = abs(filter[i]);
164  if (mag > scale) scale = mag;
165  }
166  if (scale != 0) {
167  for (int i = 0; i < flen; i++)
168  filter[i] /= scale;
169  }
170 }
171 
172 // Double the size of FFT used for equivalent SSB filter or assume FFT is half the size of the one used for SSB
173 // used with runAsym for in band / opposite band asymmetrical filtering. Can be used for vestigial sideband modulation.
174 void fftfilt::create_asym_filter(float fopp, float fin)
175 {
176  // in band
177  // initialize the filter to zero
178  memset(filter, 0, flen * sizeof(cmplx));
179 
180  for (int i = 0; i < flen2; i++) {
181  filter[i] = fsinc(fin, i, flen2);
182  filter[i] *= _blackman(i, flen2);
183  }
184 
185  fft->ComplexFFT(filter); // filter was expressed in the time domain (impulse response)
186 
187  // normalize the output filter for unity gain
188  float scale = 0, mag;
189  for (int i = 0; i < flen2; i++) {
190  mag = abs(filter[i]);
191  if (mag > scale) scale = mag;
192  }
193  if (scale != 0) {
194  for (int i = 0; i < flen; i++)
195  filter[i] /= scale;
196  }
197 
198  // opposite band
199  // initialize the filter to zero
200  memset(filterOpp, 0, flen * sizeof(cmplx));
201 
202  for (int i = 0; i < flen2; i++) {
203  filterOpp[i] = fsinc(fopp, i, flen2);
204  filterOpp[i] *= _blackman(i, flen2);
205  }
206 
207  fft->ComplexFFT(filterOpp); // filter was expressed in the time domain (impulse response)
208 
209  // normalize the output filter for unity gain
210  scale = 0;
211  for (int i = 0; i < flen2; i++) {
212  mag = abs(filterOpp[i]);
213  if (mag > scale) scale = mag;
214  }
215  if (scale != 0) {
216  for (int i = 0; i < flen; i++)
217  filterOpp[i] /= scale;
218  }
219 }
220 
221 // This filter is constructed directly from frequency domain response. Run with runFilt.
222 void fftfilt::create_rrc_filter(float fb, float a)
223 {
224  std::fill(filter, filter+flen, 0);
225 
226  for (int i = 0; i < flen; i++) {
227  filter[i] = frrc(fb, a, i, flen);
228  }
229 
230  // normalize the output filter for unity gain
231  float scale = 0, mag;
232  for (int i = 0; i < flen; i++)
233  {
234  mag = abs(filter[i]);
235  if (mag > scale) {
236  scale = mag;
237  }
238  }
239  if (scale != 0)
240  {
241  for (int i = 0; i < flen; i++) {
242  filter[i] /= scale;
243  }
244  }
245 }
246 
247 // test bypass
248 int fftfilt::noFilt(const cmplx & in, cmplx **out)
249 {
250  data[inptr++] = in;
251  if (inptr < flen2)
252  return 0;
253  inptr = 0;
254 
255  *out = data;
256  return flen2;
257 }
258 
259 // Filter with fast convolution (overlap-add algorithm).
260 int fftfilt::runFilt(const cmplx & in, cmplx **out)
261 {
262  data[inptr++] = in;
263  if (inptr < flen2)
264  return 0;
265  inptr = 0;
266 
267  fft->ComplexFFT(data);
268  for (int i = 0; i < flen; i++)
269  data[i] *= filter[i];
270 
272 
273  for (int i = 0; i < flen2; i++) {
274  output[i] = ovlbuf[i] + data[i];
275  ovlbuf[i] = data[flen2 + i];
276  }
277  memset (data, 0, flen * sizeof(cmplx));
278 
279  *out = output;
280  return flen2;
281 }
282 
283 // Second version for single sideband
284 int fftfilt::runSSB(const cmplx & in, cmplx **out, bool usb, bool getDC)
285 {
286  data[inptr++] = in;
287  if (inptr < flen2)
288  return 0;
289  inptr = 0;
290 
291  fft->ComplexFFT(data);
292 
293  // get or reject DC component
294  data[0] = getDC ? data[0]*filter[0] : 0;
295 
296  // Discard frequencies for ssb
297  if (usb)
298  {
299  for (int i = 1; i < flen2; i++) {
300  data[i] *= filter[i];
301  data[flen2 + i] = 0;
302  }
303  }
304  else
305  {
306  for (int i = 1; i < flen2; i++) {
307  data[i] = 0;
308  data[flen2 + i] *= filter[flen2 + i];
309  }
310  }
311 
312  // in-place FFT: freqdata overwritten with filtered timedata
314 
315  // overlap and add
316  for (int i = 0; i < flen2; i++) {
317  output[i] = ovlbuf[i] + data[i];
318  ovlbuf[i] = data[i+flen2];
319  }
320  memset (data, 0, flen * sizeof(cmplx));
321 
322  *out = output;
323  return flen2;
324 }
325 
326 // Version for double sideband. You have to double the FFT size used for SSB.
327 int fftfilt::runDSB(const cmplx & in, cmplx **out, bool getDC)
328 {
329  data[inptr++] = in;
330  if (inptr < flen2)
331  return 0;
332  inptr = 0;
333 
334  fft->ComplexFFT(data);
335 
336  for (int i = 0; i < flen2; i++) {
337  data[i] *= filter[i];
338  data[flen2 + i] *= filter[flen2 + i];
339  }
340 
341  // get or reject DC component
342  data[0] = getDC ? data[0] : 0;
343 
344  // in-place FFT: freqdata overwritten with filtered timedata
346 
347  // overlap and add
348  for (int i = 0; i < flen2; i++) {
349  output[i] = ovlbuf[i] + data[i];
350  ovlbuf[i] = data[i+flen2];
351  }
352 
353  memset (data, 0, flen * sizeof(cmplx));
354 
355  *out = output;
356  return flen2;
357 }
358 
359 // Version for asymmetrical sidebands. You have to double the FFT size used for SSB.
360 int fftfilt::runAsym(const cmplx & in, cmplx **out, bool usb)
361 {
362  data[inptr++] = in;
363  if (inptr < flen2)
364  return 0;
365  inptr = 0;
366 
367  fft->ComplexFFT(data);
368 
369  data[0] *= filter[0]; // always keep DC
370 
371  if (usb)
372  {
373  for (int i = 1; i < flen2; i++)
374  {
375  data[i] *= filter[i]; // usb
376  data[flen2 + i] *= filterOpp[flen2 + i]; // lsb is the opposite
377  }
378  }
379  else
380  {
381  for (int i = 1; i < flen2; i++)
382  {
383  data[i] *= filterOpp[i]; // usb is the opposite
384  data[flen2 + i] *= filter[flen2 + i]; // lsb
385  }
386  }
387 
388  // in-place FFT: freqdata overwritten with filtered timedata
390 
391  // overlap and add
392  for (int i = 0; i < flen2; i++) {
393  output[i] = ovlbuf[i] + data[i];
394  ovlbuf[i] = data[i+flen2];
395  }
396 
397  memset (data, 0, flen * sizeof(cmplx));
398 
399  *out = output;
400  return flen2;
401 }
402 
403 /* Sliding FFT from Fldigi */
404 
408 } ;
409 
410 sfft::sfft(int len)
411 {
412  vrot_bins = new vrot_bins_pair[len];
413  delay = new cmplx[len];
414  fftlen = len;
415  first = 0;
416  last = len - 1;
417  ptr = 0;
418  double phi = 0.0, tau = 2.0 * M_PI/ len;
419  k2 = 1.0;
420  for (int i = 0; i < len; i++) {
421  vrot_bins[i].vrot = cmplx( K1 * cos (phi), K1 * sin (phi) );
422  phi += tau;
423  delay[i] = vrot_bins[i].bins = 0.0;
424  k2 *= K1;
425  }
426 }
427 
429 {
430  delete [] vrot_bins;
431  delete [] delay;
432 }
433 
434 // Sliding FFT, cmplx input, cmplx output
435 // FFT is computed for each value from first to last
436 // Values are not stable until more than "len" samples have been processed.
437 void sfft::run(const cmplx& input)
438 {
439  cmplx & de = delay[ptr];
440  const cmplx z( input.real() - k2 * de.real(), input.imag() - k2 * de.imag());
441  de = input;
442 
443  if (++ptr >= fftlen)
444  ptr = 0;
445 
446  for (vrot_bins_pair *itr = vrot_bins + first, *end = vrot_bins + last; itr != end ; ++itr)
447  itr->bins = (itr->bins + z) * itr->vrot;
448 }
449 
450 // Copies the frequencies to a pointer.
451 void sfft::fetch(float *result)
452 {
453  for (vrot_bins_pair *itr = vrot_bins, *end = vrot_bins + last; itr != end; ++itr, ++result)
454  *result = itr->bins.real() * itr->bins.real()
455  + itr->bins.imag() * itr->bins.imag();
456 }
457 
int pass
Definition: fftfilt.h:49
int flen
Definition: fftfilt.h:40
Fixed< IntType, IntBits > cos(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2271
cmplx * ovlbuf
Definition: fftfilt.h:46
cmplx frrc(float fb, float a, int i, int len)
Definition: fftfilt.h:68
fftfilt(float f1, float f2, int len)
Definition: fftfilt.cpp:78
void InverseComplexFFT(std::complex< FFT_TYPE > *buf)
Definition: gfft.h:3325
int runFilt(const cmplx &in, cmplx **out)
Definition: fftfilt.cpp:260
~fftfilt()
Definition: fftfilt.cpp:96
Fixed< IntType, IntBits > abs(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2313
#define M_PI
Definition: rdsdemod.cpp:27
sfft(int len)
Definition: fftfilt.cpp:410
float _blackman(int i, int len)
Definition: fftfilt.h:59
float fsinc(float fc, int i, int len)
Definition: fftfilt.h:52
std::complex< float > cmplx
Definition: fftfilt.h:21
void ComplexFFT(std::complex< FFT_TYPE > *buf)
Definition: gfft.h:3309
#define K1
Definition: fftfilt.h:97
int noFilt(const cmplx &in, cmplx **out)
Definition: fftfilt.cpp:248
void fetch(float *result)
Definition: fftfilt.cpp:451
void run(const cmplx &input)
Definition: fftfilt.cpp:437
void create_filter(float f1, float f2)
Definition: fftfilt.cpp:107
Fixed< IntType, IntBits > sin(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2265
cmplx * filterOpp
Definition: fftfilt.h:44
int32_t i
Definition: decimators.h:244
~sfft()
Definition: fftfilt.cpp:428
void init_filter()
Definition: fftfilt.cpp:51
cmplx * filter
Definition: fftfilt.h:43
void create_rrc_filter(float fb, float a)
root raised cosine. fb is half the band pass
Definition: fftfilt.cpp:222
int runSSB(const cmplx &in, cmplx **out, bool usb, bool getDC=true)
Definition: fftfilt.cpp:284
int flen2
Definition: fftfilt.h:41
g_fft< float > * fft
Definition: fftfilt.h:42
void create_asym_filter(float fopp, float fin)
two different filters for in band and opposite band
Definition: fftfilt.cpp:174
cmplx * output
Definition: fftfilt.h:47
cmplx * data
Definition: fftfilt.h:45
std::complex< float > cmplx
Definition: fftfilt.h:99
int inptr
Definition: fftfilt.h:48
int runDSB(const cmplx &in, cmplx **out, bool getDC=true)
Definition: fftfilt.cpp:327
void create_dsb_filter(float f2)
Definition: fftfilt.cpp:148
int runAsym(const cmplx &in, cmplx **out, bool usb)
Asymmetrical fitering can be used for vestigial sideband.
Definition: fftfilt.cpp:360
int window
Definition: fftfilt.h:50