diff --git a/CMakeLists.txt b/CMakeLists.txt index e840fcb..677e5a8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,6 +3,7 @@ project(MILSTD110C) # Set C++17 standard set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) # Include directories include_directories(include) @@ -17,17 +18,43 @@ add_subdirectory(tests) # Set source files set(SOURCES main.cpp) -# Link with libsndfile +# Find required packages list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake") find_package(SndFile REQUIRED) find_package(FFTW3 REQUIRED) +find_package(fmt REQUIRED) +find_package(Gnuradio REQUIRED COMPONENTS + analog + blocks + channels + filter + fft + runtime +) + +if(NOT Gnuradio_FOUND) + message(FATAL_ERROR "GNU Radio not found!") +endif() + +# Include GNU Radio directories +include_directories(${Gnuradio_INCLUDE_DIRS}) +link_directories(${Gnuradio_LIBRARY_DIRS}) # Add executable add_executable(MILSTD110C ${SOURCES}) -# Link executable with libsndfile library -target_link_libraries(MILSTD110C SndFile::sndfile) -target_link_libraries(MILSTD110C FFTW3::fftw3) +# Link executable with required libraries +target_link_libraries(MILSTD110C + SndFile::sndfile + FFTW3::fftw3 + gnuradio::gnuradio-runtime + gnuradio::gnuradio-analog + gnuradio::gnuradio-blocks + gnuradio::gnuradio-filter + gnuradio::gnuradio-fft + gnuradio::gnuradio-channels + fmt::fmt +) # Debug and Release Build Types set(CMAKE_CONFIGURATION_TYPES "Debug;Release" CACHE STRING "" FORCE) diff --git a/include/encoder/ModemController.h b/include/encoder/ModemController.h index 48bfcca..0ec062e 100644 --- a/include/encoder/ModemController.h +++ b/include/encoder/ModemController.h @@ -57,7 +57,7 @@ public: * @return The scrambled data ready for modulation. * @note The modulated signal is generated internally but is intended to be handled externally. */ - std::vector transmit(BitStream input_data) { + std::vector transmit(const BitStream& input_data) { // Step 1: Append EOM Symbols BitStream eom_appended_data = appendEOMSymbols(input_data); diff --git a/include/encoder/SymbolFormation.h b/include/encoder/SymbolFormation.h index f5a4068..a68c893 100644 --- a/include/encoder/SymbolFormation.h +++ b/include/encoder/SymbolFormation.h @@ -43,37 +43,33 @@ class SymbolFormation { // Generate and scramble the sync preamble std::vector sync_preamble = generateSyncPreamble(); sync_preamble = scrambler.scrambleSyncPreamble(sync_preamble); - - size_t set_count = 0; - size_t symbol_count = 0; - size_t current_frame = 0; + std::vector data_stream; - size_t current_index = 0; - while (current_index < symbol_data.size()) { - // Determine the size of the current unknown data block - size_t block_size = std::min(unknown_data_block_size, symbol_data.size() - current_index); - std::vector unknown_data_block(symbol_data.begin() + current_index, symbol_data.begin() + current_index + block_size); - current_index += block_size; + if (baud_rate == 75) { + size_t set_count = 0; + for (size_t i = 0; i < symbol_data.size(); i++) { + bool is_exceptional_set = (set_count % ((interleave_setting == 1) ? 45 : 360)) == 0; + append75bpsMapping(data_stream, symbol_data[i], is_exceptional_set); + set_count++; + } + } else { + size_t symbol_count = 0; + size_t current_frame = 0; + size_t current_index = 0; - // Map the unknown data based on baud rate - if (baud_rate == 75) { - size_t set_size = (interleave_setting == 2) ? 360 : 32; - for (size_t i = 0; i < unknown_data_block.size(); i += set_size) { - bool is_exceptional_set = (set_count % ((interleave_setting == 1) ? 45 : 360)) == 0; - std::vector mapped_set = map75bpsSet(unknown_data_block, i, set_size, is_exceptional_set); - data_stream.insert(data_stream.end(), mapped_set.begin(), mapped_set.end()); - set_count++; - } - } else { - // For baud rates greater than 75 bps + while (current_index < symbol_data.size()) { + // Determine the size of the current unknown data block + size_t block_size = std::min(unknown_data_block_size, symbol_data.size() - current_index); + std::vector unknown_data_block(symbol_data.begin() + current_index, symbol_data.begin() + current_index + block_size); + current_index += block_size; + + // Map the unknown data based on baud rate std::vector mapped_unknown_data = mapUnknownData(unknown_data_block); symbol_count += mapped_unknown_data.size(); data_stream.insert(data_stream.end(), mapped_unknown_data.begin(), mapped_unknown_data.end()); - } - // Insert probe data if we are at an interleaver block boundary - if (baud_rate > 75) { + // Insert probe data if we are at an interleaver block boundary std::vector probe_data = generateProbeData(current_frame, total_frames); data_stream.insert(data_stream.end(), probe_data.begin(), probe_data.end()); current_frame = (current_frame + 1) % total_frames; @@ -134,17 +130,14 @@ class SymbolFormation { throw std::invalid_argument("Invalid channel symbol"); } - if (repeat_twice) { - // Repeat the pattern twice instead of four times for known symbols - tribit_pattern.insert(tribit_pattern.end(), tribit_pattern.begin(), tribit_pattern.end()); - } else { - // Repeat the pattern four times as per Table XIII - tribit_pattern.insert(tribit_pattern.end(), tribit_pattern.begin(), tribit_pattern.end()); - tribit_pattern.insert(tribit_pattern.end(), tribit_pattern.begin(), tribit_pattern.end()); - tribit_pattern.insert(tribit_pattern.end(), tribit_pattern.begin(), tribit_pattern.end()); + size_t repetitions = repeat_twice ? 2 : 4; + std::vector repeated_pattern; + + for (size_t i = 0; i < repetitions; i++) { + repeated_pattern.insert(repeated_pattern.end(), tribit_pattern.begin(), tribit_pattern.end()); } - return tribit_pattern; + return repeated_pattern; } std::vector generateSyncPreamble() { @@ -340,19 +333,6 @@ class SymbolFormation { } } - std::vector map75bpsSet(const std::vector& data, size_t start_index, size_t set_size, bool is_exceptional_set) { - std::vector mapped_set; - - // Make sure we do not exceed the size of the data vector - size_t end_index = std::min(start_index + set_size, data.size()); - - for (size_t i = start_index; i < end_index; ++i) { - append75bpsMapping(mapped_set, data[i], is_exceptional_set); - } - - return mapped_set; - } - std::vector mapUnknownData(const std::vector& data) { std::vector mapped_data; diff --git a/include/modulation/PSKModulator.h b/include/modulation/PSKModulator.h index 4ed7b0c..0fd1824 100644 --- a/include/modulation/PSKModulator.h +++ b/include/modulation/PSKModulator.h @@ -26,14 +26,12 @@ static constexpr double SCALE_FACTOR = 32767.0; class PSKModulator { public: PSKModulator(const double _sample_rate, const bool _is_frequency_hopping, const size_t _num_taps) - : sample_rate(validateSampleRate(_sample_rate)), gain(1.0/sqrt(2.0)), is_frequency_hopping(_is_frequency_hopping), samples_per_symbol(static_cast(sample_rate / SYMBOL_RATE)), srrc_filter(48, _sample_rate, SYMBOL_RATE, ROLLOFF_FACTOR) { + : sample_rate(validateSampleRate(_sample_rate)), gain(1.0/sqrt(2.0)), is_frequency_hopping(_is_frequency_hopping), samples_per_symbol(static_cast(sample_rate / SYMBOL_RATE)), srrc_filter(8, _sample_rate, SYMBOL_RATE, ROLLOFF_FACTOR) { initializeSymbolMap(); phase_detector = PhaseDetector(symbolMap); } std::vector modulate(const std::vector& symbols) { - std::vector baseband_I(symbols.size() * samples_per_symbol); - std::vector baseband_Q(symbols.size() * samples_per_symbol); std::vector> baseband_components(symbols.size() * samples_per_symbol); size_t symbol_index = 0; @@ -59,11 +57,11 @@ public: double carrier_phase = 0.0; double carrier_phase_increment = 2 * M_PI * CARRIER_FREQ / sample_rate; - for (const auto& sample : baseband_components) { + for (const auto& sample : filtered_components) { double carrier_cos = std::cos(carrier_phase); double carrier_sin = -std::sin(carrier_phase); double passband_value = sample.real() * carrier_cos + sample.imag() * carrier_sin; - passband_signal.emplace_back(passband_value * 32767.0); // Scale to int16_t + passband_signal.emplace_back(passband_value * SCALE_FACTOR); // Scale to int16_t carrier_phase += carrier_phase_increment; if (carrier_phase >= 2 * M_PI) carrier_phase -= 2 * M_PI; @@ -117,6 +115,7 @@ public: return processDataSymbols(detected_symbols); } + return std::vector(); } private: @@ -130,7 +129,7 @@ private: static inline double validateSampleRate(const double rate) { - if (rate <= 2 * CARRIER_FREQ) { + if (rate <= 2 * (CARRIER_FREQ + SYMBOL_RATE * (1 + ROLLOFF_FACTOR) / 2)) { throw std::out_of_range("Sample rate must be above the Nyquist frequency (PSKModulator.h)"); } return rate; diff --git a/include/utils/filters.h b/include/utils/filters.h index 96fe140..35dfb76 100644 --- a/include/utils/filters.h +++ b/include/utils/filters.h @@ -9,92 +9,64 @@ class TapGenerators { public: - std::vector generateSRRCTaps(const size_t num_taps, const double sample_rate, const double symbol_rate, const double rolloff) const { - std::vector freq_response(num_taps, 0.0); + std::vector generateSRRCTaps(size_t num_taps, double sample_rate, double symbol_rate, double rolloff) const { std::vector taps(num_taps); + double T = 1.0 / symbol_rate; // Symbol period + double dt = 1.0 / sample_rate; // Time step + double t_center = (num_taps - 1) / 2.0; - double fn = symbol_rate / 2.0; - double f_step = sample_rate / num_taps; + for (size_t i = 0; i < num_taps; ++i) { + double t = (i - t_center) * dt; + double sinc_part = (t == 0.0) ? 1.0 : std::sin(M_PI * t / T * (1 - rolloff)) / (M_PI * t / T * (1 - rolloff)); + double cos_part = (t == 0.0) ? std::cos(M_PI * t / T * (1 + rolloff)) : std::cos(M_PI * t / T * (1 + rolloff)); + double denominator = 1.0 - (4.0 * rolloff * t / T) * (4.0 * rolloff * t / T); - for (size_t i = 0; i < num_taps / 2; i++) { - double f = i * f_step; - - if (f <= fn * (1 - rolloff)) { - freq_response[i] = 1.0; - } else if (f <= fn * (1 + rolloff)) { - freq_response[i] = 0.5 * (1 - std::sin(M_PI * (f - fn * (1 - rolloff)) / (2 * rolloff * fn))); + if (std::fabs(denominator) < 1e-8) { + // Handle singularity at t = T / (4R) + taps[i] = rolloff * (std::sin(M_PI / (4.0 * rolloff)) + (1.0 / (4.0 * rolloff)) * std::cos(M_PI / (4.0 * rolloff))) / (M_PI / (4.0 * rolloff)); } else { - freq_response[i] = 0.0; + taps[i] = (4.0 * rolloff / (M_PI * std::sqrt(T))) * (cos_part / denominator); } + + taps[i] *= sinc_part; } - for (size_t i = num_taps / 2; i < num_taps; i++) { - freq_response[i] = freq_response[num_taps - i - 1]; - } - - fftw_complex* freq_domain = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * num_taps); - for (size_t i = 0; i < num_taps; i++) { - freq_domain[i][0] = freq_response[i]; - freq_domain[i][1] = 0.0; - } - - std::vector time_domain_taps(num_taps, 0.0); - - fftw_plan plan = fftw_plan_dft_c2r_1d(num_taps, freq_domain, time_domain_taps.data(), FFTW_ESTIMATE); - fftw_execute(plan); - fftw_destroy_plan(plan); - fftw_free(freq_domain); - - double norm_factor = std::sqrt(std::accumulate(time_domain_taps.begin(), time_domain_taps.end(), 0.0, [](double sum, double val) { return sum + val * val; })); - for (auto& tap : time_domain_taps) { - tap /= norm_factor; - } - - return time_domain_taps; - } - - std::vector generateLowpassTaps(const size_t num_taps, const double cutoff_freq, const double sample_rate) const { - std::vector freq_response(num_taps, 0.0); - std::vector taps(num_taps); - - double fn = cutoff_freq / 2.0; - double f_step = sample_rate / num_taps; - - // Define frequency response - for (size_t i = 0; i < num_taps / 2; i++) { - double f = i * f_step; - if (f <= fn) { - freq_response[i] = 1.0; // Passband - } else { - freq_response[i] = 0.0; // Stopband - } - } - - // Mirror the second half of the response for symmetry - for (size_t i = num_taps / 2; i < num_taps; i++) { - freq_response[i] = freq_response[num_taps - i - 1]; - } - - // Perform inverse FFT to get time-domain taps - fftw_complex* freq_domain = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * num_taps); - for (size_t i = 0; i < num_taps; i++) { - freq_domain[i][0] = freq_response[i]; - freq_domain[i][1] = 0.0; - } - - std::vector time_domain_taps(num_taps, 0.0); - fftw_plan plan = fftw_plan_dft_c2r_1d(num_taps, freq_domain, time_domain_taps.data(), FFTW_ESTIMATE); - fftw_execute(plan); - fftw_destroy_plan(plan); - fftw_free(freq_domain); - // Normalize filter taps - double norm_factor = std::sqrt(std::accumulate(time_domain_taps.begin(), time_domain_taps.end(), 0.0, [](double sum, double val) { return sum + val * val; })); - for (auto& tap : time_domain_taps) { - tap /= norm_factor; + double sum = std::accumulate(taps.begin(), taps.end(), 0.0); + for (auto& tap : taps) { + tap /= sum; } - return time_domain_taps; + return taps; + } + + std::vector generateLowpassTaps(size_t num_taps, double cutoff_freq, double sample_rate) const { + std::vector taps(num_taps); + double fc = cutoff_freq / (sample_rate / 2.0); // Normalized cutoff frequency (0 < fc < 1) + double M = num_taps - 1; + double mid = M / 2.0; + + for (size_t n = 0; n < num_taps; ++n) { + double n_minus_mid = n - mid; + double h; + if (n_minus_mid == 0.0) { + h = fc; + } else { + h = fc * (std::sin(M_PI * fc * n_minus_mid) / (M_PI * fc * n_minus_mid)); + } + + // Apply window function (e.g., Hamming window) + double window = 0.54 - 0.46 * std::cos(2.0 * M_PI * n / M); + taps[n] = h * window; + } + + // Normalize filter taps + double sum = std::accumulate(taps.begin(), taps.end(), 0.0); + for (auto& tap : taps) { + tap /= sum; + } + + return taps; } }; @@ -103,26 +75,38 @@ class Filter { Filter(const std::vector& _filter_taps) : filter_taps(_filter_taps), buffer(_filter_taps.size(), 0.0), buffer_index(0) {} double filterSample(const double sample) { - buffer[buffer_index] = std::complex(sample,0.0); + buffer[buffer_index] = std::complex(sample, 0.0); double filtered_val = 0.0; + size_t idx = buffer_index; + for (size_t j = 0; j < filter_taps.size(); j++) { - size_t signal_index = (buffer_index + j) % filter_taps.size(); - filtered_val += filter_taps[j] * buffer[signal_index].real(); + filtered_val += filter_taps[j] * buffer[idx].real(); + if (idx == 0) { + idx = buffer.size() - 1; + } else { + idx--; + } } - buffer_index = (buffer_index + 1) % filter_taps.size(); + buffer_index = (buffer_index + 1) % buffer.size(); return filtered_val; } - std::complex filterSample(const std::complex& sample) { + std::complex filterSample(const std::complex sample) { buffer[buffer_index] = sample; std::complex filtered_val = std::complex(0.0, 0.0); + size_t idx = buffer_index; + for (size_t j = 0; j < filter_taps.size(); j++) { - size_t signal_index = (buffer_index + j) % filter_taps.size(); - filtered_val += filter_taps[j] * buffer[signal_index]; + filtered_val += filter_taps[j] * buffer[idx]; + if (idx == 0) { + idx = buffer.size() - 1; + } else { + idx--; + } } - - buffer_index = (buffer_index + 1) % filter_taps.size(); + + buffer_index = (buffer_index + 1) % buffer.size(); return filtered_val; } diff --git a/include/utils/wattersonchannel.h b/include/utils/wattersonchannel.h new file mode 100644 index 0000000..51a6ab4 --- /dev/null +++ b/include/utils/wattersonchannel.h @@ -0,0 +1,410 @@ +#ifndef WATTERSONCHANNEL_H +#define WATTERSONCHANNEL_H + +#include +#include +#include +#include +#include +#include +#include +#include // FFTW library for FFT-based Hilbert transform + +constexpr double PI = 3.14159265358979323846; + +class WattersonChannel { +public: + WattersonChannel(double sampleRate, double symbolRate, double delaySpread, double fadingBandwidth, double SNRdB, int numSamples, int numpaths, bool isFading); + + // Process a block of input samples + void process(const std::vector& inputSignal, std::vector& outputSignal); + +private: + double Fs; // Sample rate + double Rs; // Symbol rate + double delaySpread; // Delay spread in seconds + std::vector delays = {0, L}; + double fadingBandwidth; // Fading bandwidth d in Hz + double SNRdB; // SNR in dB + int L; // Length of the simulated channel + std::vector f_jt; // Filter impulse response + std::vector>> h_j; // Fading tap gains over time for h_0 and h_(L-1) + double Ts; // Sample period + double k; // Normalization constant for filter + double tau; // Truncation width + double fadingSampleRate; // Sample rate for fading process + std::vector> wgnFadingReal; // WGN samples for fading (double part) + std::vector> wgnFadingImag; // WGN samples for fading (imaginary part) + std::vector> n_i; // WGN samples for noise + std::mt19937 rng; // Random number generator + int numSamples; // Number of samples in the simulation + int numFadingSamples; // Number of fading samples + + int numPaths; + bool isFading; + + void normalizeTapGains(); + void generateFilter(); + void generateFadingTapGains(); + void generateNoise(const std::vector>& x_i); + void generateWGN(std::vector& wgn, int numSamples); + void resampleFadingTapGains(); + void hilbertTransform(const std::vector& input, std::vector>& output); +}; + +WattersonChannel::WattersonChannel(double sampleRate, double symbolRate, double delaySpread, double fadingBandwidth, double SNRdB, int numSamples, int numPaths, bool isFading) + : Fs(sampleRate), Rs(symbolRate), delaySpread(delaySpread), fadingBandwidth(fadingBandwidth), SNRdB(SNRdB), numSamples(numSamples), rng(std::random_device{}()), numPaths(numPaths), isFading(isFading) +{ + Ts = 1.0 / Fs; + // Compute L + if (numPaths == 1) { + L = 1; + } else { + L = static_cast(std::round(delaySpread / Ts)); + if (L < 1) L = 1; + } + + // Compute truncation width tau + double ln100 = std::log(100.0); + tau = std::sqrt(ln100) / (PI * fadingBandwidth); + + // Initialize k (will be normalized later) + k = 1.0; + + // Fading sample rate, at least 32 times the fading bandwidth + fadingSampleRate = std::max(32.0 * fadingBandwidth, Fs); + + h_j.resize(numPaths); + wgnFadingReal.resize(numPaths); + wgnFadingImag.resize(numPaths); + + if (isFading) { + // Generate filter impulse response + generateFilter(); + + // Number of fading samples + double simulationTime = numSamples / Fs; + numFadingSamples = static_cast(std::ceil(simulationTime * fadingSampleRate)); + + // Generate WGN for fading + for (int pathIndex = 0; pathIndex < numPaths; ++pathIndex) { + generateWGN(wgnFadingReal[pathIndex], numFadingSamples); + generateWGN(wgnFadingImag[pathIndex], numFadingSamples); + } + + // Generate fading tap gains + generateFadingTapGains(); + + // Resample fading tap gains to match sample rate Fs + resampleFadingTapGains(); + } else { + // For fixed channel, set tap gains directly + generateFadingTapGains(); + } + + // Generate noise n_i +} + +void WattersonChannel::normalizeTapGains() { + double totalPower = 0.0; + int numValidSamples = h_j[0].size(); + for (int i = 0; i < numValidSamples; i++) { + for (int pathIndex = 0; pathIndex < numPaths; pathIndex++) { + totalPower += std::norm(h_j[pathIndex][i]); + } + } + totalPower /= numValidSamples; + + double normFactor = 1.0 / std::sqrt(totalPower); + for (int pathIndex = 0; pathIndex < numPaths; pathIndex++) { + for (auto& val : h_j[pathIndex]) { + val *= normFactor; + } + } +} + +void WattersonChannel::generateFilter() +{ + // Generate filter impulse response f_j(t) = k * sqrt(2) * e^{-π² * t² * d²}, -tau < t < tau + + // Number of filter samples + int numFilterSamples = static_cast(std::ceil(2 * tau * fadingSampleRate)) + 1; // Include center point + f_jt.resize(numFilterSamples); + + double dt = 1.0 / fadingSampleRate; + int halfSamples = numFilterSamples / 2; + + double totalEnergy = 0.0; + + for (int n = 0; n < numFilterSamples; ++n) { + double t = (n - halfSamples) * dt; + double val = k * std::sqrt(2.0) * std::exp(-PI * PI * t * t * fadingBandwidth * fadingBandwidth); + f_jt[n] = val; + totalEnergy += val * val * dt; + } + + // Normalize k so that total energy is 1.0 + double k_new = k / std::sqrt(totalEnergy); + for (auto& val : f_jt) { + val *= k_new; + } + k = k_new; +} + +void WattersonChannel::generateFadingTapGains() +{ + if (!isFading) { + for (int pathIndex = 0; pathIndex < numPaths; pathIndex++) { + h_j[pathIndex].assign(numSamples, std::complex(1.0, 0.0)); + } + } else { + // Prepare for FFT-based convolution + int convSize = numFadingSamples + f_jt.size() - 1; + int fftSize = 1; + while (fftSize < convSize) { + fftSize <<= 1; // Next power of two + } + + std::vector f_jtPadded(fftSize, 0.0); + std::copy(f_jt.begin(), f_jt.end(), f_jtPadded.begin()); + + fftw_complex* f_jtFFT = fftw_alloc_complex(fftSize); + fftw_plan planF_jt = fftw_plan_dft_r2c_1d(fftSize, f_jtPadded.data(), f_jtFFT, FFTW_ESTIMATE); + fftw_execute(planF_jt); + + for (int pathIndex = 0; pathIndex < numPaths; pathIndex++) { + // Zero-pad inputs + std::vector wgnRealPadded(fftSize, 0.0); + std::vector wgnImagPadded(fftSize, 0.0); + + std::copy(wgnFadingReal[pathIndex].begin(), wgnFadingReal[pathIndex].end(), wgnRealPadded.begin()); + std::copy(wgnFadingImag[pathIndex].begin(), wgnFadingImag[pathIndex].end(), wgnImagPadded.begin()); + + // Perform FFTs + fftw_complex* WGNRealFFT = fftw_alloc_complex(fftSize); + fftw_complex* WGNImagFFT = fftw_alloc_complex(fftSize); + fftw_plan planWGNReal = fftw_plan_dft_r2c_1d(fftSize, wgnRealPadded.data(), WGNRealFFT, FFTW_ESTIMATE); + fftw_plan planWGNImag = fftw_plan_dft_r2c_1d(fftSize, wgnImagPadded.data(), WGNImagFFT, FFTW_ESTIMATE); + + fftw_execute(planWGNReal); + fftw_execute(planWGNImag); + + // Multiply in frequency domain + int fftComplexSize = fftSize / 2 + 1; + for (int i = 0; i < fftComplexSize; ++i) { + // Multiply WGNRealFFT and f_jtFFT + double realPart = WGNRealFFT[i][0] * f_jtFFT[i][0] - WGNRealFFT[i][1] * f_jtFFT[i][1]; + double imagPart = WGNRealFFT[i][0] * f_jtFFT[i][1] + WGNRealFFT[i][1] * f_jtFFT[i][0]; + WGNRealFFT[i][0] = realPart; + WGNRealFFT[i][1] = imagPart; + + // Multiply WGNImagFFT and f_jtFFT + realPart = WGNImagFFT[i][0] * f_jtFFT[i][0] - WGNImagFFT[i][1] * f_jtFFT[i][1]; + imagPart = WGNImagFFT[i][0] * f_jtFFT[i][1] + WGNImagFFT[i][1] * f_jtFFT[i][0]; + WGNImagFFT[i][0] = realPart; + WGNImagFFT[i][1] = imagPart; + } + + // Perform inverse FFTs + fftw_plan planInvReal = fftw_plan_dft_c2r_1d(fftSize, WGNRealFFT, wgnRealPadded.data(), FFTW_ESTIMATE); + fftw_plan planInvImag = fftw_plan_dft_c2r_1d(fftSize, WGNImagFFT, wgnImagPadded.data(), FFTW_ESTIMATE); + + fftw_execute(planInvReal); + fftw_execute(planInvImag); + + // Normalize + double scale = 1.0 / fftSize; + for (int i = 0; i < fftSize; ++i) { + wgnRealPadded[i] *= scale; + wgnImagPadded[i] *= scale; + } + + // Assign h_j[0] and h_j[1] + int numValidSamples = numFadingSamples; + + h_j[pathIndex].resize(numValidSamples); + for (int i = 0; i < numValidSamples; i++) { + h_j[pathIndex][i] = std::complex(wgnRealPadded[i], wgnImagPadded[i]); + } + + // Clean up + fftw_destroy_plan(planWGNReal); + fftw_destroy_plan(planWGNImag); + fftw_destroy_plan(planInvReal); + fftw_destroy_plan(planInvImag); + fftw_free(WGNRealFFT); + fftw_free(WGNImagFFT); + } + + fftw_destroy_plan(planF_jt); + fftw_free(f_jtFFT); + + normalizeTapGains(); + } +} + +void WattersonChannel::resampleFadingTapGains() +{ + // Resample h_j[0] and h_j[1] from fadingSampleRate to Fs + int numOutputSamples = numSamples; + double resampleRatio = fadingSampleRate / Fs; + + for (int pathIndex = 0; pathIndex < numPaths; pathIndex++) { + std::vector> resampled_h(numOutputSamples); + for (int i = 0; i < numOutputSamples; ++i) { + double t = i * (1.0 / Fs); + double index = t * fadingSampleRate; + int idx = static_cast(index); + double frac = index - idx; + + // Simple linear interpolation + if (idx + 1 < h_j[pathIndex].size()) { + resampled_h[i] = h_j[pathIndex][idx] * (1.0 - frac) + h_j[pathIndex][idx + 1] * frac; + } + else if (idx < h_j[pathIndex].size()) { + resampled_h[i] = h_j[pathIndex][idx]; + } + else { + resampled_h[i] = std::complex(0.0, 0.0); + } + } + h_j[pathIndex] = std::move(resampled_h); + } +} + +void WattersonChannel::generateNoise(const std::vector>& x_i) +{ + // Generate WGN samples for noise n_i with appropriate power to achieve the specified SNR + n_i.resize(numSamples); + + double inputSignalPower = 0.0; + for (const auto& sample : x_i) { + inputSignalPower += std::norm(sample); + } + inputSignalPower /= x_i.size(); + + // Compute signal power (assuming average power of input signal x_i is normalized to 1.0) + double channelGainPower = 0.0; + for (int i = 0; i < numSamples; i++) { + std::complex combinedGain = std::complex(0.0, 0.0); + for (int pathIndex = 0; pathIndex < numPaths; pathIndex++) { + combinedGain += h_j[pathIndex][i]; + } + channelGainPower += std::norm(combinedGain); + } + channelGainPower /= numSamples; + + double signalPower = inputSignalPower * channelGainPower; + + // Compute noise power + double SNR_linear = std::pow(10.0, SNRdB / 10.0); + double noisePower = signalPower / SNR_linear; + + std::normal_distribution normalDist(0.0, std::sqrt(noisePower / 2.0)); // Divided by 2 for double and imag parts + + for (int i = 0; i < numSamples; ++i) { + double realPart = normalDist(rng); + double imagPart = normalDist(rng); + n_i[i] = std::complex(realPart, imagPart); + } +} + +void WattersonChannel::generateWGN(std::vector& wgn, int numSamples) +{ + wgn.resize(numSamples); + + std::normal_distribution normalDist(0.0, 1.0); // Standard normal distribution + + for (int i = 0; i < numSamples; ++i) { + wgn[i] = normalDist(rng); + } +} + +void WattersonChannel::hilbertTransform(const std::vector& input, std::vector>& output) +{ + // Implement Hilbert transform using FFT method + int N = input.size(); + + // Allocate input and output arrays for FFTW + double* in = fftw_alloc_real(N); + fftw_complex* out = fftw_alloc_complex(N); + + // Copy input signal to in array + for (int i = 0; i < N; ++i) { + in[i] = input[i]; + } + + // Create plan for forward FFT + fftw_plan plan_forward = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE); + + // Execute forward FFT + fftw_execute(plan_forward); + + // Apply the Hilbert transform in frequency domain + // For positive frequencies, multiply by 2; for zero and negative frequencies, set to zero + int N_half = N / 2 + 1; + for (int i = 0; i < N_half; ++i) { + if (i == 0 || i == N / 2) { // DC and Nyquist frequency components + out[i][0] = 0.0; + out[i][1] = 0.0; + } + else { + out[i][0] *= 2.0; + out[i][1] *= 2.0; + } + } + + // Create plan for inverse FFT + fftw_plan plan_backward = fftw_plan_dft_c2r_1d(N, out, in, FFTW_ESTIMATE); + + // Execute inverse FFT + fftw_execute(plan_backward); + + // Normalize and store result in output vector + output.resize(N); + double scale = 1.0 / N; + for (int i = 0; i < N; ++i) { + output[i] = std::complex(input[i], in[i] * scale); + } + + // Clean up + fftw_destroy_plan(plan_forward); + fftw_destroy_plan(plan_backward); + fftw_free(in); + fftw_free(out); +} + +void WattersonChannel::process(const std::vector& inputSignal, std::vector& outputSignal) +{ + // Apply Hilbert transform to input signal to get complex x_i + std::vector> x_i; + hilbertTransform(inputSignal, x_i); + + generateNoise(x_i); + + // Process the signal through the channel + std::vector> y_i(numSamples); + + // For each sample, compute y_i = h_j[0][i] * x_i + h_j[1][i] * x_{i - (L - 1)} + n_i[i] + for (int i = 0; i < numSamples; ++i) { + std::complex y = n_i[i]; + + for (int pathIndex = 0; pathIndex < numPaths; pathIndex++) { + int delay = delays[pathIndex]; + int idx = i - delay; + if (idx >= 0) { + y += h_j[pathIndex][i] * x_i[idx]; + } + } + + y_i[i] = y; + } + + // Output the double part of y_i + outputSignal.resize(numSamples); + for (int i = 0; i < numSamples; ++i) { + outputSignal[i] = y_i[i].real(); + } +} + +#endif \ No newline at end of file diff --git a/main.cpp b/main.cpp index d579b83..bb97c20 100644 --- a/main.cpp +++ b/main.cpp @@ -1,62 +1,242 @@ -#include -#include +// main.cpp + #include -#include -#include -#include #include +#include +#include +#include +#include // For WAV file handling +// GNU Radio headers +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Include your ModemController and BitStream classes #include "ModemController.h" -#include "PSKModulator.h" +#include "bitstream.h" -BitStream generateBernoulliData(size_t length, double p = 0.5) { +// Function to generate Bernoulli data +BitStream generateBernoulliData(const size_t length, const double p = 0.5, const unsigned int seed = 0) { BitStream random_data; - std::random_device rd; - std::mt19937 gen(rd()); + std::mt19937 gen(seed); std::bernoulli_distribution dist(p); for (size_t i = 0; i < length * 8; ++i) { - random_data.putBit(dist(gen)); // Generates 0 or 1 with probability p + random_data.putBit(dist(gen)); } return random_data; } -int main() { - // Convert sample data to a BitStream object - BitStream input_data = generateBernoulliData(28800); - - // Configuration for modem - size_t baud_rate = 75; - bool is_voice = false; // False indicates data mode - bool is_frequency_hopping = false; // Fixed frequency operation - size_t interleave_setting = 2; // Short interleave - - // Create ModemController instance - ModemController modem(baud_rate, is_voice, is_frequency_hopping, interleave_setting); - - const char* file_name = "modulated_signal_75bps_longinterleave.wav"; - - // Perform transmit operation to generate modulated signal - std::vector modulated_signal = modem.transmit(input_data); - - // Output modulated signal to a WAV file using libsndfile +// Function to write int16_t data to a WAV file +void writeWavFile(const std::string& filename, const std::vector& data, float sample_rate) { SF_INFO sfinfo; sfinfo.channels = 1; - sfinfo.samplerate = 48000; + sfinfo.samplerate = static_cast(sample_rate); sfinfo.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16; - SNDFILE* sndfile = sf_open(file_name, SFM_WRITE, &sfinfo); - if (sndfile == nullptr) { - std::cerr << "Unable to open WAV file for writing modulated signal: " << sf_strerror(sndfile) << "\n"; - return 1; + SNDFILE* outfile = sf_open(filename.c_str(), SFM_WRITE, &sfinfo); + if (!outfile) { + std::cerr << "Error opening output file: " << sf_strerror(nullptr) << std::endl; + return; } - sf_write_short(sndfile, modulated_signal.data(), modulated_signal.size()); - sf_close(sndfile); - std::cout << "Modulated signal written to " << file_name << '\n'; + sf_count_t frames_written = sf_write_short(outfile, data.data(), data.size()); + if (frames_written != static_cast(data.size())) { + std::cerr << "Error writing to output file: " << sf_strerror(outfile) << std::endl; + } - // Success message - std::cout << "Modem test completed successfully.\n"; + sf_close(outfile); +} + +int main() { + // Step 1: Gather parameters and variables + + // Define the preset based on your table (e.g., 4800 bps, 2 fading paths) + struct ChannelPreset { + size_t user_bit_rate; + int num_paths; + bool is_fading; + float multipath_ms; + float fading_bw_hz; + float snr_db; + double target_ber; + }; + + // For this example, let's use the second preset from your table + ChannelPreset preset = { + 4800, // user_bit_rate + 2, // num_paths + true, // is_fading + 2.0f, // multipath_ms + 0.5f, // fading_bw_hz + 27.0f, // snr_db + 1e-3 // target_ber + }; + + // Sampling rate (Hz) + double Fs = 48000.0; // Adjust to match your modem's sampling rate + double Ts = 1.0 / Fs; + + // Carrier frequency (Hz) + float carrier_freq = 1800.0f; // Adjust to match your modem's carrier frequency + + // Step 2: Initialize the modem + size_t baud_rate = preset.user_bit_rate; + bool is_voice = false; + bool is_frequency_hopping = false; + size_t interleave_setting = 2; // Adjust as necessary + + ModemController modem(baud_rate, is_voice, is_frequency_hopping, interleave_setting); + + // Step 3: Generate input modulator data + size_t data_length = 28800; // Length in bytes + unsigned int data_seed = 42; // Random seed + BitStream input_data = generateBernoulliData(data_length, 0.5, data_seed); + + // Step 4: Use the modem to modulate the input data + std::vector passband_signal = modem.transmit(input_data); + + // Write the raw passband audio to a WAV file + writeWavFile("modem_output_raw.wav", passband_signal, Fs); + + // Step 5: Process the modem output through the channel model + + // Convert passband audio to float and normalize + std::vector passband_signal_float(passband_signal.size()); + for (size_t i = 0; i < passband_signal.size(); ++i) { + passband_signal_float[i] = passband_signal[i] / 32768.0f; + } + + // Create GNU Radio top block + auto tb = gr::make_top_block("Passband to Baseband and Channel Model"); + + // Create vector source from passband signal + auto src = gr::blocks::vector_source_f::make(passband_signal_float, false); + + // Apply Hilbert Transform to get analytic signal + int hilbert_taps = 129; // Number of taps + auto hilbert = gr::filter::hilbert_fc::make(hilbert_taps); + + // Multiply by complex exponential to shift to baseband + auto freq_shift_down = gr::analog::sig_source_c::make( + Fs, gr::analog::GR_COS_WAVE, -carrier_freq, 1.0f, 0.0f); + + auto multiplier_down = gr::blocks::multiply_cc::make(); + + // Connect the blocks for downconversion + tb->connect(src, 0, hilbert, 0); + tb->connect(hilbert, 0, multiplier_down, 0); + tb->connect(freq_shift_down, 0, multiplier_down, 1); + + // At this point, multiplier_down outputs the complex baseband signal + + // Configure the channel model parameters + std::vector delays = {0.0f}; + std::vector mags = {1.0f}; + + if (preset.num_paths == 2 && preset.multipath_ms > 0.0f) { + delays.push_back(preset.multipath_ms / 1000.0f); // Convert ms to seconds + float path_gain = 1.0f / sqrtf(2.0f); // Equal average power + mags[0] = path_gain; + mags.push_back(path_gain); + } + + int N = 8; // Number of sinusoids + bool LOS = false; // Rayleigh fading + float K = 0.0f; // K-factor + unsigned int seed = 0; + int ntaps = 64; // Number of taps + + float fD = preset.fading_bw_hz; // Maximum Doppler frequency in Hz + float fDTs = fD * Ts; // Normalized Doppler frequency + + auto channel_model = gr::channels::selective_fading_model::make( + N, fDTs, LOS, K, seed, delays, mags, ntaps); + + // Add AWGN to the signal + float SNR_dB = preset.snr_db; + float SNR_linear = powf(10.0f, SNR_dB / 10.0f); + float signal_power = 0.0f; // Assume normalized + for (const auto& sample : passband_signal_float) { + signal_power += sample * sample; + } + signal_power /= passband_signal_float.size(); + float noise_power = signal_power / SNR_linear; + float noise_voltage = sqrtf(noise_power); + + auto noise_src = gr::analog::noise_source_c::make( + gr::analog::GR_GAUSSIAN, noise_voltage, seed); + + auto adder = gr::blocks::add_cc::make(); + + // Connect the blocks for channel model and noise addition + tb->connect(multiplier_down, 0, channel_model, 0); + tb->connect(channel_model, 0, adder, 0); + tb->connect(noise_src, 0, adder, 1); + + // Multiply by complex exponential to shift back to passband + auto freq_shift_up = gr::analog::sig_source_c::make( + Fs, gr::analog::GR_COS_WAVE, carrier_freq, 1.0f, 0.0f); + + auto multiplier_up = gr::blocks::multiply_cc::make(); + + // Connect the blocks for upconversion + tb->connect(adder, 0, multiplier_up, 0); + tb->connect(freq_shift_up, 0, multiplier_up, 1); + + // Convert to real signal + auto complex_to_real = gr::blocks::complex_to_real::make(); + + // Connect the blocks + tb->connect(multiplier_up, 0, complex_to_real, 0); + + // Collect the output samples + auto sink = gr::blocks::vector_sink_f::make(); + tb->connect(complex_to_real, 0, sink, 0); + + // Run the flowgraph + tb->run(); + + // Retrieve the output data + std::vector received_passband_audio = sink->data(); + + // Normalize and convert to int16_t + // Find maximum absolute value + float max_abs_value = 0.0f; + for (const auto& sample : received_passband_audio) { + if (fabs(sample) > max_abs_value) { + max_abs_value = fabs(sample); + } + } + if (max_abs_value == 0.0f) { + max_abs_value = 1.0f; + } + float scaling_factor = 0.9f / max_abs_value; // Prevent clipping at extremes + + // Apply scaling and convert to int16_t + std::vector received_passband_signal(received_passband_audio.size()); + for (size_t i = 0; i < received_passband_audio.size(); ++i) { + float sample = received_passband_audio[i] * scaling_factor; + // Ensure the sample is within [-1.0, 1.0] + if (sample > 1.0f) sample = 1.0f; + if (sample < -1.0f) sample = -1.0f; + received_passband_signal[i] = static_cast(sample * 32767.0f); + } + + // Step 6: Write the received passband audio to another WAV file + writeWavFile("modem_output_received.wav", received_passband_signal, Fs); + + std::cout << "Processing complete. Output files generated." << std::endl; return 0; -} \ No newline at end of file +} diff --git a/tests/PSKModulatorTests.cpp b/tests/PSKModulatorTests.cpp index cfce492..c77538b 100644 --- a/tests/PSKModulatorTests.cpp +++ b/tests/PSKModulatorTests.cpp @@ -34,7 +34,10 @@ TEST_F(PSKModulatorTest, DemodulationOutput) { } std::cout << std::endl; - auto decoded_symbols = modulator.demodulate(passband_signal); + size_t baud_rate; + size_t interleave_setting; + bool is_voice; + auto decoded_symbols = modulator.demodulate(passband_signal, baud_rate, interleave_setting, is_voice); // Debug: Print decoded symbols std::cout << "Decoded Symbols: ";