From 26455261ee700845243d1f3bdc7212c24ce8d40f Mon Sep 17 00:00:00 2001 From: f4exb Date: Mon, 9 Feb 2026 22:18:03 +0100 Subject: [PATCH] New FFT RRC filter --- sdrbase/CMakeLists.txt | 2 + sdrbase/dsp/fftfilterrrc.cpp | 188 +++++++++++++++++++++++++++++++++++ sdrbase/dsp/fftfilterrrc.h | 148 +++++++++++++++++++++++++++ sdrbench/CMakeLists.txt | 1 + sdrbench/mainbench.cpp | 2 + sdrbench/mainbench.h | 1 + sdrbench/parserbench.cpp | 4 +- sdrbench/parserbench.h | 3 +- sdrbench/test_fftrrc.cpp | 73 ++++++++++++++ sdrbench/test_fftrrc.py | 41 ++++++++ 10 files changed, 461 insertions(+), 2 deletions(-) create mode 100644 sdrbase/dsp/fftfilterrrc.cpp create mode 100644 sdrbase/dsp/fftfilterrrc.h create mode 100644 sdrbench/test_fftrrc.cpp create mode 100755 sdrbench/test_fftrrc.py diff --git a/sdrbase/CMakeLists.txt b/sdrbase/CMakeLists.txt index 0aaafa68f..d89f30201 100644 --- a/sdrbase/CMakeLists.txt +++ b/sdrbase/CMakeLists.txt @@ -145,6 +145,7 @@ set(sdrbase_SOURCES dsp/fftengine.cpp dsp/fftfactory.cpp dsp/fftfilt.cpp + dsp/fftfilterrrc.cpp dsp/fftnr.cpp dsp/fftwindow.cpp dsp/filterrc.cpp @@ -373,6 +374,7 @@ set(sdrbase_HEADERS dsp/fftengine.h dsp/fftfactory.h dsp/fftfilt.h + dsp/fftfilterrrc.h dsp/fftnr.h dsp/fftwengine.h dsp/fftwindow.h diff --git a/sdrbase/dsp/fftfilterrrc.cpp b/sdrbase/dsp/fftfilterrrc.cpp new file mode 100644 index 000000000..23229eda0 --- /dev/null +++ b/sdrbase/dsp/fftfilterrrc.cpp @@ -0,0 +1,188 @@ +/////////////////////////////////////////////////////////////////////////////////// +// Copyright (C) 2012 maintech GmbH, Otto-Hahn-Str. 15, 97204 Hoechberg, Germany // +// written by Christian Daniel // +// Copyright (C) 2015-2019, 2023 Edouard Griffiths, F4EXB // +// Copyright (C) 2026 SDRangel contributors // +// // +// This program 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 as version 3 of the License, or // +// (at your option) any later version. // +// // +// This program 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 V3 for more details. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +/////////////////////////////////////////////////////////////////////////////////// + +#include "dsp/fftfilterrrc.h" + +#include +#include +#include + +FFTFilterRRC::FFTFilterRRC(int len) : + m_fftLen(len), + m_fftLen2(len / 2), + m_fft(nullptr), + m_filter(nullptr), + m_fftBuffer(nullptr), + m_overlapBuffer(nullptr), + m_outputBuffer(nullptr), + m_inputIndex(0), + m_symbolRate(0.0f), + m_rolloff(0.0f) +{ + initFilter(); +} + +FFTFilterRRC::~FFTFilterRRC() +{ + delete m_fft; + delete[] m_filter; + delete[] m_fftBuffer; + delete[] m_overlapBuffer; + delete[] m_outputBuffer; +} + +void FFTFilterRRC::initFilter() +{ + m_fft = new g_fft(m_fftLen); + + m_filter = new Complex[m_fftLen]; + m_fftBuffer = new Complex[m_fftLen]; + m_overlapBuffer = new Complex[m_fftLen2]; + m_outputBuffer = new Complex[m_fftLen2]; + + // Initialize all buffers to zero + std::fill(m_filter, m_filter + m_fftLen, Complex(0.0f, 0.0f)); + std::fill(m_fftBuffer, m_fftBuffer + m_fftLen, Complex(0.0f, 0.0f)); + std::fill(m_overlapBuffer, m_overlapBuffer + m_fftLen2, Complex(0.0f, 0.0f)); + std::fill(m_outputBuffer, m_outputBuffer + m_fftLen2, Complex(0.0f, 0.0f)); + + m_inputIndex = 0; +} + +void FFTFilterRRC::create(float symbolRate, float rolloff) +{ + m_symbolRate = symbolRate; + m_rolloff = rolloff; + + // Clamp rolloff to valid range [0, 1] + if (m_rolloff < 0.0f) { + m_rolloff = 0.0f; + } else if (m_rolloff > 1.0f) { + m_rolloff = 1.0f; + } + + // Initialize filter to zero + std::fill(m_filter, m_filter + m_fftLen, Complex(0.0f, 0.0f)); + + // Compute frequency-domain RRC response for each FFT bin + for (int i = 0; i < m_fftLen; i++) { + m_filter[i] = computeRRCResponse(m_symbolRate, m_rolloff, i, m_fftLen); + } + + // Normalize for unity gain in passband + float maxMag = 0.0f; + for (int i = 0; i < m_fftLen; i++) { + float mag = std::abs(m_filter[i]); + if (mag > maxMag) { + maxMag = mag; + } + } + + if (maxMag > 0.0f) { + for (int i = 0; i < m_fftLen; i++) { + m_filter[i] /= maxMag; + } + } +} + +FFTFilterRRC::Complex FFTFilterRRC::computeRRCResponse( + float symbolRate, + float rolloff, + int binIndex, + int fftLen) const +{ + // Normalize frequency to [-0.5, 0.5] + // FFT bins: 0 to fftLen/2-1 are positive frequencies (0 to 0.5) + // fftLen/2 to fftLen-1 are negative frequencies (-0.5 to 0) + float freq = static_cast(binIndex) / static_cast(fftLen); + + // Map to [-0.5, 0.5] range + if (freq > 0.5f) { + freq -= 1.0f; + } + + // Absolute frequency + float absFreq = std::abs(freq); + + // Symbol time (inverse of symbol rate) + float T = 1.0f / symbolRate; + + // Compute frequency boundaries + float f1 = (1.0f - rolloff) / (2.0f * T); // Start of transition band + float f2 = (1.0f + rolloff) / (2.0f * T); // End of transition band + + Complex response; + + if (absFreq <= f1) { + // Passband: constant response + response = Complex(std::sqrt(T), 0.0f); + } else if (absFreq > f1 && absFreq <= f2) { + // Transition band: raised cosine roll-off + // H(f) = sqrt(T) * 0.5 * [1 + cos(pi*T/beta * (|f| - (1-beta)/(2T)))] + float arg = (M_PI * T / rolloff) * (absFreq - f1); + float amplitude = std::sqrt(T) * 0.5f * (1.0f + std::cos(arg)); + response = Complex(amplitude, 0.0f); + } else { + // Stopband: zero response + response = Complex(0.0f, 0.0f); + } + + return response; +} + +int FFTFilterRRC::process(const Complex& input, Complex** output) +{ + // Store input sample in first half of FFT buffer + m_fftBuffer[m_inputIndex] = input; + m_inputIndex++; + + // Process when first half is full + if (m_inputIndex < m_fftLen2) { + return 0; // Not enough samples yet + } + + // Second half of FFT buffer is already zero-padded from initialization + // or previous processing (overlap-add method) + + // Perform forward FFT + m_fft->ComplexFFT(m_fftBuffer); + + // Apply filter in frequency domain (element-wise multiplication) + for (int i = 0; i < m_fftLen; i++) { + m_fftBuffer[i] *= m_filter[i]; + } + + // Perform inverse FFT + m_fft->InverseComplexFFT(m_fftBuffer); + + // Overlap-add: first half gets added to overlap buffer, second half becomes new overlap + for (int i = 0; i < m_fftLen2; i++) { + m_outputBuffer[i] = m_fftBuffer[i] + m_overlapBuffer[i]; + m_overlapBuffer[i] = m_fftBuffer[i + m_fftLen2]; + } + + // Reset input index and clear first half of FFT buffer for next block + m_inputIndex = 0; + std::fill(m_fftBuffer, m_fftBuffer + m_fftLen2, Complex(0.0f, 0.0f)); + + // Return pointer to output buffer + *output = m_outputBuffer; + return m_fftLen2; +} diff --git a/sdrbase/dsp/fftfilterrrc.h b/sdrbase/dsp/fftfilterrrc.h new file mode 100644 index 000000000..78d48320e --- /dev/null +++ b/sdrbase/dsp/fftfilterrrc.h @@ -0,0 +1,148 @@ +/////////////////////////////////////////////////////////////////////////////////// +// Copyright (C) 2012 maintech GmbH, Otto-Hahn-Str. 15, 97204 Hoechberg, Germany // +// written by Christian Daniel // +// Copyright (C) 2015-2019, 2023 Edouard Griffiths, F4EXB // +// Copyright (C) 2026 SDRangel contributors // +// // +// This program 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 as version 3 of the License, or // +// (at your option) any later version. // +// // +// This program 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 V3 for more details. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +/////////////////////////////////////////////////////////////////////////////////// + +#ifndef INCLUDE_DSP_FFTFILTERRRC_H +#define INCLUDE_DSP_FFTFILTERRRC_H + +#include +#include + +#include "gfft.h" +#include "export.h" + +/** + * @brief FFT-based root raised cosine filter for complex signals + * + * Implements a frequency-domain root raised cosine (RRC) filter using overlap-add + * FFT convolution. The RRC filter is commonly used for pulse shaping in digital + * communications to minimize intersymbol interference (ISI) while satisfying + * the Nyquist criterion. + * + * The filter operates in the frequency domain for computational efficiency, + * applying the RRC frequency response directly to the FFT bins of the input signal. + */ +class SDRBASE_API FFTFilterRRC +{ +public: + using Complex = std::complex; + + /** + * @brief Construct FFT-based root raised cosine filter + * + * @param len FFT length (should be power of 2 for efficiency) + */ + explicit FFTFilterRRC(int len); + + /** + * @brief Destructor + */ + ~FFTFilterRRC(); + + // Disable copy semantics (manages raw pointers) + FFTFilterRRC(const FFTFilterRRC&) = delete; + FFTFilterRRC& operator=(const FFTFilterRRC&) = delete; + + /** + * @brief Create root raised cosine filter in frequency domain + * + * @param symbolRate Symbol rate in Hz (normalized to sample rate) + * @param rolloff Roll-off factor (beta) - typically 0.2 to 0.5 + * 0 = rectangular, 1 = full cosine roll-off + * + * The filter bandwidth extends from -symbolRate*(1+rolloff)/2 to + * +symbolRate*(1+rolloff)/2 in normalized frequency (0 to 0.5 = Nyquist). + */ + void create(float symbolRate, float rolloff); + + /** + * @brief Process one complex input sample through the filter + * + * Uses overlap-add FFT convolution. Accumulates input samples until + * a half-buffer is filled, then performs FFT, frequency-domain multiplication + * with the RRC filter, IFFT, and overlap-add to produce output samples. + * + * @param input Complex input sample + * @param output Pointer to output buffer (will point to internal buffer) + * @return Number of output samples available (0 or fftLen/2) + */ + int process(const Complex& input, Complex** output); + + /** + * @brief Get FFT length + * @return FFT length in samples + */ + int getFFTLength() const { return m_fftLen; } + + /** + * @brief Get current symbol rate + * @return Symbol rate (normalized) + */ + float getSymbolRate() const { return m_symbolRate; } + + /** + * @brief Get current rolloff factor + * @return Rolloff factor (beta) + */ + float getRolloff() const { return m_rolloff; } + + /** + * @brief Get pointer to frequency-domain filter coefficients + * + * The filter coefficients are in the frequency domain, corresponding to the FFT bins. + * + * @return Pointer to array of complex filter coefficients (length = FFT length) + */ + const Complex* getFilter() const { return m_filter; } + +private: + /** + * @brief Initialize internal buffers and FFT engine + */ + void initFilter(); + + /** + * @brief Compute RRC frequency response for a given frequency bin + * + * Implements the RRC transfer function in frequency domain: + * - For |f| <= (1-beta)/(2T): H(f) = T + * - For (1-beta)/(2T) < |f| <= (1+beta)/(2T): H(f) = T/2 * [1 + cos(pi*T/beta * (|f| - (1-beta)/(2T)))] + * - For |f| > (1+beta)/(2T): H(f) = 0 + * + * @param symbolRate Symbol rate (normalized to sample rate) + * @param rolloff Roll-off factor (beta) + * @param binIndex FFT bin index + * @param fftLen FFT length + * @return Complex frequency response at this bin + */ + Complex computeRRCResponse(float symbolRate, float rolloff, int binIndex, int fftLen) const; + + int m_fftLen; //!< FFT length + int m_fftLen2; //!< Half FFT length + g_fft* m_fft; //!< FFT engine + Complex* m_filter; //!< Frequency-domain filter coefficients + Complex* m_fftBuffer; //!< Time-domain input buffer for FFT + Complex* m_overlapBuffer; //!< Overlap buffer for overlap-add + Complex* m_outputBuffer; //!< Output buffer + int m_inputIndex; //!< Current input buffer index + float m_symbolRate; //!< Current symbol rate (normalized) + float m_rolloff; //!< Current rolloff factor +}; + +#endif // INCLUDE_DSP_FFTFILTERRRC_H diff --git a/sdrbench/CMakeLists.txt b/sdrbench/CMakeLists.txt index 9eee7766c..4cfcd5c46 100644 --- a/sdrbench/CMakeLists.txt +++ b/sdrbench/CMakeLists.txt @@ -11,6 +11,7 @@ set(sdrbench_SOURCES test_ft8.cpp test_callsign.cpp test_ft8protocols.cpp + test_fftrrc.cpp ) set(sdrbench_HEADERS diff --git a/sdrbench/mainbench.cpp b/sdrbench/mainbench.cpp index 8544ec083..06e6cd200 100644 --- a/sdrbench/mainbench.cpp +++ b/sdrbench/mainbench.cpp @@ -73,6 +73,8 @@ void MainBench::run() testCallsign(m_parser.getArgsStr()); } else if (m_parser.getTestType() == ParserBench::TestFT8Protocols) { testFT8Protocols(m_parser.getArgsStr()); + } else if (m_parser.getTestType() == ParserBench::TestFFTRRCFilter) { + testFFTRRCFilter(); } else { qDebug() << "MainBench::run: unknown test type: " << m_parser.getTestType(); } diff --git a/sdrbench/mainbench.h b/sdrbench/mainbench.h index c2417d9d8..41efa9a3b 100644 --- a/sdrbench/mainbench.h +++ b/sdrbench/mainbench.h @@ -57,6 +57,7 @@ private: void testDecimateFI(); void testDecimateFF(); void testGolay2312(); + void testFFTRRCFilter(); void testFT8(const QString& wavFile, const QString& argsStr); //!< use with sdrbench/samples/ft8/230105_091630.wav in -f option void testFT8Protocols(const QString& argsStr); void testCallsign(const QString& argsStr); diff --git a/sdrbench/parserbench.cpp b/sdrbench/parserbench.cpp index f38d1ebe2..1490f0a84 100644 --- a/sdrbench/parserbench.cpp +++ b/sdrbench/parserbench.cpp @@ -24,7 +24,7 @@ ParserBench::ParserBench() : m_testOption(QStringList() << "t" << "test", - "Test type: decimateii, decimatefi, decimateff, decimateif, decimateinfii, decimatesupii, ambe, golay2312, ft8, ft8protocols, callsign" + "Test type: decimateii, decimatefi, decimateff, decimateif, decimateinfii, decimatesupii, ambe, golay2312, ft8, ft8protocols, callsign, fftrrcfilter.", "test", "decimateii"), m_nbSamplesOption(QStringList() << "n" << "nb-samples", @@ -151,6 +151,8 @@ ParserBench::TestType ParserBench::getTestType() const return TestCallsign; } else if (m_testStr == "ft8protocols") { return TestFT8Protocols; + } else if (m_testStr == "fftrrcfilter") { + return TestFFTRRCFilter; } else { return TestDecimatorsII; } diff --git a/sdrbench/parserbench.h b/sdrbench/parserbench.h index 05476ce61..874c70dd5 100644 --- a/sdrbench/parserbench.h +++ b/sdrbench/parserbench.h @@ -40,7 +40,8 @@ public: TestGolay2312, TestFT8, TestCallsign, - TestFT8Protocols + TestFT8Protocols, + TestFFTRRCFilter } TestType; ParserBench(); diff --git a/sdrbench/test_fftrrc.cpp b/sdrbench/test_fftrrc.cpp new file mode 100644 index 000000000..ecf29c3c8 --- /dev/null +++ b/sdrbench/test_fftrrc.cpp @@ -0,0 +1,73 @@ +/////////////////////////////////////////////////////////////////////////////////// +// Copyright (C) 2026 Edouard Griffiths, F4EXB // +// // +// This program 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 as version 3 of the License, or // +// (at your option) any later version. // +// // +// This program 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 V3 for more details. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +/////////////////////////////////////////////////////////////////////////////////// + +#include "mainbench.h" +#include "dsp/fftfilterrrc.h" +#include + +void MainBench::testFFTRRCFilter() +{ + qDebug() << "MainBench::testFFTRRCFilter"; + const int RRC_FFT_SIZE = 512; + std::complex* rrcFilterOut = nullptr; + FFTFilterRRC filter(RRC_FFT_SIZE); + filter.create(0.05f, 0.35f); // 2400 baud, 0.35 rolloff + + qDebug() << "MainBench::testFFTRRCFilter: filter created"; + FILE *fd_filter = fopen("test_fftrrc_filter.txt", "w"); + + for (int i = 0; i < RRC_FFT_SIZE; i++) { + fprintf(fd_filter, "%f\n", std::abs(filter.getFilter()[i])); + } + qDebug() << "MainBench::testFFTRRCFilter: filter coefficients written to test_fftrrc_filter.txt"; + fclose(fd_filter); + + qDebug() << "MainBench::testFFTRRCFilter: running filter"; + FILE *fd = fopen("test_fftrrc.txt", "w"); + int outLen = 0; + + for (int i = 0; i < 1024; i++) + { + Real phi = i * (48000.0 / 1200.0) * (2*3.141); + Real x = sin(phi); + filter.process(Complex(x, 0.0f), &rrcFilterOut); + outLen++; + + if (outLen % (RRC_FFT_SIZE / 2) == 0) + { + // printf("\ni=%d n_out: %d\n", i, n_out); + + for (int j = 0; j < outLen; j++) + { + Complex rrc = rrcFilterOut[j]; + fprintf(fd, "%f\n", rrc.real()); + } + + outLen = 0; + } + } + + // output any remaining samples + for (int j = 0; j < outLen; j++) + { + Complex rrc = rrcFilterOut[j]; + fprintf(fd, "%f\n", rrc.real()); + } + + qDebug() << "MainBench::testFFTRRCFilter: output samples written to test_fftrrc.txt"; + fclose(fd); +} diff --git a/sdrbench/test_fftrrc.py b/sdrbench/test_fftrrc.py new file mode 100755 index 000000000..f25790bc2 --- /dev/null +++ b/sdrbench/test_fftrrc.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python3 +import matplotlib.pyplot as plt +import numpy as np + +with open('../build/test_fftrrc_filter.txt', 'r') as f: + filter_data = f.read() + +filter_out = [float(x) for x in filter_data.splitlines()] +xf = np.array(filter_out) + +with open('../build/test_fftrrc.txt', 'r') as f: + data = f.read() + +out = [float(x) for x in data.splitlines()] +# marks = [127, 254, 383, 511, 639, 767, 895, 1023] +# marks = [127, 254, 383, 511] +x = np.array(out) + +fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 6)) + +ax2.plot(xf) +ax2.set_title('RRC Filter Coefficients') +ax2.set_xlabel('Index') +ax2.set_ylabel('Amplitude') +ax2.grid(True) + +ax1.set_title('FFTRRC Filter Output (Real Part)') +ax1.plot(x) +# ax1.scatter(marks, x[marks], +# color='red', marker='.', s=100, zorder=5, edgecolors='black') +# plt.annotate('126', xy=(126, x[126]), xytext=(6, 0), +# arrowprops=dict(facecolor='black', shrink=0.05), +# fontsize=12, ha='center') +ax1.set_xlabel('Index') +ax1.set_ylabel('Real part') +ax1.grid(True) + + +plt.subplots_adjust(hspace=0.4) +plt.savefig('plot.png', dpi=150, bbox_inches='tight') +plt.close() # Prevents memory leaks