From 60d6f6db7d085f445f05038e48595bfae7d22589 Mon Sep 17 00:00:00 2001 From: RecklessAndFeckless Date: Sat, 12 Oct 2024 17:24:21 -0400 Subject: [PATCH] SRRC filter goes on the IQ components, not the bandpass signal; properly defining filter taps based on freq response --- CMakeLists.txt | 2 + cmake/FindFFTW3.cmake | 60 ++++++++ include/encoder/ModemController.h | 10 +- include/modulation/PSKModulator.h | 222 ++++++++++++++++++------------ main.cpp | 4 +- 5 files changed, 202 insertions(+), 96 deletions(-) create mode 100644 cmake/FindFFTW3.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index fde67e4..f7006d0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,12 +19,14 @@ set(SOURCES main.cpp) # Link with libsndfile list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake") find_package(SndFile REQUIRED) +find_package(FFTW3 REQUIRED) # Add executable add_executable(MILSTD110C ${SOURCES}) # Link executable with libsndfile library target_link_libraries(MILSTD110C SndFile::sndfile) +target_link_libraries(MILSTD110C FFTW3::fftw3) # Debug and Release Build Types set(CMAKE_CONFIGURATION_TYPES "Debug;Release" CACHE STRING "" FORCE) diff --git a/cmake/FindFFTW3.cmake b/cmake/FindFFTW3.cmake new file mode 100644 index 0000000..b7ab95e --- /dev/null +++ b/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) diff --git a/include/encoder/ModemController.h b/include/encoder/ModemController.h index c572d69..d602575 100644 --- a/include/encoder/ModemController.h +++ b/include/encoder/ModemController.h @@ -45,13 +45,13 @@ public: 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), + 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, 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. diff --git a/include/modulation/PSKModulator.h b/include/modulation/PSKModulator.h index 5aaa4d1..74c717e 100644 --- a/include/modulation/PSKModulator.h +++ b/include/modulation/PSKModulator.h @@ -1,136 +1,180 @@ #ifndef PSK_MODULATOR_H #define PSK_MODULATOR_H -#include -#include -#include -#include -#include + #include +#include +#include +#include +#include +#include +#include +#include + +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, 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(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(sample_rate / SYMBOL_RATE)), num_taps(_num_taps) { + initializeSymbolMap(); } - std::vector modulate(const std::vector& symbols) { - std::vector> modulated_signal; + std::vector modulate(const std::vector& symbols) const { + std::vector baseband_I(symbols.size() * samples_per_symbol); + std::vector baseband_Q(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 target_symbol = symbolMap[symbol]; + const std::complex 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_I[symbol_index * samples_per_symbol + i] = target_symbol.real(); + baseband_Q[symbol_index * samples_per_symbol + i] = target_symbol.imag(); } + + 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 + auto filter_taps = generateSRRCTaps(num_taps, sample_rate, SYMBOL_RATE, ROLLOFF_FACTOR); + auto filtered_I = applyFilter(baseband_I, filter_taps); + auto filtered_Q = applyFilter(baseband_Q, filter_taps); - // 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 and apply gain for audio output - std::vector combined_signal; - for (auto& sample : filtered_signal) { - int16_t combined_sample = static_cast(std::clamp(gain * (sample.real() + sample.imag()), -32768.0, 32767.0)); - combined_signal.push_back(combined_sample); + std::vector> baseband_components(filtered_I.size()); + for (size_t i = 0; i < filtered_I.size(); i++) { + std::complex baseband_component = {filtered_I[i], filtered_Q[i]}; + baseband_components[i] = baseband_component; } - return combined_signal; + // Combine the I and Q components + std::vector passband_signal; + passband_signal.reserve(baseband_components.size()); + + double carrier_phase = 0.0; + double carrier_phase_increment = 2 * M_PI * CARRIER_FREQ / sample_rate; + for (const auto& sample : baseband_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 + carrier_phase += carrier_phase_increment; + if (carrier_phase >= 2 * M_PI) + carrier_phase -= 2 * M_PI; + } + + std::vector final_signal; + final_signal.reserve(passband_signal.size()); + + for (const auto& sample : passband_signal) { + int16_t value = static_cast(sample); + value = std::clamp(value, (int16_t)-32768, (int16_t)32767); + final_signal.emplace_back(value); + } + + return final_signal; } - std::vector sqrtRaisedCosineFilter(size_t num_taps, double symbol_rate) { - double rolloff = 0.35; // Fixed rolloff factor as per specification - std::vector 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(num_taps - 1) / 2.0; +private: + 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. + const size_t num_taps; + std::vector> symbolMap; ///< The mapping of tribit symbols to I/Q components. + const bool is_frequency_hopping; ///< Whether to use frequency hopping methods. Not implemented (yet?) - 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; + static inline double validateSampleRate(const double rate) { + if (rate <= 2 * CARRIER_FREQ) { + throw std::out_of_range("Sample rate must be above the Nyquist frequency (PSKModulator.h)"); + } + return rate; + } + + inline void initializeSymbolMap() { + symbolMap = { + {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) + }; + } + + std::vector generateSRRCTaps(size_t num_taps,double sample_rate, double symbol_rate, double rolloff) const { + std::vector freq_response(num_taps, 0.0); + std::vector taps(num_taps); + + double fn = symbol_rate / 2.0; + double f_step = sample_rate / num_taps; + + 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))); } 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; + freq_response[i] = 0.0; } - norm_factor += filter_taps[i] * filter_taps[i]; } - 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; + 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> applyFilter(const std::vector>& signal, const std::vector& filter_taps) { - std::vector> filtered_signal(signal.size()); + std::vector applyFilter(const std::vector& signal, const std::vector& filter_taps) const { + std::vector filtered_signal(signal.size(), 0.0); 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; + double filtered_val = 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 signal_index = i + j - half_filter_length; + if (signal_index < signal.size()) { + filtered_val += filter_taps[j] * signal[signal_index]; } } - filtered_signal[i] = std::complex(filtered_i, filtered_q); + filtered_signal[i] = filtered_val; } return filtered_signal; } - -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. - size_t samples_per_symbol; ///< Number of samples per symbol, calculated to match symbol duration with cycle. - size_t symbol_rate; - std::vector> symbolMap; ///< The mapping of tribit symbols to I/Q components. - - 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) - }; - } }; #endif \ No newline at end of file diff --git a/main.cpp b/main.cpp index ecaed72..a75b919 100644 --- a/main.cpp +++ b/main.cpp @@ -16,7 +16,7 @@ int main() { BitStream input_data(sample_data, sample_data.size() * 8); // Configuration for modem - size_t baud_rate = 2400; + size_t baud_rate = 300; bool is_voice = false; // False indicates data mode bool is_frequency_hopping = false; // Fixed frequency operation size_t interleave_setting = 1; // Short interleave @@ -24,7 +24,7 @@ int main() { // Create ModemController instance ModemController modem(baud_rate, is_voice, is_frequency_hopping, interleave_setting, input_data); - const char* file_name = "modulated_signal_2400bps_voice.wav"; + const char* file_name = "modulated_signal_600bps_shortinterleave.wav"; // Perform transmit operation to generate modulated signal std::vector modulated_signal = modem.transmit();