Toying with channel models; refactoring symbol former; moved filters back to time domain

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)
@@ -11,20 +12,49 @@ include_directories(include/modulation)
 include_directories(include/utils)
 
 # Add subdirectories for organization
+enable_testing()
 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)
+# 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)

cmake/FindFFTW3.cmake

@@ -0,0 +1,60 @@
+# FindFFTW3.cmake
+# This file is used by CMake to locate the FFTW3 library on the system.
+# It sets the FFTW3_INCLUDE_DIRS and FFTW3_LIBRARIES variables.
+
+# Find the include directory for FFTW3
+find_path(FFTW3_INCLUDE_DIR fftw3.h
+  HINTS
+    ${FFTW3_DIR}/include
+    /usr/include
+    /usr/local/include
+    /opt/local/include
+)
+
+# Find the library for FFTW3
+find_library(FFTW3_LIBRARY fftw3
+  HINTS
+    ${FFTW3_DIR}/lib
+    /usr/lib
+    /usr/local/lib
+    /opt/local/lib
+)
+
+# Find the multi-threaded FFTW3 library, if available
+find_library(FFTW3_THREADS_LIBRARY fftw3_threads
+  HINTS
+    ${FFTW3_DIR}/lib
+    /usr/lib
+    /usr/local/lib
+    /opt/local/lib
+)
+
+# Check if the FFTW3 library was found
+if(FFTW3_INCLUDE_DIR AND FFTW3_LIBRARY)
+  set(FFTW3_FOUND TRUE)
+
+  # Create the FFTW3 imported target
+  add_library(FFTW3::fftw3 UNKNOWN IMPORTED)
+  set_target_properties(FFTW3::fftw3 PROPERTIES
+    IMPORTED_LOCATION ${FFTW3_LIBRARY}
+    INTERFACE_INCLUDE_DIRECTORIES ${FFTW3_INCLUDE_DIR}
+  )
+
+  # Create the FFTW3 Threads imported target, if found
+  if(FFTW3_THREADS_LIBRARY)
+    add_library(FFTW3::fftw3_threads UNKNOWN IMPORTED)
+    set_target_properties(FFTW3::fftw3_threads PROPERTIES
+      IMPORTED_LOCATION ${FFTW3_THREADS_LIBRARY}
+      INTERFACE_INCLUDE_DIRECTORIES ${FFTW3_INCLUDE_DIR}
+    )
+  endif()
+
+  message(STATUS "Found FFTW3: ${FFTW3_LIBRARY}")
+
+else()
+  set(FFTW3_FOUND FALSE)
+  message(STATUS "FFTW3 not found.")
+endif()
+
+# Mark variables as advanced to hide from the cache
+mark_as_advanced(FFTW3_INCLUDE_DIR FFTW3_LIBRARY FFTW3_THREADS_LIBRARY)

include/encoder/Interleaver.h

