fft filter

* Changed fftfilt implementation to use g_fft class vice
    Cfft class.
  * Removed multi-channel decoding when not visible or configured
    for decoding when not visible.

src/cw_rtty/rtty.cxx

@@ -190,21 +190,17 @@ void rtty::reset_filters()
 	int filter_length = 1024;
 
 	if (mark_filt) {
-		mark_filt->rtty_order(rtty_baud/samplerate, 0,
-			progdefaults.rtty_filter, 1.0);
+		mark_filt->rtty_filter(rtty_baud/samplerate);
 	} else {
 		mark_filt = new fftfilt(rtty_baud/samplerate, filter_length);
-		mark_filt->rtty_order(rtty_baud/samplerate, 0,
-			progdefaults.rtty_filter, 1.0);
+		mark_filt->rtty_filter(rtty_baud/samplerate);
      }
 
 	if (space_filt) {
-		space_filt->rtty_order(rtty_baud/samplerate, 0,
-			progdefaults.rtty_filter, 1.0);
+		space_filt->rtty_filter(rtty_baud/samplerate);
 	} else {
 		space_filt = new fftfilt(rtty_baud/samplerate, filter_length);
-		space_filt->rtty_order(rtty_baud/samplerate, 0,
-			progdefaults.rtty_filter, 1.0);
+		space_filt->rtty_filter(rtty_baud/samplerate);
      }
 }
 
@@ -584,10 +580,9 @@ int rtty::rx_process(const double *buf, int len)
 	int n_out = 0;
 	static int bitcount = 5 * nbits * symbollen;
 
