GreenFFT.

CMakeLists.txt

@@ -49,6 +49,7 @@ set(sdrbase_SOURCES
 	sdrbase/dsp/dspcommands.cpp
 	sdrbase/dsp/dspengine.cpp
 	sdrbase/dsp/fftengine.cpp
+	sdrbase/dsp/fftfilt.cxx
 	sdrbase/dsp/fftwindow.cpp
 	sdrbase/dsp/interpolator.cpp
 	sdrbase/dsp/inthalfbandfilter.cpp
@@ -107,17 +108,21 @@ set(sdrbase_HEADERS
 
 	include-gpl/dsp/channelizer.h
 	include/dsp/channelmarker.h
+	include-gpl/dsp/complex.h
 	include-gpl/dsp/dspcommands.h
 	include-gpl/dsp/dspengine.h
 	include/dsp/dsptypes.h
 	include-gpl/dsp/fftengine.h
+	include-gpl/dsp/fftfilt.h
 	include-gpl/dsp/fftwengine.h
 	include-gpl/dsp/fftwindow.h
+	include-gpl/dsp/gfft.h
 	include-gpl/dsp/interpolator.h
 	include-gpl/dsp/inthalfbandfilter.h
 	include/dsp/kissfft.h
 	include-gpl/dsp/kissengine.h
 	include-gpl/dsp/lowpass.h
+	include-gpl/dsp/misc.h
 	include-gpl/dsp/movingaverage.h
 	include-gpl/dsp/nco.h
 	sdrbase/dsp/pidcontroller.h

include-gpl/dsp/complex.h

@@ -0,0 +1,43 @@
+// ----------------------------------------------------------------------------
+// complex.h  --  Complex arithmetic
+//
+// Copyright (C) 2006-2008
+//		Dave Freese, W1HKJ
+// Copyright (C) 2008
+//		Stelios Bounanos, M0GLD
+//
+// This file is part of fldigi.
+//
+// Fldigi is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// Fldigi is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with fldigi.  If not, see <http://www.gnu.org/licenses/>.
+// ----------------------------------------------------------------------------
+
+#ifndef _COMPLEX_H
+#define _COMPLEX_H
+
+#include <cmath>
+#include <complex>
+
+typedef std::complex<float> cmplx;
+
+inline cmplx cmac (const cmplx *a, const cmplx *b, int ptr, int len) {
+	cmplx z;
+	ptr %= len;
+	for (int i = 0; i < len; i++) {
+		z += a[i] * b[ptr];
+		ptr = (ptr + 1) % len;
+		}
+	return z;
+}
+
+#endif

include-gpl/dsp/fftfilt.h

@@ -0,0 +1,54 @@
+/*
+ *    fftfilt.h  --  Fast convolution FIR filter
+*/
+
+#ifndef	_FFTFILT_H
+#define	_FFTFILT_H
+
+#include "complex.h"
+#include "gfft.h"
+
+//----------------------------------------------------------------------
+
+class fftfilt {
+enum {NONE, BLACKMAN, HAMMING, HANNING};
+
+protected:
+	int flen;
+	int flen2;
+	g_fft<float> *fft;
+	g_fft<float> *ift;
+	cmplx *ht;
+	cmplx *filter;
+	cmplx *timedata;
+	cmplx *freqdata;
+	cmplx *ovlbuf;
+	cmplx *output;
+	int inptr;
+	int pass;
+	int window;
+
+	inline float fsinc(float fc, int i, int len) {
+		return (i == len/2) ? 2.0 * fc: 
+				sin(2 * M_PI * fc * (i - len/2)) / (M_PI * (i - len/2));
+	}
+	inline float _blackman(int i, int len) {
+		return (0.42 - 
+				 0.50 * cos(2.0 * M_PI * i / len) + 
+				 0.08 * cos(4.0 * M_PI * i / len));
+	}
+	void init_filter();
+
+public:
+	fftfilt(float f1, float f2, int len);
+	fftfilt(float f, int len);
+	~fftfilt();
+// f1 < f2 ==> bandpass
+// f1 > f2 ==> band reject
+	void create_filter(float f1, float f2);
+	void rtty_filter(float);
+
+	int run(const cmplx& in, cmplx **out);
+};
+
+#endif

include-gpl/dsp/misc.h

@@ -0,0 +1,75 @@
+// ----------------------------------------------------------------------------
+// misc.h  --  Miscellaneous helper functions
+//
+// Copyright (C) 2006-2008
+//		Dave Freese, W1HKJ
+//
+// This file is part of fldigi.  These filters were adapted from code contained
+// in the gmfsk source code distribution.
+//
+// Fldigi is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// Fldigi is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with fldigi.  If not, see <http://www.gnu.org/licenses/>.
+// ----------------------------------------------------------------------------
+
+#ifndef _MISC_H
+#define _MISC_H
+
+#include <cmath>
+
+inline float sinc(float x)
+{
+	return (fabs(x) < 1e-10) ? 1.0 : (sin(M_PI * x) / (M_PI * x));
+}
+
+inline float cosc(float x)
+{
+	return (fabs(x) < 1e-10) ? 0.0 : ((1.0 - cos(M_PI * x)) / (M_PI * x));
+}
+
+inline float clamp(float x, float min, float max)
+{
+	return (x < min) ? min : ((x > max) ? max : x);
+}
+
+/// This is always called with an int weight
+inline float decayavg(float average, float input, int weight)
+{
+	if (weight <= 1) return input;
+	return ( ( input - average ) / (float)weight ) + average ;
+}
+
+// following are defined inline to provide best performance
+inline float blackman(float x)
+{
+	return (0.42 - 0.50 * cos(2 * M_PI * x) + 0.08 * cos(4 * M_PI * x));
+}
+
+inline float hamming(float x)
+{
+	return 0.54 - 0.46 * cos(2 * M_PI * x);
+}
+
+inline float hanning(float x)
+{
+	return 0.5 - 0.5 * cos(2 * M_PI * x);
+}
+
+inline float rcos( float t, float T, float alpha=1.0 )
+{
+    if( t == 0 ) return 1.0;
+    float taT = T / (2.0 * alpha);
+    if( fabs(t) == taT ) return ((alpha/2.0) * sin(M_PI/(2.0*alpha)));
+    return (sin(M_PI*t/T)/(M_PI*t/T))*cos(alpha*M_PI*t/T)/(1.0-(t/taT)*(t/taT));
+}
+
+#endif

sdrbase/dsp/fftfilt.cxx

@@ -0,0 +1,172 @@
+// ----------------------------------------------------------------------------
+//	fftfilt.cxx  --  Fast convolution Overlap-Add filter
+//
+// Filter implemented using overlap-add FFT convolution method
+// h(t) characterized by Windowed-Sinc impulse response
+//
+// Reference:
+//	 "The Scientist and Engineer's Guide to Digital Signal Processing"
+//	 by Dr. Steven W. Smith, http://www.dspguide.com
+//	 Chapters 16, 18 and 21
+//
+// Copyright (C) 2006-2008 Dave Freese, W1HKJ
+//
+// This file is part of fldigi.
+//
+// Fldigi is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// Fldigi is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with fldigi.  If not, see <http://www.gnu.org/licenses/>.
+// ----------------------------------------------------------------------------
+
+#include <memory.h>
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+#include <cmath>
+#include <typeinfo>
+
+#include <stdio.h>
+#include <sys/types.h>
+#include <unistd.h>
+#include <memory.h>
+
+#include <dsp/misc.h>
+#include <dsp/fftfilt.h>
+
+//------------------------------------------------------------------------------
+// initialize the filter
+// create forward and reverse FFTs
+//------------------------------------------------------------------------------
+
+// Only need a single instance of g_fft, used for both forward and reverse
+void fftfilt::init_filter()
+{
+	flen2	= flen >> 1;
+	fft	= new g_fft<float>(flen);
+
+	filter		= new cmplx[flen];
+	timedata	= new cmplx[flen];
+	freqdata	= new cmplx[flen];
+	output		= new cmplx[flen];
+	ovlbuf		= new cmplx[flen2];
+	ht		= new cmplx[flen];
+
+	memset(filter, 0, flen * sizeof(cmplx));
+	memset(timedata, 0, flen * sizeof(cmplx));
+	memset(freqdata, 0, flen * sizeof(cmplx));
+	memset(output, 0, flen * sizeof(cmplx));
+	memset(ovlbuf, 0, flen2 * sizeof(cmplx));
+	memset(ht, 0, flen * sizeof(cmplx));
+
+	inptr = 0;
+}
+
+//------------------------------------------------------------------------------
+// fft filter
+// f1 < f2 ==> band pass filter
+// f1 > f2 ==> band reject filter
+// f1 == 0 ==> low pass filter
+// f2 == 0 ==> high pass filter
+//------------------------------------------------------------------------------
+fftfilt::fftfilt(float f1, float f2, int len)
+{
+	flen	= len;
+	init_filter();
+	create_filter(f1, f2);
+}
+
+fftfilt::~fftfilt()
+{
+	if (fft) delete fft;
+
+	if (filter) delete [] filter;
+	if (timedata) delete [] timedata;
+	if (freqdata) delete [] freqdata;
+	if (output) delete [] output;
+	if (ovlbuf) delete [] ovlbuf;
+	if (ht) delete [] ht;
+}
+
+void fftfilt::create_filter(float f1, float f2)
+{
+// initialize the filter to zero
+	memset(ht, 0, flen * sizeof(cmplx));
+
+// create the filter shape coefficients by fft
+// filter values initialized to the ht response h(t)
+	bool b_lowpass, b_highpass;//, window;
+	b_lowpass = (f2 != 0);
+	b_highpass = (f1 != 0);
+
+	for (int i = 0; i < flen2; i++) {
+		ht[i] = 0;
+//combine lowpass / highpass
+// lowpass @ f2
+		if (b_lowpass) ht[i] += fsinc(f2, i, flen2);
+// highighpass @ f1
+		if (b_highpass) ht[i] -= fsinc(f1, i, flen2);
+	}
+// highpass is delta[flen2/2] - h(t)
+	if (b_highpass && f2 < f1) ht[flen2 / 2] += 1;
+
+	for (int i = 0; i < flen2; i++)
+		ht[i] *= _blackman(i, flen2);
+
+// this may change since green fft is in place fft
+	memcpy(filter, ht, flen * sizeof(cmplx));
+
+// ht is flen complex points with imaginary all zero
+// first half describes h(t), second half all zeros
+// perform the cmplx forward fft to obtain H(w)
+// filter is flen/2 complex values
+
+	fft->ComplexFFT(filter);
+
+// normalize the output filter for unity gain
+	float scale = 0, mag;
+	for (int i = 0; i < flen2; i++) {
+		mag = abs(filter[i]);
+		if (mag > scale) scale = mag;
+	}
+	if (scale != 0) {
+		for (int i = 0; i < flen; i++)
+			filter[i] /= scale;
+	}
+}
+
+// Filter with fast convolution (overlap-add algorithm).
+int fftfilt::run(const cmplx & in, cmplx **out)
+{
+	timedata[inptr++] = in;
+
+	if (inptr < flen2)
+		return 0;
+	inptr = 0;
+
+	memcpy(freqdata, timedata, flen * sizeof(cmplx));
+	fft->ComplexFFT(freqdata);
+
+	for (int i = 0; i < flen; i++)
+		freqdata[i] *= filter[i];
+
+	fft->InverseComplexFFT(freqdata);
+
+	// overlap and add
+	for (int i = 0; i < flen2; i++) {
+		output[i] = ovlbuf[i] + freqdata[i];
+		ovlbuf[i] = freqdata[i+flen2];
+	}
+
+	*out = output;
+	return flen2;
+}
+