@@ -180,7 +180,10 @@ private:
      * @brief Sets the matrix dimensions based on baud rate and interleave setting.
      */
     void setMatrixDimensions() {
-        if (baud_rate == 2400) {
+        if (baud_rate == 4800) {
+            rows = 0;
+            columns = 0;
+        } else if (baud_rate == 2400) {
             rows = 40;
             columns = (interleave_setting == 2) ? 576 : 72;
         } else if (baud_rate == 1200) {

include/encoder/ModemController.h

@@ -40,25 +40,24 @@ public:
      * @param data The input data stream to be transmitted. The `is_voice` parameter controls whether the modem treats it as binary file data,
      *             or a binary stream from the MELPe (or other) voice codec.
      */
-    ModemController(const size_t _baud_rate, const bool _is_voice, const bool _is_frequency_hopping, const size_t _interleave_setting, BitStream _data)
+    ModemController(const size_t _baud_rate, const bool _is_voice, const bool _is_frequency_hopping, const size_t _interleave_setting)
         : baud_rate(_baud_rate),
           is_voice(_is_voice),
           is_frequency_hopping(_is_frequency_hopping),
           interleave_setting(_interleave_setting),
-          symbol_formation(baud_rate, interleave_setting, is_voice, is_frequency_hopping),
+          symbol_formation(_baud_rate, _interleave_setting, _is_voice, _is_frequency_hopping),
           scrambler(),
-          fec_encoder(baud_rate, is_frequency_hopping),
-          interleaver(baud_rate, interleave_setting, is_frequency_hopping),
-          input_data(std::move(_data)),
-          mgd_decoder(baud_rate, is_frequency_hopping),
-          modulator(baud_rate, 48000, 0.5, is_frequency_hopping) {}
+          fec_encoder(_baud_rate, _is_frequency_hopping),
+          interleaver(_baud_rate, _interleave_setting, _is_frequency_hopping),
+          mgd_decoder(_baud_rate, _is_frequency_hopping),
+          modulator(48000, _is_frequency_hopping, 48) {}
 
     /**
      * @brief Transmits the input data by processing it through different phases like FEC encoding, interleaving, symbol formation, scrambling, and modulation.
      * @return The scrambled data ready for modulation.
      * @note The modulated signal is generated internally but is intended to be handled externally.
      */
-    std::vector<int16_t> transmit() {
+    std::vector<int16_t> transmit(const BitStream& input_data) {
         // Step 1: Append EOM Symbols
         BitStream eom_appended_data = appendEOMSymbols(input_data);
 
@@ -72,22 +71,29 @@ public:
             // Step 3: Interleaving
             processed_data = interleaver.interleaveStream(fec_encoded_data);
         }
-
+        // Step 4: MGD Decoding
         std::vector<uint8_t> mgd_decoded_data = mgd_decoder.mgdDecode(processed_data);
 
-        // Step 4: Symbol Formation. This function injects the sync preamble symbols. Scrambling is built-in.
+        // Step 5: Symbol Formation. This function injects the sync preamble symbols. Scrambling is handled internally.
         std::vector<uint8_t> symbol_stream = symbol_formation.formSymbols(mgd_decoded_data);
 
+        // Step 6. Modulation. The symbols are applied via 2400-bps 8-PSK modulation, with a 48 KHz sample rate.
         std::vector<int16_t> modulated_signal = modulator.modulate(symbol_stream);
 
         return modulated_signal;
     }
 
+    BitStream receive(const std::vector<int16_t>& passband_signal) {
+        // Step one: Demodulate the passband signal and retrieve decoded symbols
+        std::vector<uint8_t> demodulated_symbols = modulator.demodulate(passband_signal, baud_rate, interleave_setting, is_voice);
+
+        return BitStream();
+    }
+
 private:
     size_t baud_rate;                ///< The baud rate for the modem.
     bool is_voice;                   ///< Indicates if the data being transmitted is voice.
     bool is_frequency_hopping;       ///< Indicates if frequency hopping is used.
-    BitStream input_data;            ///< The input data stream.
     size_t interleave_setting;       ///< The interleave setting to be used.
     size_t sample_rate;
 
@@ -116,9 +122,10 @@ private:
         size_t fec_flush_bits = 144; // FEC encoder flush bits
         size_t interleave_flush_bits = interleaver.getFlushBits();
         size_t total_flush_bits = fec_flush_bits + ((interleave_setting == 0) ? 0 : interleave_flush_bits);
-        while ((eom_data.getMaxBitIndex() + total_flush_bits) % interleave_flush_bits)
-            total_flush_bits++;
-        
+        if (interleave_flush_bits > 0) {
+            while ((eom_data.getMaxBitIndex() + total_flush_bits) % interleave_flush_bits)
+                total_flush_bits++;
+        }
         size_t total_bytes = (total_flush_bits + 7) / 8; // Round up to ensure we have enough bytes to handle all bits.
         BitStream flush_bits(std::vector<uint8_t>(total_bytes, 0), total_flush_bits);
         eom_data += flush_bits;

include/encoder/Scrambler.h

@@ -15,14 +15,14 @@ public:
     /**
      * @brief Constructor initializes the scrambler with a predefined register value.
      */
-    Scrambler() : data_sequence_register(0x0BAD), symbol_count(0) {}
+    Scrambler() : data_sequence_register(0x0BAD), symbol_count(0), preamble_table_index(0) {}
 
     /**
      * @brief Scrambles a synchronization preamble using a fixed randomizer sequence.
      * @param preamble The synchronization preamble to scramble.
      * @return The scrambled synchronization preamble.
      */
-    std::vector<uint8_t> scrambleSyncPreamble(const std::vector<uint8_t>& preamble) const {
+    std::vector<uint8_t> scrambleSyncPreamble(const std::vector<uint8_t>& preamble) {
         static const std::array<uint8_t, 32> sync_randomizer_sequence = {
             7, 4, 3, 0, 5, 1, 5, 0, 2, 2, 1, 1,
             5, 7, 4, 3, 5, 0, 2, 6, 2, 1, 6, 2,
@@ -33,8 +33,9 @@ public:
         scrambled_preamble.reserve(preamble.size());  // Preallocate to improve efficiency
 
         for (size_t i = 0; i < preamble.size(); ++i) {
-            uint8_t scrambled_value = (preamble[i] + sync_randomizer_sequence[i % sync_randomizer_sequence.size()]) % 8;
+            uint8_t scrambled_value = (preamble[i] + sync_randomizer_sequence[preamble_table_index]) % 8;
             scrambled_preamble.push_back(scrambled_value);
+            preamble_table_index = (preamble_table_index + 1) % sync_randomizer_sequence.size();
         }
 
         return scrambled_preamble;
@@ -61,6 +62,7 @@ public:
 private:
     uint16_t data_sequence_register;
     size_t symbol_count;
+    size_t preamble_table_index;
 
     /**
      * @brief Generates the next value from the data sequence randomizing generator.

include/encoder/SymbolFormation.h

@@ -21,16 +21,10 @@ std::vector<uint8_t> baud75_normal_3 = {0, 4, 4, 0};
 
 class SymbolFormation {
     public:
-    SymbolFormation(size_t baud_rate, size_t interleave_setting, bool is_voice, bool is_frequency_hopping) : interleave_setting(interleave_setting), baud_rate(baud_rate), is_voice(is_voice), is_frequency_hopping(is_frequency_hopping) {}
-
-    std::vector<uint8_t> formSymbols(std::vector<uint8_t>& symbol_data) {
-        // Generate and scramble the sync preamble
-        std::vector<uint8_t> sync_preamble = generateSyncPreamble();
-        sync_preamble = scrambler.scrambleSyncPreamble(sync_preamble);
-
+    SymbolFormation(size_t baud_rate, size_t interleave_setting, bool is_voice, bool is_frequency_hopping) : interleave_setting(interleave_setting), baud_rate(baud_rate), is_voice(is_voice), is_frequency_hopping(is_frequency_hopping) {
         // Determine the block sizes
-        size_t unknown_data_block_size = (baud_rate >= 2400) ? 32 : 20;
-        size_t interleaver_block_size;
+        unknown_data_block_size = (baud_rate >= 2400) ? 32 : 20;
+        known_data_block_size = (baud_rate >= 2400) ? 16 : 20;
 
         if (baud_rate == 2400) {
             interleaver_block_size = (interleave_setting == 2) ? (40 * 576) : (40 * 72);
@@ -42,38 +36,43 @@ class SymbolFormation {
             interleaver_block_size = (interleave_setting == 2) ? (20 * 36) : (10 * 9);
         }
 
-        size_t set_count = 0;
-        size_t symbol_count = 0;
+        total_frames = interleaver_block_size / (unknown_data_block_size + known_data_block_size);
+    }
+
+    std::vector<uint8_t> formSymbols(std::vector<uint8_t>& symbol_data) {
+        // Generate and scramble the sync preamble
+        std::vector<uint8_t> sync_preamble = generateSyncPreamble();
+        sync_preamble = scrambler.scrambleSyncPreamble(sync_preamble);
+
         std::vector<uint8_t> 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<uint8_t> 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<uint8_t> 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<uint8_t> 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<uint8_t> 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) {
-                bool is_at_boundary = (symbol_count % interleaver_block_size) == 0;
-                std::vector<uint8_t> probe_data = generateProbeData(!is_at_boundary);
+                // Insert probe data if we are at an interleaver block boundary
+                std::vector<uint8_t> 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;
             }
         }
 
@@ -93,6 +92,10 @@ class SymbolFormation {
     int interleave_setting;
     bool is_voice;
     bool is_frequency_hopping;
+    size_t interleaver_block_size;
+    size_t unknown_data_block_size;
+    size_t known_data_block_size;
+    size_t total_frames;
     Scrambler scrambler = Scrambler();
 
     std::vector<uint8_t> mapChannelSymbolToTribitPattern(uint8_t symbol, bool repeat_twice = false) {
@@ -127,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<uint8_t> 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<uint8_t> generateSyncPreamble() {
@@ -212,59 +212,60 @@ class SymbolFormation {
         return preamble;
     }
 
-    std::vector<uint8_t> generateProbeData(bool is_inside_block) {
+    std::vector<uint8_t> generateProbeData(size_t current_frame, size_t total_frames) {
         std::vector<uint8_t> probe_data;
 
-        // Determine interleaver block size based on baud rate and interleave setting
-        size_t interleaver_block_size;
-        if (baud_rate == 2400) {
-            interleaver_block_size = (interleave_setting == 2) ? (40 * 576) : (40 * 72);
+        // Set the known symbol patterns for D1 and D2 based on Table XI
+        uint8_t D1, D2;
+        if (baud_rate == 4800) {
+            D1 = 7; D2 = 6;
+        } else if (baud_rate == 2400 && is_voice) {
+            D1 = 7; D2 = 7;
+        } else if (baud_rate == 2400) {
+            D1 = (interleave_setting <= 1) ? 6 : 4;
+            D2 = 4;
         } else if (baud_rate == 1200) {
-            interleaver_block_size = (interleave_setting == 2) ? (40 * 288) : (40 * 36);
-        } else if ((baud_rate >= 150) || (baud_rate == 75 && is_frequency_hopping)) {
-            interleaver_block_size = (interleave_setting == 2) ? (40 * 144) : (40 * 18);
+            D1 = (interleave_setting <= 1) ? 6 : 4;
+            D2 = 5;
+        } else if (baud_rate == 600) {
+            D1 = (interleave_setting <= 1) ? 6 : 4;
+            D2 = 6;
+        } else if (baud_rate == 300) {
+            D1 = (interleave_setting <= 1) ? 6 : 4;
+            D2 = 7;
+        } else if (baud_rate == 150) {
+            D1 = (interleave_setting <= 1) ? 7 : 5;
+            D2 = 4;
+        } else if (baud_rate == 75) {
+            D1 = (interleave_setting <= 1) ? 7 : 5;
+            D2 = 5;
         } else {
-            interleaver_block_size = (interleave_setting == 2) ? (20 * 36) : (10 * 9);
+            throw std::invalid_argument("Invalid baud rate for generateProbeData");
         }
 
-        // If we are inside an interleaver block, the probe data is filled with zeros
-        if (is_inside_block) {
-            probe_data.resize(interleaver_block_size, 0x00);
-        } else {
-            // Set the known symbol patterns for D1 and D2 based on Table XI
-            uint8_t D1, D2;
-            if (baud_rate == 4800) {
-                D1 = 7; D2 = 6;
-            } else if (baud_rate == 2400 && is_voice) {
-                D1 = 7; D2 = 7;
-            } else if (baud_rate == 2400) {
-                D1 = (interleave_setting <= 1) ? 6 : 4;
-                D2 = 4;
-            } else if (baud_rate == 1200) {
-                D1 = (interleave_setting <= 1) ? 6 : 4;
-                D2 = 5;
-            } else if (baud_rate == 600) {
-                D1 = (interleave_setting <= 1) ? 6 : 4;
-                D2 = 6;
-            } else if (baud_rate == 300) {
-                D1 = (interleave_setting <= 1) ? 6 : 4;
-                D2 = 7;
-            } else if (baud_rate == 150) {
-                D1 = (interleave_setting <= 1) ? 7 : 5;
-                D2 = 4;
-            } else if (baud_rate == 75) {
-                D1 = (interleave_setting <= 1) ? 7 : 5;
-                D2 = 5;
-            } else {
-                throw std::invalid_argument("Invalid baud rate for generateProbeData");
-            }
-
-            // Generate the known symbol patterns D1 and D2, repeated twice
+        // If the current frame is not the last two frames, set probe data to zeros
+        if (current_frame < total_frames - 2) {
+            probe_data.resize(known_data_block_size, 0x00);
+        }
+        // If the current frame is the second-to-last frame, set probe data to D1 pattern
+        else if (current_frame == total_frames - 2) {
             std::vector<uint8_t> d1_pattern = mapChannelSymbolToTribitPattern(D1, true);
-            std::vector<uint8_t> d2_pattern = mapChannelSymbolToTribitPattern(D2, true);
-
             probe_data.insert(probe_data.end(), d1_pattern.begin(), d1_pattern.end());
+
+            // Fill the remaining symbols with zeros if necessary
+            if (probe_data.size() < known_data_block_size) {
+                probe_data.resize(known_data_block_size, 0x00);
+            }
+        }
+        // If the current frame is the last frame, set probe data to D2 pattern
+        else if (current_frame == total_frames - 1) {
+            std::vector<uint8_t> d2_pattern = mapChannelSymbolToTribitPattern(D2, true);
             probe_data.insert(probe_data.end(), d2_pattern.begin(), d2_pattern.end());
+
+            // Fill the remaining symbols with zeros if necessary
+            if (probe_data.size() < known_data_block_size) {
+                probe_data.resize(known_data_block_size, 0x00);
+            }
         }
 
         return probe_data;
@@ -332,19 +333,6 @@ class SymbolFormation {
         }
     }
 
-    std::vector<uint8_t> map75bpsSet(const std::vector<uint8_t>& data, size_t start_index, size_t set_size, bool is_exceptional_set) {
-        std::vector<uint8_t> 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<uint8_t> mapUnknownData(const std::vector<uint8_t>& data) {
         std::vector<uint8_t> mapped_data;
 

include/modulation/FSKDemodulator.h

@@ -1,131 +0,0 @@
-#ifndef FSK_DEMODULATOR_H
-#define FSK_DEMODULATOR_H
-
-#include <cmath>
-#include <cstdint>
-#include <functional>
-#include <memory>
-#include <vector>
-
-#include "BitStreamWriter.h"
-
-class FSKDemodulatorConfig {
-    public:
-    int freq_lo;
-    int freq_hi;
-    int sample_rate;
-    int baud_rate;
-    std::shared_ptr<BitStreamWriter> bitstreamwriter;
-};
-
-namespace milstd {
-    class FSKDemodulator {
-        public:
-        FSKDemodulator(const FSKDemodulatorConfig& s) : freq_lo(s.freq_lo), freq_hi(s.freq_hi), sample_rate(s.sample_rate), baud_rate(s.baud_rate), bit_writer(s.bitstreamwriter) {
-            initialize();
-        }
-
-        void demodulate(const std::vector<int16_t>& samples) {
-            size_t nb = samples.size();
-
-            for (size_t i = 0; i < nb; i++) {
-                filter_buf[buf_ptr++] = samples[i];
-
-                if (buf_ptr == filter_buf.size()) {
-                    std::copy(filter_buf.begin() + filter_buf.size() - filter_size, filter_buf.end(), filter_buf.begin());
-                    buf_ptr = filter_size;
-                }
-
-                int corr;
-                int sum = 0;
-
-                corr = dotProduct(&filter_buf[buf_ptr - filter_size], filter_hi_i.data(), filter_size);
-                sum += corr * corr;
-
-                corr = dotProduct(&filter_buf[buf_ptr - filter_size], filter_hi_q.data(), filter_size);
-                sum += corr * corr;
-
-                corr = dotProduct(&filter_buf[buf_ptr - filter_size], filter_lo_i.data(), filter_size);
-                sum -= corr * corr;
-
-                corr = dotProduct(&filter_buf[buf_ptr - filter_size], filter_lo_q.data(), filter_size);
-                sum -= corr * corr;
-
-                int new_sample = (sum > 0) ? 1 : 0;
-
-                if (last_sample != new_sample) {
-                    last_sample = new_sample;
-                    if (baud_pll < 0.5)
-                        baud_pll += baud_pll_adj;
-                    else
-                        baud_pll -= baud_pll_adj;
-                }
-
-                baud_pll += baud_incr;
-
-                if (baud_pll >= 1.0) {
-                    baud_pll -= 1.0;
-                    bit_writer->putBit(last_sample);
-                }
-            }
-        }
-
-        private:
-        int freq_lo;
-        int freq_hi;
-        int sample_rate;
-        int baud_rate;
-        std::shared_ptr<BitStreamWriter> bit_writer;
-
-        int filter_size;
-        std::vector<double> filter_lo_i;
-        std::vector<double> filter_lo_q;
-        std::vector<double> filter_hi_i;
-        std::vector<double> filter_hi_q;
-        std::vector<double> filter_buf;
-        size_t buf_ptr;
-
-        double baud_incr;
-        double baud_pll;
-        double baud_pll_adj;
-        int last_sample;
-
-        void initialize() {
-            baud_incr = static_cast<double>(baud_rate) / sample_rate;
-            baud_pll = 0.0;
-            baud_pll_adj = baud_incr / 4;
-
-            filter_size = sample_rate / baud_rate;
-
-            filter_buf.resize(filter_size * 2, 0.0);
-            buf_ptr = filter_size;
-
-            last_sample = 0;
-
-            filter_lo_i.resize(filter_size);
-            filter_lo_q.resize(filter_size);
-            filter_hi_i.resize(filter_size);
-            filter_hi_q.resize(filter_size);
-            
-            for (int i = 0; i < filter_size; i++) {
-                double phase_lo = 2.0 * M_PI * freq_lo * i / sample_rate;
-                filter_lo_i[i] = std::cos(phase_lo);
-                filter_lo_q[i] = std::sin(phase_lo);
-                
-                double phase_hi = 2.0 * M_PI * freq_hi * i / sample_rate;
-                filter_hi_i[i] = std::cos(phase_hi);
-                filter_hi_q[i] = std::sin(phase_hi);
-            }
-        }
-
-        double dotProduct(const double* x, const double* y, size_t size) {
-            double sum = 0.0;
-            for (size_t i = 0; i < size; i++) {
-                sum += x[i] * y[i];
-            }
-            return sum;
-        }
-    };
-} // namespace milstd
-
-#endif /* FSK_DEMODULATOR_H */

include/modulation/FSKModulator.h

@@ -1,87 +0,0 @@
-#ifndef FSK_MODULATOR_H
-#define FSK_MODULATOR_H
-
-#include <cmath>
-#include <cstdint>
-#include <functional>
-#include <memory>
-#include <vector>
-
-#include "BitStreamReader.h"
-
-class FSKModulatorConfig {
-    public:
-    int freq_lo;
-    int freq_hi;
-    int sample_rate;
-    int baud_rate;
-    std::shared_ptr<BitStreamReader> bitstreamreader;
-};
-
-namespace milstd {
-    class FSKModulator {
-        public:
-        FSKModulator(const FSKModulatorConfig& s) : freq_lo(s.freq_lo), freq_hi(s.freq_hi), sample_rate(s.sample_rate), baud_rate(s.baud_rate), bit_reader(s.bitstreamreader) {
-            omega[0] = (2.0 * M_PI * freq_lo) / sample_rate;
-            omega[1] = (2.0 * M_PI * freq_hi) / sample_rate;
-            baud_incr = static_cast<double>(baud_rate) / sample_rate;
-            phase = 0.0;
-            baud_frac = 0.0;
-            current_bit = 0;
-        }
-
-        std::vector<int16_t> modulate(unsigned int num_samples) {
-            std::vector<int16_t> samples;
-            samples.reserve(num_samples);
-
-            int bit = current_bit;
-
-            for (unsigned int i = 0; i < num_samples; i++) {
-                baud_frac += baud_incr;
-                if (baud_frac >= 1.0) {
-                    baud_frac -= 1.0;
-                    try
-                    {
-                        bit = bit_reader->getNextBit();
-                    }
-                    catch(const std::out_of_range&)
-                    {
-                        bit = 0;
-                    }
-                    if (bit != 0 && bit != 1)
-                        bit = 0;
-                }
-
-                double sample = std::cos(phase);
-                int16_t sample_int = static_cast<int16_t>(sample * 32767);
-                samples.push_back(sample_int);
-
-                phase += omega[bit];
-                if (phase >= 2.0 * M_PI) {
-                    phase -= 2.0 * M_PI;
-                }
-            }
-
-            current_bit = bit;
-
-            return samples;
-        }
-
-        private:
-        // parameters
-        int freq_lo, freq_hi;
-        int sample_rate;
-        int baud_rate;
-        std::shared_ptr<BitStreamReader> bit_reader;
-
-        // state variables
-        double phase;
-        double baud_frac;
-        double baud_incr;
-        std::array<double, 2> omega;
-        int current_bit;
-    };
-
-} // namespace milstd
-
-#endif /* FSK_MODULATOR_H */

include/modulation/PSKModulator.h

@@ -1,136 +1,366 @@
 #ifndef PSK_MODULATOR_H
 #define PSK_MODULATOR_H
 
-#include <vector>
-#include <cmath>
-#include <cstdint>
-#include <stdexcept>
-#include <complex>
+
 #include <algorithm>
+#include <array>
+#include <cmath>
+#include <complex>
+#include <cstdint>
+#include <numeric>
+#include <stdexcept>
+#include <vector>
+#include <fftw3.h>
+#include <map>
+#include <tuple>
+
+#include "costasloop.h"
+#include "filters.h"
+#include "Scrambler.h"
+
+static constexpr double CARRIER_FREQ = 1800.0;
+static constexpr size_t SYMBOL_RATE = 2400;
+static constexpr double ROLLOFF_FACTOR = 0.35;
+static constexpr double SCALE_FACTOR = 32767.0;
 
 class PSKModulator {
 public:
-    PSKModulator(double baud_rate, double sample_rate, double energy_per_bit, bool is_frequency_hopping) 
-        : sample_rate(sample_rate), carrier_freq(1800), phase(0.0) {
-        initializeSymbolMap();
-        symbol_rate = 2400;  // Fixed symbol rate as per specification (2400 symbols per second)
-        samples_per_symbol = static_cast<size_t>(sample_rate / symbol_rate);
+    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<size_t>(sample_rate / SYMBOL_RATE)), srrc_filter(8, _sample_rate, SYMBOL_RATE, ROLLOFF_FACTOR) {
+            initializeSymbolMap();
+            phase_detector = PhaseDetector(symbolMap);
     }
 
     std::vector<int16_t> modulate(const std::vector<uint8_t>& symbols) {
-        std::vector<std::complex<double>> modulated_signal;
+        std::vector<std::complex<double>> baseband_components(symbols.size() * samples_per_symbol);
+        size_t symbol_index = 0;
 
-        const double phase_increment = 2 * M_PI * carrier_freq / sample_rate;
-        for (auto symbol : symbols) {
+        for (const auto& symbol : symbols) {
             if (symbol >= symbolMap.size()) {
-                throw std::out_of_range("Invalid symbol value for 8-PSK modulation");
+                throw std::out_of_range("Invalid symbol value for 8-PSK modulation. Symbol must be between 0 and 7.");
             }
-            std::complex<double> target_symbol = symbolMap[symbol];
+            const std::complex<double> target_symbol = symbolMap[symbol];
 
             for (size_t i = 0; i < samples_per_symbol; ++i) {
-                double in_phase = std::cos(phase + target_symbol.real());
-                double quadrature = std::sin(phase + target_symbol.imag());
-                modulated_signal.emplace_back(in_phase, quadrature);
-                phase = std::fmod(phase + phase_increment, 2 * M_PI);
+                baseband_components[symbol_index * samples_per_symbol + i] = target_symbol;
             }
+
+            symbol_index++;
         }
 
-        // Apply raised-cosine filter
-        auto filter_taps = sqrtRaisedCosineFilter(201, symbol_rate);  // Adjusted number of filter taps to 201 for balance
-        auto filtered_signal = applyFilter(modulated_signal, filter_taps);
+        // Filter the I/Q phase components
+        std::vector<std::complex<double>> filtered_components = srrc_filter.applyFilter(baseband_components);
 
-        // Normalize the filtered signal
-        double max_value = 0.0;
-        for (const auto& sample : filtered_signal) {
-            max_value = std::max(max_value, std::abs(sample.real()));
-            max_value = std::max(max_value, std::abs(sample.imag()));
-        }
-        double gain = (max_value > 0) ? (32767.0 / max_value) : 1.0;
+        // Combine the I and Q components
+        std::vector<double> passband_signal;
+        passband_signal.reserve(baseband_components.size());
 
-        // Combine the I and Q components and apply gain for audio output
-        std::vector<int16_t> combined_signal;
-        for (auto& sample : filtered_signal) {
-            int16_t combined_sample = static_cast<int16_t>(std::clamp(gain * (sample.real() + sample.imag()), -32768.0, 32767.0));
-            combined_signal.push_back(combined_sample);
+        double carrier_phase = 0.0;
+        double carrier_phase_increment = 2 * M_PI * CARRIER_FREQ / sample_rate;
+        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 * SCALE_FACTOR); // Scale to int16_t
+            carrier_phase += carrier_phase_increment;
+            if (carrier_phase >= 2 * M_PI)
+                carrier_phase -= 2 * M_PI;
         }
 
-        return combined_signal;
+        std::vector<int16_t> final_signal;
+        final_signal.reserve(passband_signal.size());
+
+        for (const auto& sample : passband_signal) {
+            int16_t value = static_cast<int16_t>(sample);
+            value = std::clamp(value, (int16_t)-32768, (int16_t)32767);
+            final_signal.emplace_back(value);
+        }
+
+        return final_signal;
     }
 
-    std::vector<double> sqrtRaisedCosineFilter(size_t num_taps, double symbol_rate) {
-        double rolloff = 0.35;  // Fixed rolloff factor as per specification
-        std::vector<double> filter_taps(num_taps);
-        double norm_factor = 0.0;
-        double sampling_interval = 1.0 / sample_rate;
-        double symbol_duration = 1.0 / symbol_rate;
-        double half_num_taps = static_cast<double>(num_taps - 1) / 2.0;
+    std::vector<uint8_t> demodulate(const std::vector<int16_t> passband_signal, size_t& baud_rate, size_t& interleave_setting, bool& is_voice) {
+        // Carrier recovery. initialize the Costas loop.
+        CostasLoop costas_loop(CARRIER_FREQ, sample_rate, symbolMap, 5.0, 0.05, 0.01);
 
-        for (size_t i = 0; i < num_taps; ++i) {
-            double t = (i - half_num_taps) * sampling_interval;
-            if (std::abs(t) < 1e-10) {
-                filter_taps[i] = 1.0;
-            } else {
-                double numerator = std::sin(M_PI * t / symbol_duration * (1.0 - rolloff)) +
-                                  4.0 * rolloff * t / symbol_duration * std::cos(M_PI * t / symbol_duration * (1.0 + rolloff));
-                double denominator = M_PI * t * (1.0 - std::pow(4.0 * rolloff * t / symbol_duration, 2));
-                filter_taps[i] = numerator / denominator;
-            }
-            norm_factor += filter_taps[i] * filter_taps[i];
+        // Convert passband signal to doubles.
+        std::vector<double> normalized_passband(passband_signal.size());
+        for (size_t i = 0; i < passband_signal.size(); i++) {
+            normalized_passband[i] = passband_signal[i] / 32767.0;
         }
 
-        norm_factor = std::sqrt(norm_factor);
-        std::for_each(filter_taps.begin(), filter_taps.end(), [&norm_factor](double &tap) { tap /= norm_factor; });
-        return filter_taps;
-    }
+        // Downmix passband to baseband
+        std::vector<std::complex<double>> baseband_IQ = costas_loop.process(normalized_passband);
+        std::vector<uint8_t> detected_symbols;
 
-    std::vector<std::complex<double>> applyFilter(const std::vector<std::complex<double>>& signal, const std::vector<double>& filter_taps) {
-        std::vector<std::complex<double>> filtered_signal(signal.size());
+        // Phase detection and symbol formation
+        size_t samples_per_symbol = sample_rate / SYMBOL_RATE;
+        bool sync_found = false;
+        size_t sync_segments_detected;
 
-        size_t filter_length = filter_taps.size();
-        size_t half_filter_length = filter_length / 2;
-
-        // Convolve the signal with the filter taps
-        for (size_t i = 0; i < signal.size(); ++i) {
-            double filtered_i = 0.0;
-            double filtered_q = 0.0;
-
-            for (size_t j = 0; j < filter_length; ++j) {
-                if (i >= j) {
-                    filtered_i += filter_taps[j] * signal[i - j].real();
-                    filtered_q += filter_taps[j] * signal[i - j].imag();
-                } else {
-                    // Handle edge case by zero-padding
-                    filtered_i += filter_taps[j] * 0.0;
-                    filtered_q += filter_taps[j] * 0.0;
-                }
+        size_t window_size = 32*15;
+        
+        for (size_t i = 0; i < baseband_IQ.size(); i += samples_per_symbol) {
+            std::complex<double> symbol_avg(0.0, 0.0);
+            for (size_t j = 0; j < samples_per_symbol; j++) {
+                symbol_avg += baseband_IQ[i + j];
             }
+            symbol_avg /= static_cast<double>(samples_per_symbol);
 
-            filtered_signal[i] = std::complex<double>(filtered_i, filtered_q);
+            uint8_t detected_symbol = phase_detector.getSymbol(symbol_avg);
+            detected_symbols.push_back(detected_symbol);
         }
 
-        return filtered_signal;
+        if (processSyncSegments(detected_symbols, baud_rate, interleave_setting, is_voice)) {
+            return processDataSymbols(detected_symbols);
+        }
+
+        return std::vector<uint8_t>();
     }
 
 private:
-    double sample_rate;        ///< The sample rate of the system.
-    double carrier_freq;       ///< The frequency of the carrier, set to 1800 Hz as per standard.
-    double phase;              ///< Current phase of the carrier waveform.
+    const double sample_rate;        ///< The sample rate of the system.
+    const double gain;               ///< The gain of the modulated signal.
     size_t samples_per_symbol; ///< Number of samples per symbol, calculated to match symbol duration with cycle.
-    size_t symbol_rate;
+    PhaseDetector phase_detector;
+    SRRCFilter srrc_filter;
     std::vector<std::complex<double>> symbolMap;  ///< The mapping of tribit symbols to I/Q components.
+    const bool is_frequency_hopping;                    ///< Whether to use frequency hopping methods. Not implemented (yet?)
+    
 
-    void initializeSymbolMap() {
+    static inline double validateSampleRate(const double rate) {
+        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;
+    }
+
+    inline void initializeSymbolMap() {
         symbolMap = {
-            {1.0, 0.0},                              // 0 (000) corresponds to I = 1.0, Q = 0.0
-            {std::sqrt(2.0) / 2.0, std::sqrt(2.0) / 2.0}, // 1 (001) corresponds to I = cos(45), Q = sin(45)
-            {0.0, 1.0},                              // 2 (010) corresponds to I = 0.0, Q = 1.0
-            {-std::sqrt(2.0) / 2.0, std::sqrt(2.0) / 2.0}, // 3 (011) corresponds to I = cos(135), Q = sin(135)
-            {-1.0, 0.0},                             // 4 (100) corresponds to I = -1.0, Q = 0.0
-            {-std::sqrt(2.0) / 2.0, -std::sqrt(2.0) / 2.0}, // 5 (101) corresponds to I = cos(225), Q = sin(225)
-            {0.0, -1.0},                             // 6 (110) corresponds to I = 0.0, Q = -1.0
-            {std::sqrt(2.0) / 2.0, -std::sqrt(2.0) / 2.0}  // 7 (111) corresponds to I = cos(315), Q = sin(315)
+            {gain * std::cos(2.0*M_PI*(0.0/8.0)), gain * std::sin(2.0*M_PI*(0.0/8.0))}, // 0 (000) corresponds to I = 1.0, Q = 0.0
+            {gain * std::cos(2.0*M_PI*(1.0/8.0)), gain * std::sin(2.0*M_PI*(1.0/8.0))}, // 1 (001) corresponds to I = cos(45), Q = sin(45)
+            {gain * std::cos(2.0*M_PI*(2.0/8.0)), gain * std::sin(2.0*M_PI*(2.0/8.0))}, // 2 (010) corresponds to I = 0.0, Q = 1.0
+            {gain * std::cos(2.0*M_PI*(3.0/8.0)), gain * std::sin(2.0*M_PI*(3.0/8.0))}, // 3 (011) corresponds to I = cos(135), Q = sin(135)
+            {gain * std::cos(2.0*M_PI*(4.0/8.0)), gain * std::sin(2.0*M_PI*(4.0/8.0))}, // 4 (100) corresponds to I = -1.0, Q = 0.0
+            {gain * std::cos(2.0*M_PI*(5.0/8.0)), gain * std::sin(2.0*M_PI*(5.0/8.0))}, // 5 (101) corresponds to I = cos(225), Q = sin(225)
+            {gain * std::cos(2.0*M_PI*(6.0/8.0)), gain * std::sin(2.0*M_PI*(6.0/8.0))}, // 6 (110) corresponds to I = 0.0, Q = -1.0
+            {gain * std::cos(2.0*M_PI*(7.0/8.0)), gain * std::sin(2.0*M_PI*(7.0/8.0))}  // 7 (111) corresponds to I = cos(315), Q = sin(315)
         };
     }
+
+    uint8_t extractBestTribit(const std::vector<uint8_t>& stream, const size_t start, const size_t window_size) const {
+        if (start + window_size > stream.size()) {
+            throw std::out_of_range("Window size exceeds symbol stream size.");
+        }
+
+        Scrambler scrambler;
+        std::vector<uint8_t> symbol(stream.begin() + start, stream.begin() + start + window_size);
+        std::vector<uint8_t> descrambled_symbol = scrambler.scrambleSyncPreamble(symbol);
+
+        const size_t split_len = window_size / 4;
+        std::array<uint8_t, 8> tribit_counts = {0};  // Counts for each channel symbol (000 to 111)
+
+        // Loop through each split segment (4 segments)
+        for (size_t i = 0; i < 4; ++i) {
+            // Extract the range for this split
+            size_t segment_start = start + i * split_len;
+            size_t segment_end = segment_start + split_len;
+
+            // Compare this segment to the predefined patterns from the table and map to a channel symbol
+            uint8_t tribit_value = mapSegmentToChannelSymbol(descrambled_symbol, segment_start, segment_end);
+
+            // Increment the corresponding channel symbol count
+            tribit_counts[tribit_value]++;
+        }
+
+        // Find the channel symbol with the highest count (majority vote)
+        uint8_t best_symbol = std::distance(tribit_counts.begin(), std::max_element(tribit_counts.begin(), tribit_counts.end()));
+
+        return best_symbol;
+    }
+
+    // Function to map a segment of the stream back to a channel symbol based on the repeating patterns
+    uint8_t mapSegmentToChannelSymbol(const std::vector<uint8_t>& segment, size_t start, size_t end) const {
+        std::vector<uint8_t> extracted_pattern(segment.begin() + start, segment.begin() + end);
+
+        // Compare the extracted pattern with known patterns from the table
+        if (matchesPattern(extracted_pattern, {0, 0, 0, 0, 0, 0, 0, 0})) return 0b000;
+        if (matchesPattern(extracted_pattern, {0, 4, 0, 4, 0, 4, 0, 4})) return 0b001;
+        if (matchesPattern(extracted_pattern, {0, 0, 4, 4, 0, 0, 4, 4})) return 0b010;
+        if (matchesPattern(extracted_pattern, {0, 4, 4, 0, 0, 4, 4, 0})) return 0b011;
+        if (matchesPattern(extracted_pattern, {0, 0, 0, 0, 4, 4, 4, 4})) return 0b100;
+        if (matchesPattern(extracted_pattern, {0, 4, 0, 4, 4, 0, 4, 0})) return 0b101;
+        if (matchesPattern(extracted_pattern, {0, 0, 4, 4, 4, 4, 0, 0})) return 0b110;
+        if (matchesPattern(extracted_pattern, {0, 4, 4, 0, 4, 0, 0, 4})) return 0b111;
+
+        throw std::invalid_argument("Invalid segment pattern");
+    }
+
+    // Helper function to compare two patterns
+    bool matchesPattern(const std::vector<uint8_t>& segment, const std::vector<uint8_t>& pattern) const {
+        return std::equal(segment.begin(), segment.end(), pattern.begin());
+    }
+
+    bool configureModem(uint8_t D1, uint8_t D2, size_t& baud_rate, size_t& interleave_setting, bool& is_voice) {
+        // Predefine all the valid combinations in a lookup map
+        static const std::map<std::pair<uint8_t, uint8_t>, std::tuple<size_t, size_t, bool>> modemConfig = {
+            {{7, 6}, {4800, 1, false}},  // 4800 bps
+            {{7, 7}, {2400, 1, true}},   // 2400 bps, voice
+            {{6, 4}, {2400, 1, false}},  // 2400 bps, data
+            {{6, 5}, {1200, 1, false}},  // 1200 bps
+            {{6, 6}, {600,  1, false}},  // 600 bps
+            {{6, 7}, {300,  1, false}},  // 300 bps
+            {{7, 4}, {150,  1, false}},  // 150 bps
+            {{7, 5}, {75,   1, false}},  // 75 bps
+            {{4, 4}, {2400, 2, false}},  // 2400 bps, long interleave
+            {{4, 5}, {1200, 2, false}},  // 1200 bps, long interleave
+            {{4, 6}, {600,  2, false}},  // 600 bps, long interleave
+            {{4, 7}, {300,  2, false}},  // 300 bps, long interleave
+            {{5, 4}, {150,  2, false}},  // 150 bps, long interleave
+            {{5, 5}, {75,   2, false}},  // 75 bps, long interleave
+        };
+
+        // Use D1 and D2 to look up the correct configuration
+        auto it = modemConfig.find({D1, D2});
+        if (it != modemConfig.end()) {
+            // Set the parameters if found
+            std::tie(baud_rate, interleave_setting, is_voice) = it->second;
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+    uint8_t calculateSegmentCount(const uint8_t C1, const uint8_t C2, const uint8_t C3) {
+        uint8_t extracted_C1 = C1 & 0b11;
+        uint8_t extracted_C2 = C2 & 0b11;
+        uint8_t extracted_C3 = C3 & 0b11;
+        
+        uint8_t segment_count = (extracted_C1 << 4) | (extracted_C2 << 2) | extracted_C3;
+        return segment_count;
+    }
+
+    bool processSegment(const std::vector<uint8_t>& detected_symbols, size_t& start, size_t symbol_size, size_t& segment_count, uint8_t& D1, uint8_t& D2) {
+        size_t sync_pattern_length = 9;
+
+        if (start + symbol_size * sync_pattern_length > detected_symbols.size()) {
+            start = detected_symbols.size();
+            return false;
+        }
+
+        std::vector<uint8_t> window(detected_symbols.begin() + start, detected_symbols.begin() + start + sync_pattern_length * symbol_size);
+        std::vector<uint8_t> extracted_window;
+        for (size_t i = 0; i < sync_pattern_length; i++) {
+            extracted_window.push_back(extractBestTribit(window, i * symbol_size, symbol_size));
+        }
+
+        if (!matchesPattern(extracted_window, {0, 1, 3, 0, 1, 3, 1, 2, 0})) {
+            start += symbol_size;
+            return false;
+        }
+
+        start += sync_pattern_length * symbol_size;
+        size_t D1_index = start + symbol_size;
+        size_t D2_index = D1_index + symbol_size;
+
+        if (D2_index + symbol_size > detected_symbols.size()) {
+            start = detected_symbols.size();
+            return false;
+        }
+
+        D1 = extractBestTribit(detected_symbols, D1_index, symbol_size);
+        D2 = extractBestTribit(detected_symbols, D2_index, symbol_size);
+
+        // Process the count symbols (C1, C2, C3)
+        size_t C1_index = D2_index + symbol_size;
+        size_t C2_index = C1_index + symbol_size;
+        size_t C3_index = C2_index + symbol_size;
+
+        if (C3_index + symbol_size > detected_symbols.size()) {
+            start = detected_symbols.size();
+            return false;
+        }
+
+        uint8_t C1 = extractBestTribit(detected_symbols, C1_index, symbol_size);
+        uint8_t C2 = extractBestTribit(detected_symbols, C2_index, symbol_size);
+        uint8_t C3 = extractBestTribit(detected_symbols, C3_index, symbol_size);
+
+        segment_count = calculateSegmentCount(C1, C2, C3);
+
+        // Check for the constant zero pattern
+        size_t constant_zero_index = C3_index + symbol_size;
+
+        if (constant_zero_index + symbol_size > detected_symbols.size()) {
+            start = detected_symbols.size();
+            return false;
+        }
+        uint8_t constant_zero = extractBestTribit(detected_symbols, constant_zero_index, symbol_size);
+
+        if (constant_zero != 0) {
+            start = constant_zero_index + symbol_size;
+            return false; // Failed zero check, resync
+        }
+
+        // Successfully processed the segment
+        start = constant_zero_index + symbol_size; // Move start to next segment
+        return true;
+    }
+
+    bool processSyncSegments(const std::vector<uint8_t>& detected_symbols, size_t& baud_rate, size_t& interleave_setting, bool& is_voice) {
+        size_t symbol_size = 32;
+        size_t start = 0;
+        size_t segment_count = 0;
+        std::map<std::pair<uint8_t, uint8_t>, int> vote_map;
+        const int short_interleave_threshold = 2;
+        const int long_interleave_threshold = 5;
+
+        // Attempt to detect interleave setting dynamically
+        bool interleave_detected = false;
+        int current_threshold = short_interleave_threshold;  // Start by assuming short interleave
+
+        while (start + symbol_size * 15 < detected_symbols.size()) {
+            uint8_t D1 = 0, D2 = 0;
+            if (processSegment(detected_symbols, start, symbol_size, segment_count, D1, D2)) {
+                vote_map[{D1, D2}]++;
+
+                // Check if we have enough votes to make a decision based on current interleave assumption
+                if (vote_map.size() >= current_threshold) {
+                    auto majority_vote = std::max_element(vote_map.begin(), vote_map.end(), [](const auto& a, const auto& b) { return a.second < b.second; });
+
+                    if (configureModem(majority_vote->first.first, majority_vote->first.second, baud_rate, interleave_setting, is_voice)) {
+                        interleave_detected = true;
+                        break;  // Successfully configured modem, exit loop
+                    } else {
+                        // If configuration fails, retry with the other interleave type
+                        if (current_threshold == short_interleave_threshold) {
+                            current_threshold = long_interleave_threshold;  // Switch to long interleave
+                            vote_map.clear();  // Clear the vote map and start fresh
+                            start = 0;  // Restart segment processing
+                        } else {
+                            continue;  // Both short and long interleave attempts failed, signal is not usable
+                        }
+                    }
+                }
+
+                if (segment_count > 0) {
+                    while (segment_count > 0 && start < detected_symbols.size()) {
+                        uint8_t dummy_D1, dummy_D2;
+                        if (!processSegment(detected_symbols, start, symbol_size, segment_count, dummy_D1, dummy_D2)) {
+                            continue;
+                        }
+                    }
+                }
+            } else {
+                start += symbol_size;  // Move to the next segment
+            }
+        }
+
+        return interleave_detected;
+    }
+
+    std::vector<uint8_t> processDataSymbols(const std::vector<uint8_t>& detected_symbols) {
+        return std::vector<uint8_t>();
+    }
 };
 
 #endif

include/utils/bitstream.h

@@ -19,7 +19,7 @@ public:
     /**
      * @brief Default constructor.
      */
-    BitStream() : bit_index(0), max_bit_idx(0) {}
+    BitStream() : std::vector<uint8_t>(), bit_index(0), max_bit_idx(0) {}
 
     /**
      * @brief Constructs a BitStream from an existing vector of bytes.

include/utils/costasloop.h

@@ -0,0 +1,113 @@
+#include <complex>
+#include <cmath>
+#include <vector>
+#include <iostream>
+
+#include "filters.h"
+
+class PhaseDetector {
+public:
+    PhaseDetector() {}
+    PhaseDetector(const std::vector<std::complex<double>>& _symbolMap) : symbolMap(_symbolMap) {}
+
+    uint8_t getSymbol(const std::complex<double>& input) {
+        double phase = std::atan2(input.imag(), input.real());
+        return symbolFromPhase(phase);
+    }
+
+private:
+    std::vector<std::complex<double>> symbolMap;
+
+    uint8_t symbolFromPhase(const double phase) {
+        // Calculate the closest symbol based on phase difference
+        double min_distance = 2 * M_PI; // Maximum possible phase difference
+        uint8_t closest_symbol = 0;
+
+        for (uint8_t i = 0; i < symbolMap.size(); ++i) {
+            double symbol_phase = std::atan2(symbolMap[i].imag(), symbolMap[i].real());
+            double distance = std::abs(symbol_phase - phase);
+
+            if (distance < min_distance) {
+                min_distance = distance;
+                closest_symbol = i;
+            }
+        }
+
+        return closest_symbol;
+    }
+};
+
+class CostasLoop {
+public:
+    CostasLoop(const double _carrier_freq, const double _sample_rate, const std::vector<std::complex<double>>& _symbolMap, const double _vco_gain, const double _alpha, const double _beta) 
+        : carrier_freq(_carrier_freq), sample_rate(_sample_rate), vco_gain(_vco_gain), alpha(_alpha), beta(_beta), freq_error(0.0), k_factor(-1 / (_sample_rate * _vco_gain)),
+          prev_in_iir(0), prev_out_iir(0), prev_in_vco(0), feedback(1.0, 0.0),
+          error_total(0), out_iir_total(0), in_vco_total(0),
+          srrc_filter(SRRCFilter(48, _sample_rate, 2400, 0.35)) {}
+
+    std::vector<std::complex<double>> process(const std::vector<double>& input_signal) {
+        std::vector<std::complex<double>> output_signal(input_signal.size());
+        double current_phase = 0.0;
+        
+        error_total = 0;
+        out_iir_total = 0;
+        in_vco_total = 0;
+
+        for (size_t i = 0; i < input_signal.size(); ++i) {
+            // Multiply input by feedback signal
+            std::complex<double> multiplied = input_signal[i] * feedback;
+
+            // Filter signal
+            std::complex<double> filtered = srrc_filter.filterSample(multiplied);
+
+            // Output best-guess corrected I/Q components
+            output_signal[i] = filtered;
+
+            // Generate limited components
+            std::complex<double> limited = limiter(filtered);
+            
+            // IIR Filter
+            double error_real = (limited.real() > 0 ? 1.0 : -1.0) * limited.imag();
+            double error_imag = (limited.imag() > 0 ? 1.0 : -1.0) * limited.real();
+            double phase_error = error_real - error_imag;
+            phase_error = 0.5 * (std::abs(phase_error + 1) - std::abs(phase_error - 1));
+
+            freq_error += beta * phase_error;
+            double phase_adjust = alpha * phase_error + freq_error;
+
+            current_phase += 2 * M_PI * carrier_freq / sample_rate + k_factor * phase_adjust;
+            if (current_phase > M_PI) current_phase -= 2 * M_PI;
+            else if (current_phase < -M_PI) current_phase += 2 * M_PI;
+
+            // Generate feedback signal for next iteration
+            double feedback_real = std::cos(current_phase);
+            double feedback_imag = -std::sin(current_phase);
+            feedback = std::complex<double>(feedback_real, feedback_imag);
+        }
+
+        return output_signal;
+    }
+
+private:
+    double carrier_freq;
+    double sample_rate;
+    double k_factor;
+    double prev_in_iir;
+    double prev_out_iir;
+    double prev_in_vco;
+    std::complex<double> feedback;
+    double error_total;
+    double out_iir_total;
+    double in_vco_total;
+    SRRCFilter srrc_filter;
+    double vco_gain;
+    double alpha;
+    double beta;
+    double freq_error;
+
+    std::complex<double> limiter(const std::complex<double>& sample) const {
+        double limited_I = std::clamp(sample.real(), -1.0, 1.0);
+        double limited_Q = std::clamp(sample.imag(), -1.0, 1.0);
+        return std::complex<double>(limited_I, limited_Q);
+    }
+};

include/utils/filters.h

@@ -0,0 +1,151 @@
+#ifndef FILTERS_H
+#define FILTERS_H
+
+#include <cmath>
+#include <cstdint>
+#include <fftw3.h>
+#include <numeric>
+#include <vector>
+
+class TapGenerators {
+    public:
+    std::vector<double> generateSRRCTaps(size_t num_taps, double sample_rate, double symbol_rate, double rolloff) const {
+        std::vector<double> 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;
+
+        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);
+
+            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 {
+                taps[i] = (4.0 * rolloff / (M_PI * std::sqrt(T))) * (cos_part / denominator);
+            }
+
+            taps[i] *= sinc_part;
+        }
+
+        // Normalize filter taps
+        double sum = std::accumulate(taps.begin(), taps.end(), 0.0);
+        for (auto& tap : taps) {
+            tap /= sum;
+        }
+
+        return taps;
+    }
+
+    std::vector<double> generateLowpassTaps(size_t num_taps, double cutoff_freq, double sample_rate) const {
+        std::vector<double> 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;
+    }
+};
+
+class Filter {
+    public:
+    Filter(const std::vector<double>& _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<double>(sample, 0.0);
+        double filtered_val = 0.0;
+        size_t idx = buffer_index;
+
+        for (size_t j = 0; j < filter_taps.size(); j++) {
+            filtered_val += filter_taps[j] * buffer[idx].real();
+            if (idx == 0) {
+                idx = buffer.size() - 1;
+            } else {
+                idx--;
+            }
+        }
+
+        buffer_index = (buffer_index + 1) % buffer.size();
+        return filtered_val;
+    }
+
+    std::complex<double> filterSample(const std::complex<double> sample) {
+        buffer[buffer_index] = sample;
+        std::complex<double> filtered_val = std::complex<double>(0.0, 0.0);
+        size_t idx = buffer_index;
+    
+        for (size_t j = 0; j < filter_taps.size(); j++) {
+            filtered_val += filter_taps[j] * buffer[idx];
+            if (idx == 0) {
+                idx = buffer.size() - 1;
+            } else {
+                idx--;
+            }
+        }
+    
+        buffer_index = (buffer_index + 1) % buffer.size();
+        return filtered_val;
+    }
+
+    std::vector<double> applyFilter(const std::vector<double>& signal) {
+        std::vector<double> filtered_signal(signal.size(), 0.0);
+
+        // Convolve the signal with the filter taps
+        for (size_t i = 0; i < signal.size(); ++i) {
+            filtered_signal[i] = filterSample(signal[i]);
+        }
+
+        return filtered_signal;
+    }
+
+    std::vector<std::complex<double>> applyFilter(const std::vector<std::complex<double>>& signal) {
+        std::vector<std::complex<double>> filtered_signal(signal.size(), std::complex<double>(0.0, 0.0));
+
+        // Convolve the signal with the filter taps
+        for (size_t i = 0; i < signal.size(); ++i) {
+            filtered_signal[i] = filterSample(signal[i]);
+        }
+
+        return filtered_signal;
+    }
+
+    private:
+    std::vector<double> filter_taps;
+    std::vector<std::complex<double>> buffer;
+    size_t buffer_index;
+};
+
+class SRRCFilter : public Filter {
+    public:
+    SRRCFilter(const size_t num_taps, const double sample_rate, const double symbol_rate, const double rolloff) : Filter(TapGenerators().generateSRRCTaps(num_taps, sample_rate, symbol_rate, rolloff)) {}
+};
+
+class LowPassFilter : public Filter {
+    public:
+    LowPassFilter(const size_t num_taps, const double cutoff_freq, const double sample_rate) : Filter(TapGenerators().generateLowpassTaps(num_taps, cutoff_freq, sample_rate)) {}
+};
+
+#endif /* FILTERS_H */

include/utils/wattersonchannel.h

@@ -0,0 +1,410 @@
+#ifndef WATTERSONCHANNEL_H
+#define WATTERSONCHANNEL_H
+
+#include <iostream>
+#include <complex>
+#include <vector>
+#include <cmath>
+#include <random>
+#include <algorithm>
+#include <functional>
+#include <fftw3.h> // 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<double>& inputSignal, std::vector<double>& outputSignal);
+
+private:
+    double Fs; // Sample rate
+    double Rs; // Symbol rate
+    double delaySpread; // Delay spread in seconds
+    std::vector<int> delays = {0, L};
+    double fadingBandwidth; // Fading bandwidth d in Hz
+    double SNRdB; // SNR in dB
+    int L; // Length of the simulated channel
+    std::vector<double> f_jt; // Filter impulse response
+    std::vector<std::vector<std::complex<double>>> 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<std::vector<double>> wgnFadingReal; // WGN samples for fading (double part)
+    std::vector<std::vector<double>> wgnFadingImag; // WGN samples for fading (imaginary part)
+    std::vector<std::complex<double>> 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<std::complex<double>>& x_i);
+    void generateWGN(std::vector<double>& wgn, int numSamples);
+    void resampleFadingTapGains();
+    void hilbertTransform(const std::vector<double>& input, std::vector<std::complex<double>>& 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<int>(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<int>(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<int>(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<double>(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<double> 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<double> wgnRealPadded(fftSize, 0.0);
+            std::vector<double> 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<double>(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<std::complex<double>> resampled_h(numOutputSamples);
+        for (int i = 0; i < numOutputSamples; ++i) {
+            double t = i * (1.0 / Fs);
+            double index = t * fadingSampleRate;
+            int idx = static_cast<int>(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<double>(0.0, 0.0);
+            }
+        }
+        h_j[pathIndex] = std::move(resampled_h);
+    }
+}
+
+void WattersonChannel::generateNoise(const std::vector<std::complex<double>>& 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<double> combinedGain = std::complex<double>(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<double> 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<double>(realPart, imagPart);
+    }
+}
+
+void WattersonChannel::generateWGN(std::vector<double>& wgn, int numSamples)
+{
+    wgn.resize(numSamples);
+
+    std::normal_distribution<double> 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<double>& input, std::vector<std::complex<double>>& 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<double>(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<double>& inputSignal, std::vector<double>& outputSignal)
+{
+    // Apply Hilbert transform to input signal to get complex x_i
+    std::vector<std::complex<double>> x_i;
+    hilbertTransform(inputSignal, x_i);
+
+    generateNoise(x_i);
+
+    // Process the signal through the channel
+    std::vector<std::complex<double>> 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<double> 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

main.cpp

@@ -1,52 +1,242 @@
-#include <bitset>
-#include <fstream>
+// main.cpp
+
 #include <iostream>
-#include <string>
-#include <sndfile.h>
 #include <vector>
+#include <cmath>
+#include <complex>
+#include <random>
+#include <sndfile.h> // For WAV file handling
 
+// GNU Radio headers
+#include <gnuradio/top_block.h>
+#include <gnuradio/blocks/vector_source.h>
+#include <gnuradio/blocks/vector_sink.h>
+#include <gnuradio/blocks/wavfile_sink.h>
+#include <gnuradio/blocks/wavfile_source.h>
+#include <gnuradio/blocks/multiply.h>
+#include <gnuradio/blocks/complex_to_real.h>
+#include <gnuradio/blocks/add_blk.h>
+#include <gnuradio/analog/sig_source.h>
+#include <gnuradio/analog/noise_source.h>
+#include <gnuradio/filter/hilbert_fc.h>
+#include <gnuradio/channels/selective_fading_model.h>
+
+// Include your ModemController and BitStream classes
 #include "ModemController.h"
+#include "bitstream.h"
 
-int main() {
-    // Sample test data
-    std::string sample_string = "The quick brown fox jumps over the lazy dog 1234567890";
-    std::vector<uint8_t> sample_data(sample_string.begin(), sample_string.end());
+// 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::mt19937 gen(seed);
+    std::bernoulli_distribution dist(p);
 
-    // Convert sample data to a BitStream object
-    BitStream bitstream(sample_data, sample_data.size() * 8);
+    for (size_t i = 0; i < length * 8; ++i) {
+        random_data.putBit(dist(gen));
+    }
+    return random_data;
+}
 
-    // Configuration for modem
-    size_t baud_rate = 150;
-    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, bitstream);
-
-    const char* file_name = "modulated_signal_150bps_longinterleave.wav";
-
-    // Perform transmit operation to generate modulated signal
-    std::vector<int16_t> modulated_signal = modem.transmit();
-
-    // 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<int16_t>& data, float sample_rate) {
     SF_INFO sfinfo;
     sfinfo.channels = 1;
-    sfinfo.samplerate = 48000;
+    sfinfo.samplerate = static_cast<int>(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<sf_count_t>(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<int16_t> 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<float> 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<float> delays = {0.0f};
+    std::vector<float> 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<float> 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<int16_t> 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<int16_t>(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;
-}
+}

tests/CMakeLists.txt

@@ -0,0 +1,19 @@
+list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake")
+
+# Find the installed gtest package
+find_package(GTest REQUIRED)
+find_package(SndFile REQUIRED)
+find_package(FFTW3 REQUIRED)
+
+# Add test executable
+add_executable(PSKModulatorTest PSKModulatorTests.cpp)
+
+# Link the test executable with the GTest libraries
+target_link_libraries(PSKModulatorTest GTest::GTest GTest::Main FFTW3::fftw3 SndFile::sndfile)
+
+# Enable C++17 standard
+set_target_properties(PSKModulatorTest PROPERTIES CXX_STANDARD 17)
+
+# Add test cases
+include(GoogleTest)
+gtest_discover_tests(PSKModulatorTest)

tests/FSKModulatorTests.cpp

@@ -1,23 +0,0 @@
-#include "gtest/gtest.h"
-#include "FSKModulator.h"
-#include <vector>
-
-TEST(FSKModulatorTest, SignalLength) {
-    using namespace milstd;
-
-    // Parameters
-    FSKModulator modulator(FSKModulator::ShiftType::NarrowShift, 75.0, 8000.0);
-
-    // Input data bits
-    std::vector<uint8_t> dataBits = {1, 0, 1, 1, 0};
-
-    // Modulate the data
-    std::vector<double> signal = modulator.modulate(dataBits);
-
-    // Calculate expected number of samples
-    size_t samplesPerSymbol = static_cast<size_t>(modulator.getSampleRate() * modulator.getSymbolDuration());
-    size_t expectedSamples = dataBits.size() * samplesPerSymbol;
-
-    // Verify signal length
-    EXPECT_EQ(signal.size(), expectedSamples);
-}

tests/PSKModulatorTests.cpp

@@ -0,0 +1,68 @@
+#include <gtest/gtest.h>
+#include "modulation/PSKModulator.h"
+
+// Fixture for PSK Modulator tests
+class PSKModulatorTest : public ::testing::Test {
+protected:
+    double sample_rate = 48000;
+    size_t num_taps = 48;
+    bool is_frequency_hopping = false;
+    PSKModulator modulator{sample_rate, is_frequency_hopping, num_taps};
+
+    std::vector<uint8_t> symbols = {0, 3, 5, 7};
+};
+
+TEST_F(PSKModulatorTest, ModulationOutputLength) {
+    auto signal = modulator.modulate(symbols);
+
+    size_t expected_length = symbols.size() * (sample_rate / SYMBOL_RATE);
+    ASSERT_EQ(signal.size(), expected_length);
+
+    for (auto& sample : signal) {
+        EXPECT_GE(sample, -32768);
+        EXPECT_LE(sample, 32767);
+    }
+}
+
+TEST_F(PSKModulatorTest, DemodulationOutput) {
+    auto passband_signal = modulator.modulate(symbols);
+
+    // Debug: Print modulated passband signal
+    std::cout << "Modulated Passband Signal: ";
+    for (const auto& sample : passband_signal) {
+        std::cout << sample << " ";
+    }
+    std::cout << std::endl;
+
+    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: ";
+    for (const auto& symbol : decoded_symbols) {
+        std::cout << (int)symbol << " ";
+    }
+    std::cout << std::endl;
+
+    // Debug: Print expected symbols
+    std::cout << "Expected Symbols: ";
+    for (const auto& symbol : symbols) {
+        std::cout << (int)symbol << " ";
+    }
+    std::cout << std::endl;
+
+    ASSERT_EQ(symbols.size(), decoded_symbols.size());
+
+    for (size_t i = 0; i < symbols.size(); i++) {
+        EXPECT_EQ(symbols[i], decoded_symbols[i]) << " at index " << i;
+    }
+}
+
+
+TEST_F(PSKModulatorTest, InvalidSymbolInput) {
+    std::vector<uint8_t> invalid_symbols = {0, 8, 9};
+
+    EXPECT_THROW(modulator.modulate(invalid_symbols), std::out_of_range);
+}