-	if (rttyviewer && !bHistory && 
-		(!progdefaults.report_when_visible ||
-		 (progdefaults.report_when_visible && (dlgViewer->visible() || progStatus.show_channels))))
-		rttyviewer->rx_process(buf, len);
+	if ( !progdefaults.report_when_visible ||
+		 dlgViewer->visible() || progStatus.show_channels )
+		if (!bHistory && rttyviewer) rttyviewer->rx_process(buf, len);
 
 	if (progStatus.rtty_filter_changed) {
 		progStatus.rtty_filter_changed = false;

src/cw_rtty/view_rtty.cxx

@@ -105,21 +105,19 @@ view_rtty::~view_rtty()
 
 void view_rtty::reset_filters(int ch)
 {
-// filter_length = 512 / 1024 / 2048
-// use extended rcos(0) filter implementation
 	int filter_length = 1024;
 	if (channel[ch].mark_filt) {
-		channel[ch].mark_filt->rtty_order(rtty_baud/samplerate, 0);
+		channel[ch].mark_filt->rtty_filter(rtty_baud/samplerate);
 	} else {
 		channel[ch].mark_filt = new fftfilt(rtty_baud/samplerate, filter_length);
-		channel[ch].mark_filt->rtty_order(rtty_baud/samplerate, 0);
+		channel[ch].mark_filt->rtty_filter(rtty_baud/samplerate);
 	}
 
 	if (channel[ch].space_filt) {
-		channel[ch].space_filt->rtty_order(rtty_baud/samplerate, 0);
+		channel[ch].space_filt->rtty_filter(rtty_baud/samplerate);
 	} else {
 		channel[ch].space_filt = new fftfilt(rtty_baud/samplerate, filter_length);
-		channel[ch].space_filt->rtty_order(rtty_baud/samplerate, 0);
+		channel[ch].space_filt->rtty_filter(rtty_baud/samplerate);
 	}
 }
 

src/filters/fftfilt.cxx

@@ -27,471 +27,274 @@
 // along with fldigi.  If not, see <http://www.gnu.org/licenses/>.
 // ----------------------------------------------------------------------------
 
-#include <fstream>
-
 #include <config.h>
+
 #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 <cmath>
 #include "misc.h"
-
 #include "fftfilt.h"
 
+//------------------------------------------------------------------------------
+// initialize the filter
+// create forward and reverse FFTs
+//------------------------------------------------------------------------------
 
-fftfilt::fftfilt(double f1, double f2, int len)
+// probably only need a single instance of g_fft !!
+// use for both forward and reverse
+
+void fftfilt::init_filter()
 {
-	filterlen = len;
-	fft = new Cfft(filterlen);
-	ift = new Cfft(filterlen);
+	flen2 = flen >> 1;
+	fft			= new g_fft<double>(flen);
 
-	ovlbuf		= new cmplx[filterlen/2];
-	filter		= new cmplx[filterlen];
-	filtdata	= new cmplx[filterlen];
-	ht			= new cmplx[filterlen];
+	filter		= new cmplx[flen];
+	timedata	= new cmplx[flen];
+	freqdata	= new cmplx[flen];
+	output		= new cmplx[flen];
+	ovlbuf		= new cmplx[flen2];
+	ht			= new cmplx[flen];
 
-	for (int i = 0; i < filterlen; i++)
-		filter[i].real() = filter[i].imag() =
-		filtdata[i].real() = filtdata[i].imag() = 0.0;
-	for (int i = 0; i < filterlen/2; i++)
-		ovlbuf[i].real() = ovlbuf[i].imag() = 0.0;
+	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(double f1, double f2, int len)
+{
+	flen	= len;
+	init_filter();
 	create_filter(f1, f2);
 }
 
+//------------------------------------------------------------------------------
+// low pass filter
+//------------------------------------------------------------------------------
 fftfilt::fftfilt(double f, int len)
 {
-	filterlen = len;
-	fft = new Cfft(filterlen);
-	ift = new Cfft(filterlen);
-
-	ovlbuf		= new cmplx[filterlen/2];
-	filter		= new cmplx[filterlen];
-	filtdata	= new cmplx[filterlen];
-	ht			= new cmplx[filterlen];
-
-	for (int i = 0; i < filterlen; i++)
-		filter[i].real() = filter[i].imag() =
-		filtdata[i].real() = filtdata[i].imag() = 0.0;
-	for (int i = 0; i < filterlen/2; i++)
-		ovlbuf[i].real() = ovlbuf[i].imag() = 0.0;
-
-	inptr = 0;
-
+	flen	= len;
+	init_filter();
 	create_lpf(f);
 }
 
 fftfilt::~fftfilt()
 {
 	if (fft) delete fft;
-	if (ift) delete ift;
-	if (ovlbuf) delete [] ovlbuf;
+
 	if (filter) delete [] filter;
-	if (filtdata) delete [] filtdata;
+	if (timedata) delete [] timedata;
+	if (freqdata) delete [] freqdata;
+	if (output) delete [] output;
+	if (ovlbuf) delete [] ovlbuf;
 	if (ht) delete [] ht;
 }
 
+void fftfilt::create_filter(double f1, double 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);
+//	fft->transform(ht, filter);
+
+// normalize the output filter for unity gain
+	double 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;
+	}
+
+// perform the reverse fft to obtain h(t)
+// for testing
+// uncomment to obtain filter characteristics
+/*
+	cmplx *revht = new cmplx[flen];
+	memcpy(revht, filter, flen * sizeof(cmplx));
+
+	fft->InverseComplexFFT(revht);
+
+	std::fstream fspec;
+	fspec.open("fspec.csv", std::ios::out);
+	fspec << "i,imp.re,imp.im,filt.re,filt.im,filt.abs,revimp.re,revimp.im\n";
+	for (int i = 0; i < flen2; i++)
+		fspec
+			<< i << "," << ht[i].real() << "," << ht[i].imag() << ","
+			<< filter[i].real() << "," << filter[i].imag() << ","
+			<< abs(filter[i]) << ","
+			<< revht[i].real() << "," << revht[i].imag() << ","
+			<< std::endl;
+	fspec.close();
+	delete [] revht;
+*/
+	pass = 2;
+}
+
 /*
  * Filter with fast convolution (overlap-add algorithm).
  */
