Determining Error Rate of Phase Shift Keying via Monte Carlo Simulation

Question

My program calculates the symbol error rate for various signal-to-noise ratio (SNR) and modulation orders using Monte Carlo simulation on Phase Shift Keying (PSK). How can it be improved?

main.cpp

#include <iostream>
#include "psksimulator.h"
int main()
{
 constexpr unsigned int order = 16;
 constexpr double snr = 10.0;
 constexpr int threshold = 1000;
 std::cout << PSKSimulator::simulate(order, snr, threshold) << "\n";
 return 0;
}

psksimulator.h

#ifndef PSKSIMULATOR_H
#define PSKSIMULATOR_H
#include <complex>
#include <vector>
using complex_d = std::complex<double>;
using vector_d = std::vector<complex_d>;
// Calculates symbol error rate for PSK
class PSKSimulator
{
public:
 PSKSimulator() = default;
 static double simulate(unsigned int order, double snr_dB, int threshold = 100);
private:
 static bool check_order(unsigned int order);
 static void generate_symbols(std::vector<int>& symbols, unsigned int order);
 static void symbols_to_iq(vector_d& symbols_iq, const std::vector<int>& symbols, unsigned int order);
 static void generate_noise(vector_d& noise, double snr_dB);
 static int demodulate_and_check(const vector_d& received, const std::vector<int>& symbols, unsigned int order);
 static double wrap_angle(double angle);
};
#endif // PSKSIMULATOR_H

psksimulator.cpp

#include "psksimulator.h"
#include <cstddef>
#include <iostream>
#include <algorithm>
#include <random>
#include <cmath>
double PSKSimulator::simulate(const unsigned int order, double snr_dB, const int threshold)
{
 if(!check_order(order))
 {
 return 1.0;
 }
 constexpr size_t vector_size = 10'000;
 std::vector<int> symbols(vector_size);
 std::vector<std::complex<double>> symbols_iq(vector_size);
 std::vector<std::complex<double>> noise(vector_size);
 std::vector<std::complex<double>> received(vector_size);
 int number_of_errors = 0;
 int number_of_iter = 0;
 while(number_of_errors < threshold)
 {
 generate_symbols(symbols, order);
 symbols_to_iq(symbols_iq, symbols, order);
 generate_noise(noise, snr_dB);
 std::transform(symbols_iq.begin(), symbols_iq.end(), noise.begin(), received.begin(), std::plus<complex_d>());
 number_of_errors += demodulate_and_check(received, symbols, order);
 number_of_iter++;
 }
 return static_cast<double>(number_of_errors) / (number_of_iter* vector_size);
}
bool PSKSimulator::check_order(unsigned int order)
{
 if((order < 2) || order > 256)
 {
 std::cerr << "Modulation order must be between 2 and 256!\n";
 return false;
 }
 else if(!(order & (order - 1)) == 0) // http://www.graphics.stanford.edu/~seander/bithacks.html
 {
 std::cerr << "Modulation order must be power of 2!\n";
 return false;
 }
 else
 {
 return true;
 }
}
void PSKSimulator::generate_symbols(std::vector<int>& symbols, unsigned int order)
{
 std::random_device random_device;
 std::mt19937 generator(random_device());
 std::uniform_int_distribution<int> distr(0, order - 1);
 for(auto & item : symbols)
 {
 item = distr(generator);
 }
}
void PSKSimulator::symbols_to_iq(vector_d& symbols_iq, const std::vector<int>& symbols, unsigned int order)
{
 const int nonbinary = (order == 2) ? 0 : 1;
 for(int idx = 0; idx < symbols.size(); idx++)
 {
 symbols_iq[idx] = std::polar<double>(1, nonbinary * (M_PI / order) + symbols[idx] * (2 * M_PI / order));
 }
}
void PSKSimulator::generate_noise(vector_d& noise, double snr_dB)
{
 constexpr double mean = 0.0;
 const double variance = (1 / std::sqrt(2)) * std::pow(10.0, -snr_dB / 20);
 std::random_device random_device{};
 std::mt19937 generator{random_device()};
 std::normal_distribution distr{mean, variance};
 for(auto & item : noise)
 {
 item = std::complex<double>(distr(generator), distr(generator));
 }
}
int PSKSimulator::demodulate_and_check(const vector_d& received, const std::vector<int>& symbols, unsigned int order)
{
 int error = 0;
 const bool is_bpsk = (order == 2);
 for(int idx = 0; idx < received.size(); idx++)
 {
 int decided_symbol;
 if(is_bpsk)
 {
 if(received[idx].real() > 0)
 {
 decided_symbol = 0;
 }
 else
 {
 decided_symbol = 1;
 }
 }
 else
 {
 const auto angle = wrap_angle(std::arg(received[idx]));
 for(int angle_idx = 0; angle_idx < order; angle_idx++)
 {
 if((angle > (angle_idx * 2 * M_PI / order)) && (angle < ((angle_idx + 1) * 2 * M_PI / order)))
 {
 decided_symbol = angle_idx;
 }
 }
 }
 if(decided_symbol != symbols[idx])
 {
 error++;
 }
 }
 return error;
}
double PSKSimulator::wrap_angle(double angle)
{
 return (angle > 0) ? angle : (2 * M_PI + angle);
}

