From 75653763057bb9bf018dfde1a3cda8681c8b8e34 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Sat, 19 Aug 2017 11:05:01 +0200 Subject: Add WIP CFR implementation to OfdmGenerator --- src/OfdmGenerator.cpp | 134 +++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 127 insertions(+), 7 deletions(-) (limited to 'src/OfdmGenerator.cpp') diff --git a/src/OfdmGenerator.cpp b/src/OfdmGenerator.cpp index 20ba2b7..3d3dfd7 100644 --- a/src/OfdmGenerator.cpp +++ b/src/OfdmGenerator.cpp @@ -2,7 +2,7 @@ Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011 Her Majesty the Queen in Right of Canada (Communications Research Center Canada) - Copyright (C) 2016 + Copyright (C) 2017 Matthias P. Braendli, matthias.braendli@mpb.li http://opendigitalradio.org @@ -34,19 +34,23 @@ #include #include #include -typedef std::complex complexf; +#include OfdmGenerator::OfdmGenerator(size_t nbSymbols, size_t nbCarriers, size_t spacing, bool inverse) : - ModCodec(), - myFftPlan(NULL), - myFftIn(NULL), myFftOut(NULL), + ModCodec(), RemoteControllable("ofdm"), + myFftPlan(nullptr), + myFftIn(nullptr), myFftOut(nullptr), myNbSymbols(nbSymbols), myNbCarriers(nbCarriers), - mySpacing(spacing) + mySpacing(spacing), + myCfr(false), + myCfrClip(1.0f), + myCfrErrorClip(1.0f), + myCfrFft(nullptr) { PDEBUG("OfdmGenerator::OfdmGenerator(%zu, %zu, %zu, %s) @ %p\n", nbSymbols, nbCarriers, spacing, inverse ? "true" : "false", this); @@ -56,6 +60,11 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols, "OfdmGenerator::OfdmGenerator nbCarriers > spacing!"); } + /* register the parameters that can be remote controlled */ + RC_ADD_PARAMETER(cfr, "Enable crest factor reduction"); + RC_ADD_PARAMETER(cfrclip, "CFR: Clip to amplitude"); + RC_ADD_PARAMETER(cfrerrorclip, "CFR: Limit error"); + if (inverse) { myPosDst = (nbCarriers & 1 ? 0 : 1); myPosSrc = 0; @@ -91,6 +100,13 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols, myFftIn, myFftOut, FFTW_BACKWARD, FFTW_MEASURE); + myCfrPostClip.resize(N); + myCfrPostFft.resize(N); + myCfrFft = fftwf_plan_dft_1d(N, + reinterpret_cast(myCfrPostClip.data()), + reinterpret_cast(myCfrPostFft.data()), + FFTW_FORWARD, FFTW_MEASURE); + if (sizeof(complexf) != sizeof(FFT_TYPE)) { printf("sizeof(complexf) %zu\n", sizeof(complexf)); printf("sizeof(FFT_TYPE) %zu\n", sizeof(FFT_TYPE)); @@ -115,6 +131,10 @@ OfdmGenerator::~OfdmGenerator() if (myFftPlan) { fftwf_destroy_plan(myFftPlan); } + + if (myCfrFft) { + fftwf_destroy_plan(myCfrFft); + } } int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) @@ -157,7 +177,11 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) memcpy(&myFftIn[myNegDst], &in[myNegSrc], myNegSize * sizeof(FFT_TYPE)); - fftwf_execute(myFftPlan); + fftwf_execute(myFftPlan); // IFFT from myFftIn to myFftOut + + if (myCfr) { + cfr_one_iteration(dataOut, i); + } memcpy(out, myFftOut, mySpacing * sizeof(FFT_TYPE)); @@ -167,3 +191,99 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) return sizeOut; } +void OfdmGenerator::cfr_one_iteration(Buffer *symbols, size_t symbol_ix) +{ + complexf* out = reinterpret_cast(symbols->getData()); + + out += (symbol_ix * mySpacing); + + // use std::norm instead of std::abs to avoid calculating the + // square roots + const float clip_squared = myCfrClip * myCfrClip; + + // Clip + for (size_t i = 0; i < mySpacing; i++) { + const float mag_squared = std::norm(out[i]); + if (mag_squared > clip_squared) { + // normalise absolute value to myCfrClip: + // x_clipped = x * clip / |x| + // = x * sqrt(clip_squared) / sqrt(mag_squared) + // = x * sqrt(clip_squared / mag_squared) + out[i] *= std::sqrt(clip_squared / mag_squared); + } + } + + // Take FFT of our clipped signal + std::copy(out, out + mySpacing, myCfrPostClip.begin()); + fftwf_execute(myCfrFft); // FFT from myCfrPostClip to myCfrPostFft + + // Calculate the error in frequency domain by subtracting our reference + // and clip it to myCfrErrorClip. By adding this clipped error signal + // to our FFT output, we compensate the introduced error to some + // extent. + const float err_clip_squared = myCfrErrorClip * myCfrErrorClip; + + complexf *reference = reinterpret_cast(myFftIn); + for (size_t i = 0; i < mySpacing; i++) { + complexf error = myCfrPostFft[i] - reference[i]; + + const float mag_squared = std::norm(error); + if (mag_squared > err_clip_squared) { + error *= std::sqrt(err_clip_squared / mag_squared); + } + + // Update the reference in-place to avoid another copy for the + // subsequence IFFT + reference[i] += error; + } + + // Run our error-compensated symbol through the IFFT again + fftwf_execute(myFftPlan); // IFFT from myFftIn to myFftOut + +} + + +void OfdmGenerator::set_parameter(const std::string& parameter, + const std::string& value) +{ + using namespace std; + stringstream ss(value); + ss.exceptions ( stringstream::failbit | stringstream::badbit ); + + if (parameter == "cfr") { + ss >> myCfr; + } + else if (parameter == "clip") { + ss >> myCfrClip; + } + else if (parameter == "errorclip") { + ss >> myCfrErrorClip; + } + else { + stringstream ss; + ss << "Parameter '" << parameter + << "' is not exported by controllable " << get_rc_name(); + throw ParameterError(ss.str()); + } +} + +const std::string OfdmGenerator::get_parameter(const std::string& parameter) const +{ + using namespace std; + stringstream ss; + if (parameter == "cfr") { + ss << myCfr; + } + else if (parameter == "clip") { + ss << std::fixed << myCfrClip; + } + else if (parameter == "errorclip") { + ss << std::fixed << myCfrErrorClip; + } + else { + ss << "Parameter '" << parameter << + "' is not exported by controllable " << get_rc_name(); + throw ParameterError(ss.str()); + } + return ss.str(); +} -- cgit v1.2.3 From 571d57c591b23d212934f0bb969c33d15c729982 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Mon, 21 Aug 2017 18:14:55 +0200 Subject: Fix some CFR things, still WIP --- doc/example.ini | 11 ++++++++ src/ConfigParser.cpp | 7 +++++ src/ConfigParser.h | 6 +++++ src/DabModulator.cpp | 9 ++++++- src/OfdmGenerator.cpp | 73 ++++++++++++++++++++++++++++++++------------------- src/OfdmGenerator.h | 12 ++++++--- 6 files changed, 87 insertions(+), 31 deletions(-) (limited to 'src/OfdmGenerator.cpp') diff --git a/doc/example.ini b/doc/example.ini index 425dfa4..512e196 100644 --- a/doc/example.ini +++ b/doc/example.ini @@ -127,6 +127,17 @@ dac_clk_rate=0 ; and ;dac_clk_rate=128000000 +; Settings for crest factor reduction +[cfr] +enable=0 + +; At what amplitude the signal should be clipped +clip=1.0 + +; How much to clip the error signal used to compensate the effect +; of clipping +error_clip=1.0 + [firfilter] ; The FIR Filter can be used to create a better spectral quality. enabled=1 diff --git a/src/ConfigParser.cpp b/src/ConfigParser.cpp index 9ac1280..1cc94c0 100644 --- a/src/ConfigParser.cpp +++ b/src/ConfigParser.cpp @@ -177,6 +177,13 @@ static void parse_configfile( pt.get("poly.num_threads", 0); } + // Crest factor reduction + if (pt.get("cfr.enabled", 0) == 1) { + mod_settings.enableCfr = true; + mod_settings.cfrClip = pt.get("cfr.clip"); + mod_settings.cfrErrorClip = pt.get("cfr.error_clip"); + } + // Output options std::string output_selected; try { diff --git a/src/ConfigParser.h b/src/ConfigParser.h index 89f0fb7..a8d7837 100644 --- a/src/ConfigParser.h +++ b/src/ConfigParser.h @@ -77,6 +77,12 @@ struct mod_settings_t { std::string polyCoefFilename = ""; unsigned polyNumThreads = 0; + // Settings for crest factor reduction + bool enableCfr = false; + float cfrClip = 1.0f; + float cfrErrorClip = 1.0f; + + #if defined(HAVE_OUTPUT_UHD) OutputUHDConfig outputuhd_conf; #endif diff --git a/src/DabModulator.cpp b/src/DabModulator.cpp index cc2642a..0914469 100644 --- a/src/DabModulator.cpp +++ b/src/DabModulator.cpp @@ -183,7 +183,14 @@ int DabModulator::process(Buffer* dataOut) } auto cifOfdm = make_shared( - (1 + myNbSymbols), myNbCarriers, mySpacing); + (1 + myNbSymbols), + myNbCarriers, + mySpacing, + m_settings.enableCfr, + m_settings.cfrClip, + m_settings.cfrErrorClip); + + rcs.enrol(cifOfdm.get()); auto cifGain = make_shared( mySpacing, diff --git a/src/OfdmGenerator.cpp b/src/OfdmGenerator.cpp index 3d3dfd7..41c1c85 100644 --- a/src/OfdmGenerator.cpp +++ b/src/OfdmGenerator.cpp @@ -26,6 +26,8 @@ #include "OfdmGenerator.h" #include "PcDebug.h" + +#include #include "fftw3.h" #define FFT_TYPE fftwf_complex @@ -33,23 +35,25 @@ #include #include #include -#include #include OfdmGenerator::OfdmGenerator(size_t nbSymbols, - size_t nbCarriers, - size_t spacing, - bool inverse) : + size_t nbCarriers, + size_t spacing, + bool enableCfr, + float cfrClip, + float cfrErrorClip, + bool inverse) : ModCodec(), RemoteControllable("ofdm"), myFftPlan(nullptr), myFftIn(nullptr), myFftOut(nullptr), myNbSymbols(nbSymbols), myNbCarriers(nbCarriers), mySpacing(spacing), - myCfr(false), - myCfrClip(1.0f), - myCfrErrorClip(1.0f), + myCfr(enableCfr), + myCfrClip(cfrClip), + myCfrErrorClip(cfrErrorClip), myCfrFft(nullptr) { PDEBUG("OfdmGenerator::OfdmGenerator(%zu, %zu, %zu, %s) @ %p\n", @@ -62,8 +66,8 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols, /* register the parameters that can be remote controlled */ RC_ADD_PARAMETER(cfr, "Enable crest factor reduction"); - RC_ADD_PARAMETER(cfrclip, "CFR: Clip to amplitude"); - RC_ADD_PARAMETER(cfrerrorclip, "CFR: Limit error"); + RC_ADD_PARAMETER(clip, "CFR: Clip to amplitude"); + RC_ADD_PARAMETER(errorclip, "CFR: Limit error"); if (inverse) { myPosDst = (nbCarriers & 1 ? 0 : 1); @@ -100,11 +104,10 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols, myFftIn, myFftOut, FFTW_BACKWARD, FFTW_MEASURE); - myCfrPostClip.resize(N); - myCfrPostFft.resize(N); + myCfrPostClip = (FFT_TYPE*)fftwf_malloc(sizeof(FFT_TYPE) * N); + myCfrPostFft = (FFT_TYPE*)fftwf_malloc(sizeof(FFT_TYPE) * N); myCfrFft = fftwf_plan_dft_1d(N, - reinterpret_cast(myCfrPostClip.data()), - reinterpret_cast(myCfrPostFft.data()), + myCfrPostClip, myCfrPostFft, FFTW_FORWARD, FFTW_MEASURE); if (sizeof(complexf) != sizeof(FFT_TYPE)) { @@ -167,6 +170,9 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) "OfdmGenerator::process output size not valid!"); } + myNumClip = 0; + myNumErrorClip = 0; + for (size_t i = 0; i < myNbSymbols; ++i) { myFftIn[0][0] = 0; myFftIn[0][1] = 0; @@ -177,10 +183,17 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) memcpy(&myFftIn[myNegDst], &in[myNegSrc], myNegSize * sizeof(FFT_TYPE)); + std::vector reference; + if (myCfr) { + reference.resize(mySpacing); + memcpy(reference.data(), myFftIn, mySpacing * sizeof(FFT_TYPE)); + } + fftwf_execute(myFftPlan); // IFFT from myFftIn to myFftOut if (myCfr) { - cfr_one_iteration(dataOut, i); + complexf *symbol = reinterpret_cast(myFftOut); + cfr_one_iteration(symbol, reference.data()); } memcpy(out, myFftOut, mySpacing * sizeof(FFT_TYPE)); @@ -188,33 +201,36 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) in += myNbCarriers; out += mySpacing; } + + if (myCfr) { + etiLog.level(debug) << "CFR: " << myNumClip << " clipped, " << + myNumErrorClip << " err clipped"; + } + return sizeOut; } -void OfdmGenerator::cfr_one_iteration(Buffer *symbols, size_t symbol_ix) +void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *reference) { - complexf* out = reinterpret_cast(symbols->getData()); - - out += (symbol_ix * mySpacing); - // use std::norm instead of std::abs to avoid calculating the // square roots const float clip_squared = myCfrClip * myCfrClip; // Clip for (size_t i = 0; i < mySpacing; i++) { - const float mag_squared = std::norm(out[i]); + const float mag_squared = std::norm(symbol[i]); if (mag_squared > clip_squared) { // normalise absolute value to myCfrClip: // x_clipped = x * clip / |x| // = x * sqrt(clip_squared) / sqrt(mag_squared) // = x * sqrt(clip_squared / mag_squared) - out[i] *= std::sqrt(clip_squared / mag_squared); + symbol[i] *= std::sqrt(clip_squared / mag_squared); + myNumClip++; } } // Take FFT of our clipped signal - std::copy(out, out + mySpacing, myCfrPostClip.begin()); + memcpy(myCfrPostClip, symbol, mySpacing * sizeof(FFT_TYPE)); fftwf_execute(myCfrFft); // FFT from myCfrPostClip to myCfrPostFft // Calculate the error in frequency domain by subtracting our reference @@ -223,23 +239,26 @@ void OfdmGenerator::cfr_one_iteration(Buffer *symbols, size_t symbol_ix) // extent. const float err_clip_squared = myCfrErrorClip * myCfrErrorClip; - complexf *reference = reinterpret_cast(myFftIn); for (size_t i = 0; i < mySpacing; i++) { - complexf error = myCfrPostFft[i] - reference[i]; + const complexf constellation_point = + reinterpret_cast(myCfrPostFft)[i]; + + complexf error = reference[i] - constellation_point; const float mag_squared = std::norm(error); if (mag_squared > err_clip_squared) { error *= std::sqrt(err_clip_squared / mag_squared); + myNumErrorClip++; } - // Update the reference in-place to avoid another copy for the + // Update the input to the FFT directl to avoid another copy for the // subsequence IFFT - reference[i] += error; + complexf *fft_in = reinterpret_cast(myFftIn); + fft_in[i] = constellation_point + error; } // Run our error-compensated symbol through the IFFT again fftwf_execute(myFftPlan); // IFFT from myFftIn to myFftOut - } diff --git a/src/OfdmGenerator.h b/src/OfdmGenerator.h index 9a6fa01..a41d583 100644 --- a/src/OfdmGenerator.h +++ b/src/OfdmGenerator.h @@ -46,6 +46,9 @@ class OfdmGenerator : public ModCodec, public RemoteControllable OfdmGenerator(size_t nbSymbols, size_t nbCarriers, size_t spacing, + bool enableCfr, + float cfrClip, + float cfrErrorClip, bool inverse = true); virtual ~OfdmGenerator(); OfdmGenerator(const OfdmGenerator&) = delete; @@ -65,7 +68,7 @@ class OfdmGenerator : public ModCodec, public RemoteControllable const std::string& parameter) const override; protected: - void cfr_one_iteration(Buffer *symbols, size_t symbol_ix); + void cfr_one_iteration(complexf *symbol, const complexf *reference); fftwf_plan myFftPlan; fftwf_complex *myFftIn, *myFftOut; @@ -85,8 +88,11 @@ class OfdmGenerator : public ModCodec, public RemoteControllable float myCfrClip; float myCfrErrorClip; fftwf_plan myCfrFft; - std::vector myCfrPostClip; - std::vector myCfrPostFft; + fftwf_complex *myCfrPostClip; + fftwf_complex *myCfrPostFft; + + size_t myNumClip; + size_t myNumErrorClip; }; -- cgit v1.2.3 From 63ae20781af39e9cb6673e910e99e28f86c6adb0 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Mon, 21 Aug 2017 18:56:25 +0200 Subject: Normalise the CFR FFT --- src/OfdmGenerator.cpp | 12 ++++++++++-- src/OfdmGenerator.h | 6 +++--- 2 files changed, 13 insertions(+), 5 deletions(-) (limited to 'src/OfdmGenerator.cpp') diff --git a/src/OfdmGenerator.cpp b/src/OfdmGenerator.cpp index 41c1c85..8d1a560 100644 --- a/src/OfdmGenerator.cpp +++ b/src/OfdmGenerator.cpp @@ -173,6 +173,10 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) myNumClip = 0; myNumErrorClip = 0; + // It is not guaranteed that fftw keeps the FFT input vector intact. + // That's why we copy it to the reference. + std::vector reference; + for (size_t i = 0; i < myNbSymbols; ++i) { myFftIn[0][0] = 0; myFftIn[0][1] = 0; @@ -183,7 +187,6 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) memcpy(&myFftIn[myNegDst], &in[myNegSrc], myNegSize * sizeof(FFT_TYPE)); - std::vector reference; if (myCfr) { reference.resize(mySpacing); memcpy(reference.data(), myFftIn, mySpacing * sizeof(FFT_TYPE)); @@ -240,8 +243,13 @@ void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *referenc const float err_clip_squared = myCfrErrorClip * myCfrErrorClip; for (size_t i = 0; i < mySpacing; i++) { + // FFTW computes an unnormalised trasform, i.e. a FFT-IFFT pair + // or vice-versa give back the original vector scaled by a factor + // FFT-size. Because we're comparing our constellation point + // (calculated with IFFT-clip-FFT) against reference (input to + // the IFFT), we need to divide by our FFT size. const complexf constellation_point = - reinterpret_cast(myCfrPostFft)[i]; + reinterpret_cast(myCfrPostFft)[i] / (float)mySpacing; complexf error = reference[i] - constellation_point; diff --git a/src/OfdmGenerator.h b/src/OfdmGenerator.h index a41d583..ad39f5b 100644 --- a/src/OfdmGenerator.h +++ b/src/OfdmGenerator.h @@ -72,9 +72,9 @@ class OfdmGenerator : public ModCodec, public RemoteControllable fftwf_plan myFftPlan; fftwf_complex *myFftIn, *myFftOut; - size_t myNbSymbols; - size_t myNbCarriers; - size_t mySpacing; + const size_t myNbSymbols; + const size_t myNbCarriers; + const size_t mySpacing; unsigned myPosSrc; unsigned myPosDst; unsigned myPosSize; -- cgit v1.2.3 From 16c4bcb085457514438d35bbbe11ef979e36bb85 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Wed, 23 Aug 2017 08:35:26 +0200 Subject: Show CFR error min and max --- src/OfdmGenerator.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) (limited to 'src/OfdmGenerator.cpp') diff --git a/src/OfdmGenerator.cpp b/src/OfdmGenerator.cpp index 8d1a560..8f5d8a0 100644 --- a/src/OfdmGenerator.cpp +++ b/src/OfdmGenerator.cpp @@ -242,6 +242,8 @@ void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *referenc // extent. const float err_clip_squared = myCfrErrorClip * myCfrErrorClip; + std::vector error_norm(mySpacing); + for (size_t i = 0; i < mySpacing; i++) { // FFTW computes an unnormalised trasform, i.e. a FFT-IFFT pair // or vice-versa give back the original vector scaled by a factor @@ -254,17 +256,23 @@ void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *referenc complexf error = reference[i] - constellation_point; const float mag_squared = std::norm(error); + error_norm[i] = mag_squared; + if (mag_squared > err_clip_squared) { error *= std::sqrt(err_clip_squared / mag_squared); myNumErrorClip++; } - // Update the input to the FFT directl to avoid another copy for the + // Update the input to the FFT directly to avoid another copy for the // subsequence IFFT complexf *fft_in = reinterpret_cast(myFftIn); fft_in[i] = constellation_point + error; } + auto minmax = std::minmax_element(error_norm.begin(), error_norm.end()); + etiLog.level(debug) << "err min: " << std::sqrt(*minmax.first) + << " max: " << std::sqrt(*minmax.second); + // Run our error-compensated symbol through the IFFT again fftwf_execute(myFftPlan); // IFFT from myFftIn to myFftOut } -- cgit v1.2.3 From 4574756f50fd34355aa9d6bfc8a5662f87e88661 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Wed, 27 Sep 2017 20:29:14 +0200 Subject: Add CFR statistics to RC --- doc/example.ini | 7 ++--- src/OfdmGenerator.cpp | 71 ++++++++++++++++++++++++++++++++++++++++----------- src/OfdmGenerator.h | 14 +++++++--- 3 files changed, 71 insertions(+), 21 deletions(-) (limited to 'src/OfdmGenerator.cpp') diff --git a/doc/example.ini b/doc/example.ini index 512e196..21f15ea 100644 --- a/doc/example.ini +++ b/doc/example.ini @@ -127,16 +127,17 @@ dac_clk_rate=0 ; and ;dac_clk_rate=128000000 -; Settings for crest factor reduction +; Settings for crest factor reduction. Statistics for ratio of +; samples that were clipped are available through the RC. [cfr] enable=0 ; At what amplitude the signal should be clipped -clip=1.0 +clip=70.0 ; How much to clip the error signal used to compensate the effect ; of clipping -error_clip=1.0 +error_clip=0.05 [firfilter] ; The FIR Filter can be used to create a better spectral quality. diff --git a/src/OfdmGenerator.cpp b/src/OfdmGenerator.cpp index 8f5d8a0..4d76dbc 100644 --- a/src/OfdmGenerator.cpp +++ b/src/OfdmGenerator.cpp @@ -36,7 +36,9 @@ #include #include #include +#include +static const size_t MAX_CLIP_STATS = 10; OfdmGenerator::OfdmGenerator(size_t nbSymbols, size_t nbCarriers, @@ -68,6 +70,7 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols, RC_ADD_PARAMETER(cfr, "Enable crest factor reduction"); RC_ADD_PARAMETER(clip, "CFR: Clip to amplitude"); RC_ADD_PARAMETER(errorclip, "CFR: Limit error"); + RC_ADD_PARAMETER(clip_stats, "CFR: statistics (clip ratio, errorclip ratio)"); if (inverse) { myPosDst = (nbCarriers & 1 ? 0 : 1); @@ -170,13 +173,13 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) "OfdmGenerator::process output size not valid!"); } - myNumClip = 0; - myNumErrorClip = 0; - // It is not guaranteed that fftw keeps the FFT input vector intact. // That's why we copy it to the reference. std::vector reference; + size_t num_clip = 0; + size_t num_error_clip = 0; + for (size_t i = 0; i < myNbSymbols; ++i) { myFftIn[0][0] = 0; myFftIn[0][1] = 0; @@ -196,7 +199,9 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) if (myCfr) { complexf *symbol = reinterpret_cast(myFftOut); - cfr_one_iteration(symbol, reference.data()); + const auto stat = cfr_one_iteration(symbol, reference.data()); + num_clip += stat.clip_count; + num_error_clip += stat.errclip_count; } memcpy(out, myFftOut, mySpacing * sizeof(FFT_TYPE)); @@ -206,19 +211,35 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) } if (myCfr) { - etiLog.level(debug) << "CFR: " << myNumClip << " clipped, " << - myNumErrorClip << " err clipped"; + std::lock_guard lock(myCfrRcMutex); + + const double num_samps = myNbSymbols * mySpacing; + const double clip_ratio = (double)num_clip / num_samps; + + myClipRatios.push_back(clip_ratio); + while (myClipRatios.size() > MAX_CLIP_STATS) { + myClipRatios.pop_front(); + } + + const double errclip_ratio = (double)num_error_clip / num_samps; + myErrorClipRatios.push_back(errclip_ratio); + while (myErrorClipRatios.size() > MAX_CLIP_STATS) { + myErrorClipRatios.pop_front(); + } } return sizeOut; } -void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *reference) +OfdmGenerator::cfr_iter_stat_t OfdmGenerator::cfr_one_iteration( + complexf *symbol, const complexf *reference) { // use std::norm instead of std::abs to avoid calculating the // square roots const float clip_squared = myCfrClip * myCfrClip; + OfdmGenerator::cfr_iter_stat_t ret; + // Clip for (size_t i = 0; i < mySpacing; i++) { const float mag_squared = std::norm(symbol[i]); @@ -228,7 +249,7 @@ void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *referenc // = x * sqrt(clip_squared) / sqrt(mag_squared) // = x * sqrt(clip_squared / mag_squared) symbol[i] *= std::sqrt(clip_squared / mag_squared); - myNumClip++; + ret.clip_count++; } } @@ -245,8 +266,8 @@ void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *referenc std::vector error_norm(mySpacing); for (size_t i = 0; i < mySpacing; i++) { - // FFTW computes an unnormalised trasform, i.e. a FFT-IFFT pair - // or vice-versa give back the original vector scaled by a factor + // FFTW computes an unnormalised transform, i.e. a FFT-IFFT pair + // or vice-versa gives back the original vector scaled by a factor // FFT-size. Because we're comparing our constellation point // (calculated with IFFT-clip-FFT) against reference (input to // the IFFT), we need to divide by our FFT size. @@ -260,7 +281,7 @@ void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *referenc if (mag_squared > err_clip_squared) { error *= std::sqrt(err_clip_squared / mag_squared); - myNumErrorClip++; + ret.errclip_count++; } // Update the input to the FFT directly to avoid another copy for the @@ -269,12 +290,10 @@ void OfdmGenerator::cfr_one_iteration(complexf *symbol, const complexf *referenc fft_in[i] = constellation_point + error; } - auto minmax = std::minmax_element(error_norm.begin(), error_norm.end()); - etiLog.level(debug) << "err min: " << std::sqrt(*minmax.first) - << " max: " << std::sqrt(*minmax.second); - // Run our error-compensated symbol through the IFFT again fftwf_execute(myFftPlan); // IFFT from myFftIn to myFftOut + + return ret; } @@ -294,6 +313,9 @@ void OfdmGenerator::set_parameter(const std::string& parameter, else if (parameter == "errorclip") { ss >> myCfrErrorClip; } + else if (parameter == "clip_stats") { + throw ParameterError("Parameter 'clip_stats' is read-only"); + } else { stringstream ss; ss << "Parameter '" << parameter @@ -315,6 +337,25 @@ const std::string OfdmGenerator::get_parameter(const std::string& parameter) con else if (parameter == "errorclip") { ss << std::fixed << myCfrErrorClip; } + else if (parameter == "clip_stats") { + std::lock_guard lock(myCfrRcMutex); + if (myClipRatios.empty() or myErrorClipRatios.empty()) { + ss << "No stats available"; + } + else { + const double avg_clip_ratio = + std::accumulate(myClipRatios.begin(), myClipRatios.end(), 0.0) / + myClipRatios.size(); + + const double avg_errclip_ratio = + std::accumulate(myErrorClipRatios.begin(), myErrorClipRatios.end(), 0.0) / + myErrorClipRatios.size(); + + ss << "Statistics : " << std::fixed << + avg_clip_ratio * 100 << "%"" samples clipped, " << + avg_errclip_ratio * 100 << "%"" errors clipped"; + } + } else { ss << "Parameter '" << parameter << "' is not exported by controllable " << get_rc_name(); diff --git a/src/OfdmGenerator.h b/src/OfdmGenerator.h index ad39f5b..9ac0387 100644 --- a/src/OfdmGenerator.h +++ b/src/OfdmGenerator.h @@ -68,7 +68,13 @@ class OfdmGenerator : public ModCodec, public RemoteControllable const std::string& parameter) const override; protected: - void cfr_one_iteration(complexf *symbol, const complexf *reference); + struct cfr_iter_stat_t { + size_t clip_count = 0; + size_t errclip_count = 0; + }; + + cfr_iter_stat_t cfr_one_iteration( + complexf *symbol, const complexf *reference); fftwf_plan myFftPlan; fftwf_complex *myFftIn, *myFftOut; @@ -85,14 +91,16 @@ class OfdmGenerator : public ModCodec, public RemoteControllable unsigned myZeroSize; bool myCfr; // Whether to enable crest factor reduction + mutable std::mutex myCfrRcMutex; float myCfrClip; float myCfrErrorClip; fftwf_plan myCfrFft; fftwf_complex *myCfrPostClip; fftwf_complex *myCfrPostFft; - size_t myNumClip; - size_t myNumErrorClip; + // Statistics for CFR + std::deque myClipRatios; + std::deque myErrorClipRatios; }; -- cgit v1.2.3 From 0cfca99b6f0a9b03f148e30c7c44e1f32b82baa1 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Thu, 28 Sep 2017 13:34:53 +0200 Subject: CFR: add some sort of MER calculation that is not correct yet --- src/OfdmGenerator.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++++++++-- src/OfdmGenerator.h | 3 +++ 2 files changed, 51 insertions(+), 2 deletions(-) (limited to 'src/OfdmGenerator.cpp') diff --git a/src/OfdmGenerator.cpp b/src/OfdmGenerator.cpp index 4d76dbc..6a5044e 100644 --- a/src/OfdmGenerator.cpp +++ b/src/OfdmGenerator.cpp @@ -177,9 +177,14 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) // That's why we copy it to the reference. std::vector reference; + std::vector before_cfr; + size_t num_clip = 0; size_t num_error_clip = 0; + // For performance reasons, do not calculate MER for every symbols. + myLastMERCalc = (myLastMERCalc + 1) % myNbSymbols; + for (size_t i = 0; i < myNbSymbols; ++i) { myFftIn[0][0] = 0; myFftIn[0][1] = 0; @@ -198,8 +203,40 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) fftwf_execute(myFftPlan); // IFFT from myFftIn to myFftOut if (myCfr) { + if (myLastMERCalc == i) { + before_cfr.resize(mySpacing); + memcpy(before_cfr.data(), myFftOut, mySpacing * sizeof(FFT_TYPE)); + } + complexf *symbol = reinterpret_cast(myFftOut); + /* cfr_one_iteration runs the myFftPlan again at the end, and + * therefore writes the output data to myFftOut. + */ const auto stat = cfr_one_iteration(symbol, reference.data()); + + if (myLastMERCalc == i) { + /* MER definition, ETSI ETR 290, Annex C + * + * \sum I^2 Q^2 + * MER[dB] = 10 log_10( -------------- ) + * \sum dI^2 dQ^2 + * Where I and Q are the ideal coordinates, and dI and dQ are the errors + * in the received datapoints. + * + * In our case, we consider the constellation points given to the + * OfdmGenerator as "ideal", and we compare the CFR output to it. + */ + + double sum_iq = 0; + double sum_delta = 0; + for (size_t i = 0; i < mySpacing; i++) { + sum_iq += std::norm(before_cfr[i]); + sum_delta += std::norm(symbol[i] - before_cfr[i]); + } + const double mer = 10.0 * std::log10(sum_iq / sum_delta); + myMERs.push_back(mer); + } + num_clip += stat.clip_count; num_error_clip += stat.errclip_count; } @@ -226,6 +263,10 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) while (myErrorClipRatios.size() > MAX_CLIP_STATS) { myErrorClipRatios.pop_front(); } + + while (myMERs.size() > MAX_CLIP_STATS) { + myMERs.pop_front(); + } } return sizeOut; @@ -339,7 +380,7 @@ const std::string OfdmGenerator::get_parameter(const std::string& parameter) con } else if (parameter == "clip_stats") { std::lock_guard lock(myCfrRcMutex); - if (myClipRatios.empty() or myErrorClipRatios.empty()) { + if (myClipRatios.empty() or myErrorClipRatios.empty() or myMERs.empty()) { ss << "No stats available"; } else { @@ -351,9 +392,14 @@ const std::string OfdmGenerator::get_parameter(const std::string& parameter) con std::accumulate(myErrorClipRatios.begin(), myErrorClipRatios.end(), 0.0) / myErrorClipRatios.size(); + const double avg_mer = + std::accumulate(myMERs.begin(), myMERs.end(), 0.0) / + myMERs.size(); + ss << "Statistics : " << std::fixed << avg_clip_ratio * 100 << "%"" samples clipped, " << - avg_errclip_ratio * 100 << "%"" errors clipped"; + avg_errclip_ratio * 100 << "%"" errors clipped. " << + "MER after CFR: " << avg_mer << " dB"; } } else { diff --git a/src/OfdmGenerator.h b/src/OfdmGenerator.h index 9ac0387..f357fa6 100644 --- a/src/OfdmGenerator.h +++ b/src/OfdmGenerator.h @@ -101,6 +101,9 @@ class OfdmGenerator : public ModCodec, public RemoteControllable // Statistics for CFR std::deque myClipRatios; std::deque myErrorClipRatios; + + size_t myLastMERCalc = 0; + std::deque myMERs; }; -- cgit v1.2.3