-int fftfilt::run(const cmplx& in, cmplx **out)
-{
-// collect filterlen/2 input samples
-	const int filterlen_div2 = filterlen / 2 ;
-	filtdata[inptr++] = in;
 
-	if (inptr < filterlen_div2)
+int fftfilt::run(const cmplx & in, cmplx **out)
+{
+// collect flen/2 input samples
+	timedata[inptr++] = in;
+
+	if (inptr < flen2)
 		return 0;
 	if (pass) --pass; // filter output is not stable until 2 passes
 
-// zero the rest of the input data
-	for (int i = filterlen_div2 ; i < filterlen; i++)
-		filtdata[i].real() = filtdata[i].imag() = 0.0;
-
 // FFT transpose to the frequency domain
-	fft->cdft(filtdata);
+	memcpy(freqdata, timedata, flen * sizeof(cmplx));
+	fft->ComplexFFT(freqdata);
 
 // multiply with the filter shape
-	for (int i = 0; i < filterlen; i++)
-		filtdata[i] *= filter[i];
+	for (int i = 0; i < flen; i++)
+		freqdata[i] *= filter[i];
 
-// IFFT transpose back to the time domain
-	ift->icdft(filtdata);
+// transform back to time domain
+	fft->InverseComplexFFT(freqdata);
 
 // overlap and add
-	for (int i = 0; i < filterlen_div2; i++) {
-		filtdata[i] += ovlbuf[i];
+// save the second half for overlapping next inverse FFT
+	for (int i = 0; i < flen2; i++) {
+		output[i] = ovlbuf[i] + freqdata[i];
+		ovlbuf[i] = freqdata[i+flen2];
 	}
-	*out = filtdata;
-
-// save the second half for overlapping
-	// Memcpy is allowed because cmplx are POD objects.
-	memcpy( ovlbuf, filtdata + filterlen_div2, sizeof( ovlbuf[0] ) * filterlen_div2 );
 
 // clear inbuf pointer
 	inptr = 0;
 
-// signal the caller there is filterlen/2 samples ready
+// signal the caller there is flen/2 samples ready
 	if (pass) return 0;
 
-	return filterlen_div2;
+	*out = output;
+	return flen2;
 }
 
-void fftfilt::create_filter(double f1, double f2)
-{
-	int len = filterlen / 2 + 1;
-	double t, h, x, it;
-	Cfft *tmpfft;
-	tmpfft = new Cfft(filterlen);
-
-// initialize the filter to zero
-	for (int i = 0; i < filterlen; i++)
-		filter[i].real()   = filter[i].imag()   = 0.0;
-
-// create the filter shape coefficients by fft
-// filter values initialized to the impulse response h(t)
-	for (int i = 0; i < len; i++) {
-		it = (double) i;
-		t = it - (len - 1) / 2.0;
-		h = it / (len - 1);
-
-		x = f2 * sinc(2 * f2 * t) - f1 * sinc(2 * f1 * t); // sinc(x)
-//		x *= hamming(t);
-//		x *= hanning(h);
-		x *= blackman(h);	// windowed by Blackman function
-		x *= filterlen;		// scaled for unity in passband
-		filter[i].real() = x;
-	}
-// perform the cmplx forward fft to obtain H(w)
-	tmpfft->cdft(filter);
-// start outputs after 2 full passes are complete
-	pass = 2;
-	delete tmpfft;
-}
-
-void fftfilt::create_lpf(double f)
-{
-	int len = filterlen / 2 + 1;
-	double t, h, x, it;
-	Cfft *tmpfft;
-	tmpfft = new Cfft(filterlen);
-
-// initialize the filter to zero
-	for (int i = 0; i < filterlen; i++)
-		filter[i].real()   = filter[i].imag()   = 0.0;
-
-// create the filter shape coefficients by fft
-// filter values initialized to the impulse response h(t)
-	for (int i = 0; i < len; i++) {
-		it = (double) i;
-		t = it - (len - 1) / 2.0;
-		h = it / (len - 1);
-
-		x = f * sinc(2 * f * t);
-		x *= blackman(h);	// windowed by Blackman function
-		x *= filterlen;		// scaled for unity in passband
-		filter[i].real() = x;
-	}
-// perform the cmplx forward fft to obtain H(w)
-	tmpfft->cdft(filter);
-// start outputs after 2 full passes are complete
-	pass = 2;
-	delete tmpfft;
-}
+//------------------------------------------------------------------------------
+// rtty filter
+//------------------------------------------------------------------------------
 
 //bool print_filter = true; // flag to inhibit printing multiple copies
 
-void fftfilt::create_rttyfilt(double f)
+void fftfilt::rtty_filter(double f)
 {
-	int len = filterlen / 2 + 1;
-	double t, h, it;
-	Cfft *tmpfft;
-	tmpfft = new Cfft(filterlen);
+// Raised cosine filter designed iaw Section 1.2.6 of
+// Telecommunications Measurements, Analysis, and Instrumentation
+// by Dr. Kamilo Feher / Engineers of Hewlett-Packard
+//
+// Frequency scaling factor determined hueristically by testing various values 
+// and measuring resulting decoder CER with input s/n = - 9 dB
+//
+//    K     CER
+//   1.0   .0244
+//   1.1   .0117
+//   1.2   .0081
+//   1.3   .0062
+//   1.4   .0054
+//   1.5   .0062
+//   1.6   .0076
 
-	// initialize the filter to zero
-	for (int i = 0; i < filterlen; i++)
-		filter[i].real()   = filter[i].imag()   = 0.0;
+	f *= 1.4;
 
-	// get an array to hold the sinc-respose
-	double* sinc_array = new double[ len ];
+	double dht;
+	for( int i = 0; i < flen2; ++i ) {
+		double x = (double)i/(double)(flen2);	
 
-	// create the impulse-response in it
-	for (int i = 0; i < len; ++i) {
-		it = (double)i;
-		t  = it - ( (double)len - 1.0) / 2.0;
-		h  = it / ( (double)len - 1.0);
+// raised cosine response (changed for -1.0...+1.0 times Nyquist-f
+// instead of books versions ranging from -1..+1 times samplerate)
 
-		// create the filter impulses with an additional zero at 1.5f
-		// remark: sinc(..) is scaled by 2, see misc.h
+		dht =
+			x <= 0 ? 1.0 :
+			x > 2.0 * f ? 0.0 :
+			cos((M_PI * x) / (f * 4.0));
 
-// Modified Lanzcos filter see http://en.wikipedia.org/wiki/Lanczos_resampling
-		sinc_array[i] = 
-			( sinc( 3.0 * f * t			 ) +
-			  sinc( 3.0 * f * t - 1.0	   ) * 0.8 +
-			  sinc( 3.0 * f * t + 1.0	   ) * 0.8 ) *
-			( sinc( 4.0 * f * t / 3.0	   ) +
-			  sinc( 4.0 * f * t / 3.0 - 1.0 ) * 0.8 +
-			  sinc( 4.0 * f * t / 3.0 + 1.0 ) * 0.8 );
+		dht *= dht; // cos^2
+
+// amplitude equalized nyquist-channel response
+		dht /= sinc(2.0 * i * f);
+
+		filter[i].real() = dht*cos((double)i* - 0.5*M_PI);
+		filter[i].imag() = dht*sin((double)i* - 0.5*M_PI);
+
+		filter[(flen-i)%flen].real() = dht*cos((double)i*+0.5*M_PI);
+		filter[(flen-i)%flen].imag() = dht*sin((double)i*+0.5*M_PI);
 	}
 
-	// normalize the impulse-responses
-	double sum = 0.0;
-	for (int i = 0; i < len; ++i) {
-		sum += sinc_array[i];
-		}
-	for (int i = 0; i < len; ++i) {
-		sinc_array[i] /= 8*sum;
-		}
-
-	// setup windowed-filter
-	for (int i = 0; i < len; ++i) {
-		it = (double)i;
-		t  = it - ( (double)len - 1.0) / 2.0;
-		h  = it / ( (double)len - 1.0);
-
-		filter[i].real() = ( sinc_array[i] ) * (double)filterlen * blackman(h);
-		sinc_array[i] = filter[i].real();
-		}
-
+// perform the reverse fft to obtain h(t)
+// for testing
+// uncomment to obtain filter characteristics
 /*
-// create an identical filter impulse response for testing
-// ht_B should be identical to ht_A within limits of math processing
-// Hw is the frequency response of filter created using ht_A impulse
-// response
-	Cfft test_fft(filterlen);
-	cmplx ht_A[filterlen]; // original impulse response
-	cmplx ht_B[filterlen]; // computed impulse response
-	cmplx Hw[filterlen];   // computed H(w)
+	cmplx *revht = new cmplx[flen];
+	memcpy(revht, filter, flen * sizeof(cmplx));
 
-// ht_A retains the original normalized impulse response
-// ht_B used for forward / reverse FFT
-	for (int i = 0; i < len; ++i)
-		ht_B[i] = ht_A[i] = filter[i];
+	fft->InverseComplexFFT(revht);
 
-// perform the cmplx forward fft to obtain H(w)
-	test_fft.cdft(ht_B);
-	for (int i = 0; i < len; ++i)
-		Hw[i] = ht_B[i];
-
-// perform the cmplx reverse fft to obtain h(t) again
-	test_fft.icdft(ht_B);
-
-// ht_B should be equal to ht_A
-	std::fstream file1("filter_debug.csv", std::ios::out );
-	for (int i = 0; i < len; ++i)
-		file1 << ht_A[i].real() << "," << ht_A[i].imag() << "," << 
-				ht_B[i].real() << "," << ht_B[i].imag() << "," <<
-				Hw[i].real() << "," << Hw[i].imag() << "," << Hw[i].mag() << "\n";
-	file1.close();
+	std::fstream fspec;
+	fspec.open("rtty_filter.csv", std::ios::out);
+	fspec << "i,filt.re,filt.im,filt.abs,,revimp.re,revimp.im\n";
+	for (int i = 0; i < flen; i++)
+		fspec
+			<< i << ","
+			<< filter[i].real() << "," << filter[i].imag() << "," << abs(filter[i]) 
+			<< ",," << revht[i].real() << "," << revht[i].imag()
+			<< std::endl;
+	fspec.close();
+	delete [] revht;
 */
-
-// perform the cmplx forward fft to obtain H(w)
-	tmpfft->cdft(filter);
-/*
-	if (print_filter) {
-		std::fstream file2("filter_response.csv", std::ios::out );
-		file2 << "Modified Lanzcos 1.5 stop bit filter\n\n";
-		file2 << "h(t), |H(w)|, dB\n\n";
-		double dc = 20*log10(filter[0].mag());
-		for (int i = 0; i < len; i++)
-			file2 << sinc_array[i] << "," << filter[i].mag() << ","
-					20*log10(filter[i].mag()) - dc << "\n";
-		file2.close();
-		print_filter = false;
-	}
-*/
-
 // start outputs after 2 full passes are complete
 	pass = 2;
-	delete tmpfft;
-	delete [] sinc_array;
-
 }
 
-double xrcos(double t, double T, int order, double alpha = 1)
-{
-	if (order == 1) return rcos(t, T, alpha);
-	order--;
-	return xrcos(2*t - T/2, T, order, alpha) + xrcos(2*t + T/2, T, order, alpha);
-}
-
-double stefan(double t, double T)
-{
-// Stefan implementation
-	double a=.7;
-	double h = rcos( t		, T/4.0, a );
-	h += rcos( t - T/4.0, T/4.0, a );
-	h += rcos( t + T/4.0, T/4.0, a );
-	return h;
-}
-
-double matched(double t, double T)
-{
-	if (t > -T/2 && t < T/2) return 1;
-	return 0;
-}
-
-double sinc_filter(double t, double T)
-{
-	return sinc(t / T);
-}
-
-void fftfilt::rtty_order(double f, int N, double twarp, double alpha)
-{
-	int len = filterlen / 2 + 1;
-	double ft;
-	Cfft tmpfft(filterlen);
-
-	// create the impulse-response
-	for (int i = 0; i < filterlen; ++i) {
-		if (i > len) {
-			ht[i].real() = ht[i].imag() = 0.0;
-			continue;
-		}
-		ft = f * (1.0* i - len / 2.0);
-		switch(N) {
-		default:
-		case 0:
-			ft *= twarp; // compromise filter CPFSK vs SHAPED_AFSK
-			ht[i] = xrcos( ft, 1.0, 1, alpha);
-			break;
-		case 1:
-			ft *= 1.1;
-			ht[i] = xrcos( ft, 1.0, 2, alpha );
-			break;
-		case 2:
-//			ft *= 1.0;
-			ht[i] = xrcos( ft, 1.0, 3, alpha );
-			break;
-		case 3:
-			ft *= 1.5;
-			ht[i] = rcos( ft, 1.0 );
-			break;
-		case 4:
-			ft *= (1.0 + M_PI/2.0);
-			ht[i]  = rcos( ft - 0.5, 1.0 );
-			ht[i] += rcos( ft + 0.5, 1.0 );
-			break;
-		case 5:
-			ft *= (3.0 + M_PI/2.0);
-			ht[i]  = rcos( ft - 1.5, 1.0 );
-			ht[i] += rcos( ft - 0.5, 1.0 );
-			ht[i] += rcos( ft + 0.5, 1.0 );
-			ht[i] += rcos( ft + 1.5, 1.0 );
-			break;
-		case 6:
-			ft *= M_PI / 2.0;
-			ht[i] = sinc_filter(ft, 1.0 );
-			break;
-		case 7:
-			ft *= (2.0 + M_PI/2.0);
-			ht[i]  = rcos( ft - 1.0, 1.0 );
-			ht[i] += rcos( ft, 1.0 );
-			ht[i] += rcos( ft + 1.0, 1.0 );
-			break;
-		case 8:
-//			ft *= 1.0+0.57079/10.0E10; // simulating inf
-			ht[i] = matched(ft, 1.0);
-			break;
-		}
-	}
-
-// normalize the impulse-response
-	double sum = 0.0;
-	for (int i = 0; i <= len; ++i) {
-		sum += ht[i].real();
-	}
-	for (int i = 0; i < filterlen; ++i) {
-		ht[i].real() *= filterlen/sum;
-		filter[i] = ht[i];
-	}
-
-/*
-// create an identical filter impulse response for testing
-// ht_B should be identical to ht_A within limits of math processing
-// Hw is the frequency response of filter created using ht_A impulse
-// response
-	Cfft test_fft(filterlen);
-	cmplx ht_A[filterlen]; // original impulse response
-	cmplx ht_B[filterlen]; // computed impulse response
-	cmplx Hw[filterlen];   // computed H(w)
-
-// ht_A retains the original normalized impulse response
-// ht_B used for forward / reverse FFT
-	for (int i = 0; i < filterlen; i++)
-		ht_B[i] = ht_A[i] = filter[i];
-
-// perform the cmplx forward fft to obtain H(w)
-	test_fft.cdft(ht_B);
-	for (int i = 0; i < filterlen; i++)
-		Hw[i] = ht_B[i];
-
-// perform the cmplx reverse fft to obtain h(t) again
-	test_fft.icdft(ht_B);
-
-// ht_B should be equal to ht_A
-	std::fstream file1("filter_debug.csv", std::ios::out );
-	for (int i = 0; i < filterlen; i++)//len; ++i)
-		file1 << ht_A[i].real() << "," << ht_B[i].real() << "," 
-			  << ht_A[i].real() - ht_B[i].real() << ","
-			  << Hw[i].mag() << "\n";
-	file1.close();
-*/
-
-// perform the cmplx forward fft to obtain H(w)
-//	tmpfft->cdft(filter);
-	tmpfft.cdft(filter);
-
-// start outputs after 2 full passes are complete
-	pass = 2;
-//	delete tmpfft;
-
-// Stefan's latest
-
-	f*=1.275; // This factor is ominous to me. I can't explain it. It shouldn't
-			  // be there. But if I leave it out ht(f) differs inbetween the 
-			  // raised cosine from above and this one. And if left out the error
-			  // rate increases... So, this is an unsolved mystery for now.
-
-	for( int i = 0; i < filterlen/2; ++i ) {
-		double a = 1.0; 
-		double x = (double)i/(double)(filterlen/2);		
-
-	// raised cosine response (changed for -1.0...+1.0 times Nyquist-f
-	// instead of books versions ranging from -1..+1 times samplerate)
-
-		double ht =
-			fabs(x) <= (1.0 - a)/(1.0/f) ? 1.0:
-			fabs(x) >  (1.0 + a)/(1.0/f) ? 0.0:
-			cos(M_PI/(f*4.0*a)*(fabs(x)-(1.0-a)/(1.0/f)));
-		ht *= ht; // cos^2
-
-		// equalized nyquist-channel response
-		double eq = 1.0/sinc((double)i*f*2);  
-
-		// compensate for "awkward" FFT-implementation. For every other imple-
-		// mentation of a FFT this would have been just...
-
-		filter[i].real() = eq*ht*sin((double)i* - 0.5*M_PI);
-		filter[i].imag() = eq*ht*cos((double)i* - 0.5*M_PI);
-
-		filter[(filterlen-i)%filterlen].real() = eq*ht*sin((double)i*+0.5*M_PI);
-		filter[(filterlen-i)%filterlen].imag() = eq*ht*cos((double)i*+0.5*M_PI);
-
-		// ... this (caused most headache):
-		//filter[i].real() = eq*ht*0.7071;
-		//filter[i].imag() = eq*ht*0.7071;
-		//filter[(filterlen-i)%filterlen].real() = eq*ht*0.7071;
-		//filter[(filterlen-i)%filterlen].imag() = eq*ht*0.7071;
-
-	}
-	std::fstream file1("filter_debug.csv", std::ios::out );
-	for (int i = 0; i < filterlen/2; i++)
-		file1 << filter[i].real() << "," << filter[i].imag() << ","
-			  << abs(filter[i]) << "\n";
-	file1.close();
-
-}

src/include/complex.h

@@ -40,131 +40,4 @@ inline cmplx cmac (const cmplx *a, const cmplx *b, int ptr, int len) {
 	return z;
 }
 
-// conjprod = conj(x) * y;
-// previously %
-//inline cmplx conjprod(cmplx &x, cmplx& y) {
-//	cmplx z;
-//	z = (x.real() * y.real() + x.imag() * y.imag(),
-//		 x.real() * y.imag() - x.imag() * y.real());
-//	return z;
-//}
-
-/*
-class complex {
-public:
-	double re;
-	double im;
-	complex(double r = 0.0, double i = 0.0)
-	    : re(r), im(i) { }
-
-	double real() { return re; };
-	void real(double R) {re = R;};
-	double imag() { return im; };
-	void imag(double I) {im = I;};
-
-// Z = X * Y
-	complex& operator*=(const complex& y) {
-		double temp = re * y.re - im * y.im;
-		im = re * y.im + im * y.re;
-		re = temp;
-		return *this;
-	}
-	complex operator*(const complex& y) const {
-		return complex(re * y.re - im * y.im,  re * y.im + im * y.re);
-	}
-
-// Z = X * y
-	complex& operator*=(double y) {
-		re *= y;
-		im *= y;
-		return *this;
-	}
-	complex operator*(double y) const {
-		return complex(re * y, im * y);
-	}
-
-// Z = X + Y
-	complex& operator+=(const complex& y) {
-		re += y.re;
-		im += y.im;
-		return *this;
-        }
-	complex operator+(const complex& y) const {
-		return complex(re + y.re,  im + y.im);
-	}
-
-// Z = X - Y
-	complex& operator-=(const complex& y) {
-		re -= y.re;
-		im -= y.im;
-		return *this;
-	}
-	complex operator-(const complex& y) const {
-		return complex(re - y.re,  im - y.im);
-	}
-
-// Z = X / Y
-	complex& operator/=(const complex& y) {
-		double temp, denom = y.re*y.re + y.im*y.im;
-		if (denom == 0.0) denom = 1e-10;
-		temp = (re * y.re + im * y.im) / denom;
-		im = (im * y.re - re * y.im) / denom;
-		re = temp;
-		return *this;
-	}
-	complex operator/(const complex& y) const {
-		double denom = y.re*y.re + y.im*y.im;
-		if (denom == 0.0) denom = 1e-10;
-		return complex((re * y.re + im * y.im) / denom,  (im * y.re - re * y.im) / denom);
-	}
-	
-// Z = (complex conjugate of X) * Y
-// Z1 = x1 + jy1, or Z1 = |Z1|exp(jP1)
-// Z2 = x2 + jy2, or Z2 = |Z2|exp(jP2)
-// Z = (x1 - jy1) * (x2 + jy2)
-// or Z = |Z1|*|Z2| exp (j (P2 - P1))
-	complex& operator%=(const complex& y) {
-		double temp = re * y.re + im * y.im;
-		im = re * y.im - im * y.re;
-		re = temp;
-		return *this;
-	}
-	complex operator%(const complex& y) const {
-		complex z;
-		z.re = re * y.re + im * y.im;
-		z.im = re * y.im - im * y.re;
-		return z;
-	}
-
-// n = |Z| * |Z| 	
-	double norm() const {
-		return (re * re + im * im);
-	}
-
-// n = |Z|
-	double mag() const {
-		return sqrt(norm());
-	}
-
-// Z = x + jy
-// Z = |Z|exp(jP)
-// arg returns P
-	double arg() const {
-		return atan2(im, re);
-	}
-
-};
-
-inline 	complex cmac (const complex *a, const complex *b, int ptr, int len) {
-		complex z;
-		ptr %= len;
-		for (int i = 0; i < len; i++) {
-			z += a[i] * b[ptr];
-			ptr = (ptr + 1) % len;
-		}
-		return z;
-	}
-	
-*/
-
 #endif

src/include/fftfilt.h

@@ -6,7 +6,7 @@
 #define	_FFTFILT_H
 
 #include "complex.h"
-#include "fft.h"
+#include "gfft.h"
 
 //----------------------------------------------------------------------
 
@@ -14,24 +14,45 @@ class fftfilt {
 enum {NONE, BLACKMAN, HAMMING, HANNING};
 
 protected:
-	int filterlen;
-    Cfft *fft;
-    Cfft *ift;
-	cmplx *filter;
-	cmplx *filtdata;
-	cmplx *ovlbuf;
+	int flen;
+	int flen2;
+	g_fft<double> *fft;
+	g_fft<double> *ift;
 	cmplx *ht;
+	cmplx *filter;
+	cmplx *timedata;
+	cmplx *freqdata;
+	cmplx *ovlbuf;
+	cmplx *output;
 	int inptr;
 	int pass;
 	int window;
+
+	inline double fsinc(double 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 double _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(double f1, double f2, int len);
 	fftfilt(double f, int len);
 	~fftfilt();
+// f1 < f2 ==> bandpass
+// f1 > f2 ==> band reject
 	void create_filter(double f1, double f2);
-	void create_lpf(double f);
-	void create_rttyfilt(double f);
-	void rtty_order(double, int, double twarp = 1.275, double alpha = 1.0);
+	void create_lpf(double f) {
+		create_filter(0, f);
+	}
+	void create_hpf(double f) {
+		create_filter(f, 0);
+	}
+	void rtty_filter(double);
 
 	int run(const cmplx& in, cmplx **out);
 };

src/include/misc.h

@@ -86,8 +86,9 @@ inline double hanning(double x)
 inline double rcos( double t, double T, double alpha=1.0 )
 {
     if( t == 0 ) return 1.0;
-    if( fabs(t) == ( T/(2.0*alpha) ) ) 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-(2.0*alpha*t/T)*(2.0*alpha*t/T));
+    double 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));
 }
 
 // Rectangular - no pre filtering of data array

src/psk/psk.cxx

@@ -1182,10 +1182,9 @@ int psk::rx_process(const double *buf, int len)
 	bool can_rx_symbol = false;
 
 	if (numcarriers == 1) {
-		if (pskviewer && !bHistory && 
-			(!progdefaults.report_when_visible ||
-			 (progdefaults.report_when_visible && (dlgViewer->visible() || progStatus.show_channels))))
-			pskviewer->rx_process(buf, len);
+		if (!progdefaults.report_when_visible ||
+			 dlgViewer->visible() || progStatus.show_channels )
+			if (pskviewer && !bHistory) pskviewer->rx_process(buf, len);
 		if (evalpsk)
 			evalpsk->sigdensity();
 	}