Answer 1

Add a way to pass simulation parameters via the command line

You are hardcoding the test parameters in main(). That means you have to change that file and recompile the code before you can simulate with different parameters. This gets old really quick. A simple solution is to just use int main(int argc, char* argv[]) and use std::stoi() and std::stod() to convert the parameters to integers and doubles.

Declare type aliases inside PSKSimulator

Using type aliases can help avoid typing long type names everywhere, and also makes it easier to change the types. However, complex_d and vector_d are only used inside PSKSimulator, so declare those type aliases inside that class. This avoids poluting the global namespace.

Unnecessary use of a class

class PSKSimulator has no member variables and only static member functions. (Well, technically you declared a constructor which is a non-static member, but it's not doing anything so it can be removed.) That means there is no real benefit to making this a class. Instead, the static member functions can just become regular functions. If you want, you can put them in a namespace instead.

Unnecessary temporary storage

You are creating a lot of std::vectors, each holding the result of one stage in the simulation. But this is totally unnecessary; you could just make each function process a single symbol, which is equivalent to setting vector_size to 1. This avoids a lot of memory bandwidth, and keeps everything in the lowest level cache of the CPU, and it will also simplify the code.

This problem can be solved analytically

Monte Carlo simulations are very useful, but the problem is often that you need a lot of iterations to get a reasonably accurate result. If the noise level is very small, it will take a long time to reach the threshold. In this case though, there is a much better way to solve the problem; you can just analytically calculate what the error rate is. This will give you an exact answer in one step.

Return a proper error if the order is incorrect

If check_order() returns false, your simulate() returns 1.0. However, this could also be a valid result for a simulation with a very high noise level. You do print an error message to std::cerr, but this is not always visible (consider your program being called from a script running in the background, maybe on a remote cluster). Make sure you also exit with a proper error code, or equally effective, throw one of the standard exceptions.

Why limit the order?

You are limiting the order to be a power of two between 2 and 256. However, there is no reason why it has to be this strict. Any integer larger than 1 can be simulated by your code.

Naming things

Make sure the names of functions, classes and variables are clear, precise and consistent. This will avoid errors and makes it easier for everyone, including you in the future, to understand what the code is doing. Here are some issues I see:

• snr: this is inconsistent, elsewhere you use snr_dB.
• threshold: threshold of what? This is not a precise name. Something better might be error_threshold or minimum_error_count.
• distr: avoid unnessary abbreviations, just write distribution.
• vector_d, complex_d: the _d is an unfortunate abbreviation. You really have to look up the definition to find out what it means. And if you expand it, these type aliases are not so useful anymore. But then you are not using them consistently anyway. I recommend just avoiding them altogether.
• is_bpsk: again an unnecessary abbreviation, and inconsistent with nonbinary used elsewhere.

Answer 2

Sounds like a cool project!

I just have some basic things to point out:

1. Consider removing the space between auto and & in your for loops for typing consistency. Everywhere else, you have no white space before the &.
2. Consider moving your call to received.size() (and similar calls) out of the for loop header and saving the output into a variable. If you don't, that method will be called every iteration.
3. You never change the value of snr_dB or order, so consider making them const.
4. For consistency and general good practice, consider initializing all objects. You initialize random_device in one spot, but not the other.

I will edit this answer if I find anything else to bring up.

• unique: Thank you, I had a misunderstanding that compiler optimizes the size function in loop condition if the vector passed as const. Now, I learned that is not the case. Link