From 9acafbbee999d397f219beb1536a77d7b63975d0 Mon Sep 17 00:00:00 2001 From: AlexandreRouma Date: Wed, 30 Jul 2025 15:57:07 +0200 Subject: [PATCH] progress --- .gitignore | 3 +- CMakeLists.txt | 38 +++++++- dsp/block.cpp | 32 +++++-- dsp/block.h | 77 +++++++++++----- dsp/complex.h | 8 +- dsp/demod/fm.cpp | 94 +++++++++++++++++++ dsp/demod/fm.h | 100 ++++++++++++++++++++ dsp/filter/fir.cpp | 0 dsp/filter/fir.h | 77 ++++++++++++++++ dsp/mailbox.h | 72 --------------- dsp/processor.h | 156 +++++++++++++++++++++++++++++++- dsp/stereo.h | 90 ------------------ dsp/stream.cpp | 92 ++++++++++++++----- dsp/stream.h | 79 +++++++--------- dsp/taps.cpp | 28 +++--- dsp/taps.h | 26 ++++-- dsp/taps/low_pass.h | 10 +- dsp/window.cpp | 16 ++++ dsp/window.h | 34 +++++++ dsp/window/rectangular.h | 14 +++ license | 1 + main.cpp | 48 ++++++++++ nuttall.txt | 3 + src/main.cpp | 14 --- todo.txt | 9 ++ tools/filt_spec_reg.py | 191 +++++++++++++++++++++++++++++++++++++++ tools/spec_fit.py | 33 +++++++ 27 files changed, 1040 insertions(+), 305 deletions(-) create mode 100644 dsp/demod/fm.cpp create mode 100644 dsp/demod/fm.h create mode 100644 dsp/filter/fir.cpp create mode 100644 dsp/filter/fir.h delete mode 100644 dsp/mailbox.h delete mode 100644 dsp/stereo.h create mode 100644 dsp/window.cpp create mode 100644 dsp/window.h create mode 100644 dsp/window/rectangular.h create mode 100644 license create mode 100644 main.cpp create mode 100644 nuttall.txt create mode 100644 todo.txt create mode 100644 tools/filt_spec_reg.py create mode 100644 tools/spec_fit.py diff --git a/.gitignore b/.gitignore index 608b1d4..c30f610 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ .vscode/ -build/ \ No newline at end of file +build/ +.old/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 5def151..c52d7bf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,6 +5,40 @@ file(GLOB_RECURSE SRC "src/*.cpp" "dsp/*.cpp") add_executable(${PROJECT_NAME} ${SRC}) -target_include_directories(${PROJECT_NAME} PRIVATE "dsp/") +target_include_directories(${PROJECT_NAME} PRIVATE "/") -target_compile_options(${PROJECT_NAME} PRIVATE /std:c++20) \ No newline at end of file +target_compile_options(${PROJECT_NAME} PRIVATE /std:c++20) + + +if (MSVC) + # Lib path + target_link_directories(${PROJECT_NAME} PUBLIC "C:/Program Files/PothosSDR/lib/") + + # Misc headers + target_include_directories(${PROJECT_NAME} PUBLIC "C:/Program Files/PothosSDR/include/") + + # Volk + target_link_libraries(${PROJECT_NAME} PUBLIC volk) + +else() + find_package(PkgConfig) + + pkg_check_modules(VOLK REQUIRED volk) + + target_include_directories(${PROJECT_NAME} PUBLIC + ${VOLK_INCLUDE_DIRS} + ) + + target_link_directories(${PROJECT_NAME} PUBLIC + ${VOLK_LIBRARY_DIRS} + ) + + target_link_libraries(${PROJECT_NAME} PUBLIC + ${VOLK_LIBRARIES} + ) + + if (${CMAKE_SYSTEM_NAME} MATCHES "Linux") + target_link_libraries(${PROJECT_NAME} PUBLIC stdc++fs pthread) + endif () + +endif () \ No newline at end of file diff --git a/dsp/block.cpp b/dsp/block.cpp index e9f360c..3a86c6c 100644 --- a/dsp/block.cpp +++ b/dsp/block.cpp @@ -1,8 +1,12 @@ #include "block.h" +#include namespace dsp { - Block::Block() { + Block::Block() {} + Block::~Block() { + // Stop the worker + stop(); } void Block::start() { @@ -53,7 +57,7 @@ namespace dsp { _running = false; } - bool Block::running() { + bool Block::running() const { // Acquire worker variables std::lock_guard lck(workerMtx); @@ -61,7 +65,7 @@ namespace dsp { return _running; } - void Block::registerInput(Signaler* input) { + void Block::registerInput(StopNotifier* input) { // Acquire worker variables std::lock_guard lck(workerMtx); @@ -69,14 +73,19 @@ namespace dsp { inputs.push_back(input); } - void Block::unregisterInput(Signaler* input) { + void Block::unregisterInput(StopNotifier* input) { // Acquire worker variables std::lock_guard lck(workerMtx); - // TODO + // Find the notifier + auto it = std::find(inputs.begin(), inputs.end(), input); + if (it == inputs.end()) { return; } + + // Remove it from the list + inputs.erase(it); } - void Block::registerOutput(Signaler* output) { + void Block::registerOutput(StopNotifier* output) { // Acquire worker variables std::lock_guard lck(workerMtx); @@ -84,15 +93,20 @@ namespace dsp { outputs.push_back(output); } - void Block::unregisterOutput(Signaler* output) { + void Block::unregisterOutput(StopNotifier* output) { // Acquire worker variables std::lock_guard lck(workerMtx); - // TODO + // Find the notifier + auto it = std::find(outputs.begin(), outputs.end(), output); + if (it == outputs.end()) { return; } + + // Remove it from the list + inputs.erase(it); } void Block::worker() { // Call the run function repeatedly - while (!run()); + while (run()); } } \ No newline at end of file diff --git a/dsp/block.h b/dsp/block.h index b726946..7b8bdd4 100644 --- a/dsp/block.h +++ b/dsp/block.h @@ -1,16 +1,46 @@ #pragma once #include +#include #include -#include "stream.h" namespace dsp { + /** + * Interface to be used by any blocking class (stream, mailbox, etc...) to allow cancelling an operation. + */ + class StopNotifier { + public: + /** + * Notify the sending thread that it should stop. + */ + virtual void stopSender() = 0; + + /** + * Clear the sender stop flag to allow restarting the sender thread. + */ + virtual void clearSendStop() = 0; + + /** + * Notify the receiving thread that it should stop. + */ + virtual void stopReceiver() = 0; + + /** + * Clear the receiver stop flag to allow restarting the sender thread. + */ + virtual void clearRecvStop() = 0; + }; + /** * General DSP block class handling the worker thread start/stop operations. - * All input and output streams of the derived blocks must be registered using the appropriate functions. + * All input and output stop notifiers (usually streams) of the derived blocks must be registered using the appropriate functions. + * This class is thread-safe. */ class Block { public: + // Default constructor Block(); + + // Destructor virtual ~Block(); /** @@ -26,43 +56,40 @@ namespace dsp { /** * Check wether or not the block's worker thread is running. */ - bool running(); + bool running() const; protected: /** - * Register an input stream. - * @param input Input stream to register. + * Register an input stop notifier. + * @param input Input stop notifier to register. */ - void registerInput(Signaler* input); + void registerInput(StopNotifier* input); /** - * Unregister an input stream. - * @param input Input stream to unregister. + * Unregister an input stop notifier. + * @param input Input stop notifier to unregister. */ - void unregisterInput(Signaler* input); + void unregisterInput(StopNotifier* input); /** - * Register an output stream. - * @param input Output stream to register. + * Register an output stop notifier. + * @param input Output stop notifier to register. */ - void registerOutput(Signaler* output); + void registerOutput(StopNotifier* output); /** - * Unregister an output stream. - * @param input Output stream to unregister. + * Unregister an output stop notifier. + * @param input Output stop notifier to unregister. */ - void unregisterOutput(Signaler* output); + void unregisterOutput(StopNotifier* output); + + // TODO: Pause/Resume for when inputs/outputs change to avoid needing to restart a thread /** * Run the DSP code. - * @return True if the worker thread should stop. + * @return False if the worker thread should stop. True otherwise. */ - virtual bool run(); - - /** - * Mutex to be used for block settings. - */ - std::recursive_mutex settingsMtx; + virtual bool run() = 0; private: /** @@ -71,10 +98,10 @@ namespace dsp { void worker(); // Worker variables - std::mutex workerMtx; + mutable std::mutex workerMtx; std::thread workerThread; - std::vector inputs; - std::vector outputs; + std::vector inputs; + std::vector outputs; bool _running = false; }; } \ No newline at end of file diff --git a/dsp/complex.h b/dsp/complex.h index f105bd5..fc92d85 100644 --- a/dsp/complex.h +++ b/dsp/complex.h @@ -4,6 +4,7 @@ namespace dsp { /** * Complex 32bit floating-point number. + * This struct is NOT thread-safe. */ struct Complex { // TODO: Allow construction from a float @@ -29,7 +30,8 @@ namespace dsp { return sqrtf(re*re + im*im); } - void operator=(float b) { + // Assignment operator + inline void operator=(float b) { re = b; } @@ -95,10 +97,10 @@ inline constexpr dsp::Complex operator/(const dsp::Complex& a, const dsp::Comple return dsp::Complex{ (a.re*b.re + a.im*b.im) / denom, (a.im*b.re - a.re*b.im) / denom }; } -inline constexpr dsp::Complex operator""_j(unsigned long long value) { +inline constexpr dsp::Complex operator""j(unsigned long long value) { return dsp::Complex{ 0.0f, (float)value }; } -inline constexpr dsp::Complex operator""_j(long double value) { +inline constexpr dsp::Complex operator""j(long double value) { return dsp::Complex{ 0.0f, (float)value }; } \ No newline at end of file diff --git a/dsp/demod/fm.cpp b/dsp/demod/fm.cpp new file mode 100644 index 0000000..e123bcb --- /dev/null +++ b/dsp/demod/fm.cpp @@ -0,0 +1,94 @@ +#include "fm.h" +#define _USE_MATH_DEFINES +#include + +namespace dsp::demod { + FM_s::FM_s() {} + + FM_s::FM_s(float deviation, float samplerate) { + // Save the settings + this->deviation = deviation; + this->samplerate = samplerate; + + // Update the normalization factor + normFact = samplerate / (2.0f * M_PI * deviation); + } + + float FM_s::getDeviation() const { + // Acquire the settings mutex + std::lock_guard lck(settingsMtx); + + // Return the deviation + return deviation; + } + + void FM_s::setDeviation(float deviation) { + // Acquire the settings mutex + std::lock_guard lck(settingsMtx); + + // Update the deviation + this->deviation = deviation; + + // Update the normalization factor + normFact = samplerate / (2.0f * M_PI * deviation); + } + + float FM_s::getSamplerate() const { + // Acquire the settings mutex + std::lock_guard lck(settingsMtx); + + // Return the samplerate + return samplerate; + } + + void FM_s::setSamplerate(float samplerate) { + // Acquire the settings mutex + std::lock_guard lck(settingsMtx); + + // Update the samplerate + this->samplerate = samplerate; + + // Update the normalization factor + normFact = samplerate / (2.0f * M_PI * deviation); + } + + void FM_s::reset() { + // Acquire the settings mutex + std::lock_guard lck(settingsMtx); + + // Set the current phase to zero + phase = 0.0f; + } + + int FM_s::process(const Complex* in, float* out, int count) { + for (int i = 0; i < count; i++) { + // Get the current phase + float cphase = in[i].phase(); + + // Compute the difference with the last phase + // TODO + //out[i] = math::normalizePhase(cphase - phase) * normFact; + + // Save the current phase for the next iteration + phase = cphase; + } + return count; + } + + FM::FM() {} + + FM::FM(float deviation, float samplerate) : + FM_s(deviation, samplerate), + ProcessorAsync(this) + {} + + FM::FM(Stream* in, float deviation, float samplerate) : + FM_s(deviation, samplerate), + ProcessorAsync(this, in) + {} + + FM::FM(Stream* in, Stream* out, float deviation, float samplerate) : + FM_s(deviation, samplerate), + ProcessorAsync(this, in, out) + {} +} \ No newline at end of file diff --git a/dsp/demod/fm.h b/dsp/demod/fm.h new file mode 100644 index 0000000..9870afd --- /dev/null +++ b/dsp/demod/fm.h @@ -0,0 +1,100 @@ +#pragma once +#include "../processor.h" +#include "../complex.h" + +namespace dsp::demod { + /** + * FM demodulator (Synchronous). + * This class is thread-safe except for `process()`. + */ + class FM_s : public ProcessorSync { + public: + // Default constructor + FM_s(); + + /** + * Create an FM demodulator. + * @param deviation Deviation of the FM signal in Hz. + * @param samplerate Samplerate of the signal in Hz. + */ + FM_s(float deviation, float samplerate); + + /** + * Get the deviation. + * @return Deviation of the FM signal in Hz. + */ + float getDeviation() const; + + /** + * Set the deviation. + * @param deviation Deviation of the FM signal in Hz. + */ + void setDeviation(float deviation); + + /** + * Get the samplerate. + * @return Samplerate of the signal in Hz. + */ + float getSamplerate() const; + + /** + * Set the samplerate. + * @param deviation Samplerate of the signal in Hz. + */ + void setSamplerate(float samplerate); + + /** + * Reset the state of the demodulator. + */ + void reset(); + + /** + * Demodulate a FM-modulated signal. A lock must be acquired using `getLock()` prior to invoking if settings could be set from another thread. + * @param in Modulated signal buffer. + * @param out Demodulated signal buffer. + * @param count Number of samples in the input buffer. + * @return Number of samples in the output buffer. Will always be equal to the number of samples in the input buffer. + */ + int process(const Complex* in, float* out, int count); + + private: + float deviation = 0.0f; + float samplerate = 0.0f; + float normFact = 1.0f; + float phase = 0.0f; + }; + + /** + * FM demodulator. + * This class is thread-safe. + */ + class FM : public FM_s, public ProcessorAsync { + public: + // Default constructor + FM(); + + /** + * Create an FM demodulator. + * @param deviation Deviation of the FM signal in Hz. + * @param samplerate Samplerate of the signal in Hz. + */ + FM(float deviation, float samplerate); + + /** + * Create an FM demodulator. + * @param in Modulated signal stream. + * @param deviation Deviation of the FM signal in Hz. + * @param samplerate Samplerate of the signal in Hz. + */ + FM(Stream* in, float deviation, float samplerate); + + /** + * Create an FM demodulator. + * @param in Modulated signal stream. + * @param out Demodulated signal stream. + * @param deviation Deviation of the FM signal in Hz. + * @param samplerate Samplerate of the signal in Hz. + */ + FM(Stream* in, Stream* out, float deviation, float samplerate); + }; +} \ No newline at end of file diff --git a/dsp/filter/fir.cpp b/dsp/filter/fir.cpp new file mode 100644 index 0000000..e69de29 diff --git a/dsp/filter/fir.h b/dsp/filter/fir.h new file mode 100644 index 0000000..8c0435b --- /dev/null +++ b/dsp/filter/fir.h @@ -0,0 +1,77 @@ +#pragma once +#include "../processor.h" +#include "../taps.h" + +namespace dsp::filter { + /** + * Finite Inpulse Response filter (Synchronous). + * This class is thread-safe except for `process()`. + */ + template + class FIR_s : public ProcessorSync { + public: + /** + * Create a FIR filter. + * @param taps Filter taps. + */ + FIR_s(const Taps& taps); + + /** + * Get the filter taps. + * @return Filter taps. + */ + Taps getTaps(); + + /** + * Set the filter taps. + * @param taps Filter taps. + */ + void setTaps(const Taps& taps); + + /** + * Reset the state of the filter. + */ + void reset(); + + /** + * Filter samples. + * @param in Input sample buffer. + * @param out Filtered sample buffer. + * @param count Number of samples in the input buffer. + * @return Number of samples in the output buffer. Will always be equal to the number of samples in the input buffer. + */ + int process(const SAMP_T* in, SAMP_T* out, int count); + + private: + Taps taps; + }; + + /** + * Finite Inpulse Response filter. + * This class is thread-safe except for `process()`. + */ + template + class FIR : public FIR_s, public ProcessorAsync { + public: + /** + * Create a FIR filter. + * @param taps Filter taps. + */ + FIR(const Taps& taps); + + /** + * Create a FIR filter. + * @param in Input samples. + * @param taps Filter taps. + */ + FIR(Stream* in, const Taps& taps); + + /** + * Create a FIR filter. + * @param in Input samples. + * @param in Filtered samples. + * @param taps Filter taps. + */ + FIR(Stream* in, Stream* out, const Taps& taps); + }; +} \ No newline at end of file diff --git a/dsp/mailbox.h b/dsp/mailbox.h deleted file mode 100644 index f7098c2..0000000 --- a/dsp/mailbox.h +++ /dev/null @@ -1,72 +0,0 @@ -#pragma once -#include "stream.h" - -namespace gui { - /** - * Mailboxes allow to exchange objects between two threads. - */ - template - class Mailbox : public Signaler { - public: - /** - * Create a mailbox object. - */ - Mailbox(); - - ~Mailbox(); - - /** - * Notify the sending thread that it should stop. - */ - void stopSender(); - - /** - * Notify the receiving thread that it should stop. - */ - void stopReceiver(); - - /** - * Send an object. - * @return True if the sender thread must exist, false otherwise. - */ - bool send(); - - /** - * Wait for an object. May also return in case of a signal to exit. ack() must be called as soon as the received object has been processed. - * @return True if the receiver thread must exist, false otherwise. - */ - bool recv(); - - /** - * Acknowledge reception and processing of the samples. Allows sender thread to send a new buffer. - */ - void ack(); - - /** - * Get the sending object. - * @return Sending object. - */ - T* getSendObject() { return sendObj; } - - /** - * Get the receiving object. - * @return Sending object. - */ - const T* getRecvObject() { return recvObj; } - - private: - // Sender variables - std::condition_variable sendCV; - std::mutex sendMtx; - bool canSend = true; - bool stopSend = false; - T* sendObj; - - // Receiver variables - std::condition_variable recvCV; - std::mutex recvMtx; - bool available = false; - bool stopRecv = false; - T* recvObj; - }; -} \ No newline at end of file diff --git a/dsp/processor.h b/dsp/processor.h index 86b3161..bef6de0 100644 --- a/dsp/processor.h +++ b/dsp/processor.h @@ -1,12 +1,164 @@ #pragma once +#include #include "block.h" +#include "stream.h" namespace dsp { - template - class Processor : public Block { + // TODO: Deal with the fact that always locking will be slow in hier blocks... + + /** + * Class representing a DSP kernel taking one input and one output. + * Dervied classes must be thread-safe by locking the provided `settingsMtx` mutex in any function changing settings. + */ + template + class ProcessorKernel { public: + // Destructor + virtual ~ProcessorKernel() {} + // TODO: Copy/Move Constructor/Operator + + /** + * Acquire a lock to the settings of the kernel. Mandatory if settings are changed in a different thread than the processing function. + * @return Lock to the settings of the block. + */ + inline std::lock_guard getLock() const { + return std::lock_guard(settingsMtx); + } + + /** + * Process samples. This function is NOT thread-safe. + * @param in Input buffer. + * @param out Output buffer. + * @param count Number of samples in the input buffer. + * @return Number of samples in the output buffer. + */ + virtual int process(const IN* in, OUT* out, int count) = 0; + + protected: + /** + * Mutex to be used for kernel settings. + */ + mutable std::recursive_mutex settingsMtx; + }; + + template + class ProcessorBlock : public Block { + public: + // Default constructor + ProcessorBlock() { + // TODO: Maybe something to do to prevent bad shit if started? + } + + /** + * TODO + */ + ProcessorBlock(ProcessorSync* proc) { + // Save the pointer to the processor + this->proc = proc; + + // Set the streams + setInput(NULL); + setOutput(NULL); + } + + /** + * TODO + */ + ProcessorBlock(ProcessorSync* proc, Stream* in) { + // Save the pointer to the processor + this->proc = proc; + + // Set the streams + setInput(in); + setOutput(NULL); + } + + /** + * TODO + */ + ProcessorBlock(ProcessorSync* proc, Stream* in, Stream* out) { + // Save the pointer to the processor + this->proc = proc; + + // Set the streams + setInput(in); + setOutput(out); + } + + // Destructor + virtual ~ProcessorBlock() {} + + /** + * Set the input stream. + * @param in Input stream. + */ + void setInput(Stream* in) { + // TODO: Lock + + // If the current input if it already was registered + unregisterInput(_in); + + // TODO: Manage allocating and freeing streams + + // Update the input + _in = in; + + // Register the new input + registerInput(_in); + } + + /** + * Set the output stream. + * @param in Output stream. + */ + void setOutput(Stream* out) { + // TODO: Lock + + // If the current output if it already was registered + unregisterOutput(_out); + + // TODO: Manage allocating and freeing streams + + // Update the output + _out = out; + + // Register the new output + registerOutput(_out); + } + + Stream* in() const { + // TODO: Lock + return _in; + } + + Stream* out() const { + // TODO: Lock + return _out; + } + private: + bool run() { + // Read samples + auto bufSet = _in->recv(); + if (!bufSet.samples) { return true; } + // Process samples + { + auto lck = proc->getLock(); + proc->process(NULL/*TODO*/, NULL/*TODO*/, 0/*TODO*/); + } + + // TODO: Write samples + + // TODO + return false; + } + + ProcessorKernel* proc; + Stream* _in = NULL; + Stream* _out = NULL; + bool ownInput = false; + bool ownOutput = false; }; } \ No newline at end of file diff --git a/dsp/stereo.h b/dsp/stereo.h deleted file mode 100644 index 33409fd..0000000 --- a/dsp/stereo.h +++ /dev/null @@ -1,90 +0,0 @@ -#pragma once -#include - -namespace dsp { - /** - * Two 32bit floating-point number representing the left and right channels. - */ - struct Stereo { - // TODO: Allow construction from a float - - void operator=(float b) { - l = b; - r = b; - } - - /** - * Left channel. - */ - float l; - - /** - * Right channel. - */ - float r; - }; - - inline constexpr dsp::Stereo operator+(const dsp::Stereo& a, float b) { - return dsp::Stereo{ a.l + b, a.r + b}; - } - - inline constexpr dsp::Stereo operator+(float a, const dsp::Stereo& b) { - return dsp::Stereo{ a + b.l, a + b.r }; - } - - inline constexpr dsp::Stereo operator-(const dsp::Stereo& a, float b) { - return dsp::Stereo{ a.l - b, a.r - b }; - } - - inline constexpr dsp::Stereo operator-(float a, const dsp::Stereo& b) { - return dsp::Stereo{ a - b.l, a - b.r }; - } - - inline constexpr dsp::Stereo operator*(const dsp::Stereo& a, float b) { - return dsp::Stereo{ a.l*b, a.r*b }; - } - - inline constexpr dsp::Stereo operator*(float a, const dsp::Stereo& b) { - return dsp::Stereo{ a*b.l, a*b.r }; - } - - inline constexpr dsp::Stereo operator/(const dsp::Stereo& a, float b) { - return dsp::Stereo{ a.l/b, a.r/b }; - } - - inline constexpr dsp::Stereo operator/(float a, const dsp::Stereo& b) { - return dsp::Stereo{ a/b.l, a/b.r }; - } - - inline constexpr dsp::Stereo operator+(const dsp::Stereo& a, const dsp::Stereo& b) { - return dsp::Stereo{ a.l + b.l, a.r + b.r }; - } - - inline constexpr dsp::Stereo operator-(const dsp::Stereo& a, const dsp::Stereo& b) { - return dsp::Stereo{ a.l - b.l, a.r - b.r }; - } - - inline constexpr dsp::Stereo operator*(const dsp::Stereo& a, const dsp::Stereo& b) { - return dsp::Stereo{ a.l*b.l, a.r*b.r }; - } - - inline constexpr dsp::Stereo operator/(const dsp::Stereo& a, const dsp::Stereo& b) { - return dsp::Stereo{ a.l/b.l, a.r/b.r }; - } - - inline constexpr dsp::Stereo operator""_L(unsigned long long value) { - return dsp::Stereo{ (float)value, 0.0f }; - } - - inline constexpr dsp::Stereo operator""_L(long double value) { - return dsp::Stereo{ (float)value, 0.0f }; - } - - inline constexpr dsp::Stereo operator""_R(unsigned long long value) { - return dsp::Stereo{ 0.0f, (float)value }; - } - - inline constexpr dsp::Stereo operator""_R(long double value) { - return dsp::Stereo{ 0.0f, (float)value }; - } -} \ No newline at end of file diff --git a/dsp/stream.cpp b/dsp/stream.cpp index 4d48e62..f31507f 100644 --- a/dsp/stream.cpp +++ b/dsp/stream.cpp @@ -1,20 +1,29 @@ -#include "Stream.h" +#include "stream.h" #include "./complex.h" -#include "stereo.h" namespace dsp { template - Stream::Stream(int channels, int bufferSize) { - // Allocate both send and receive buffers aligned by the size of the type - sendBuf = new T[bufferSize, sizeof(T)]; - recvBuf = new T[bufferSize, sizeof(T)]; + Stream::Stream() { + // Allocate both send and receive buffer sets + sendSet = new BufferSet; + recvSet = new BufferSet; + + // Initialize them to zero + sendSet->buffer = new T*; + *sendSet->buffer = NULL; + recvSet->buffer = new T*; + *recvSet->buffer = NULL; + sendSet->capacity = 0; + recvSet->capacity = 0; + sendSet->samples = 0; + recvSet->samples = 0; } template Stream::~Stream() { - // Free both send and receive buffers - delete[] sendBuf; - delete[] recvBuf; + // Free both send and receive buffer sets + delete sendSet; + delete recvSet; } template @@ -66,15 +75,46 @@ namespace dsp { } template - bool Stream::send(int count) { + const BufferSet& Stream::reserve(size_t bufferSize, size_t channels) { // Acquire the sender variables std::unique_lock slck(sendMtx); + // If the capacity is too small or too large, reallocate + if (bufferSize > sendSet->capacity || bufferSize < (sendSet->capacity >> 1) || sendSet->channels != channels) { + // Free all buffers + delete[] sendSet->buffer[0]; + delete[] sendSet->buffer; + + // TODO: Use volk instead + + // Reallocate buffers + sendSet->buffer = new T*[channels]; + T* base = new T[channels * bufferSize]; + for (size_t i = 0; i < channels; i++) { + sendSet->buffer[i] = &base[bufferSize * i]; + } + } + + // Return the send buffer set + return *sendSet; + } + + template + bool Stream::send(size_t count) { + // Acquire the sender variables + std::unique_lock slck(sendMtx); + + // Update the sample count + sendSet->samples = count; + // Wait until the sender can send or is notified it should stop sendCV.wait(slck, [=](){ return canSend || stopSend; }); // If asked to stop, return true - if (stopSend) { return true; } + if (stopSend) { return false; } + + // If trying to send no samples, do nothing + if (!count) { return true; } // Mark that data can no longer be sent canSend = false; @@ -82,38 +122,47 @@ namespace dsp { // Acquire the receiver variables std::unique_lock rlck(recvMtx); - // Swap buffers - T* tmp = sendBuf; - sendBuf = recvBuf; - recvBuf = tmp; + // Swap the buffers sets + BufferSet* tmp = sendSet; + sendSet = recvSet; + recvSet = tmp; // Release the sender variables slck.unlock(); - // Set the number of items that are readable - available = count; + // Set the available flag + available = true; // Release the receiver variables rlck.unlock(); // Notify the receiver thread that there are items available recvCV.notify_all(); + + // Return successfully + return true; } template - int Stream::recv() { + const BufferSet& Stream::recv() { // Acquire the receiver variables std::unique_lock rlck(recvMtx); // Wait until there are items that are readable or the receiver is notified that it should stop recvCV.wait(rlck, [=](){ return available || stopRecv; }); - // Return the number of readable items or -1 if the receiver should stop - return stopRecv ? -1 : available; + // Reset the available flag + available = false; + + // Zero out the number of samples if asked to stop + if (stopRecv) { recvSet->samples = 0; } + + // Return the buffer set + return *recvSet; } template - void Stream::ack() { + void Stream::flush() { // Acquire the sender variables std::unique_lock slck(sendMtx); @@ -139,5 +188,4 @@ namespace dsp { template class Stream; template class Stream; template class Stream; - template class Stream; } \ No newline at end of file diff --git a/dsp/stream.h b/dsp/stream.h index cd68c76..a7bfd95 100644 --- a/dsp/stream.h +++ b/dsp/stream.h @@ -1,49 +1,47 @@ #pragma once #include #include - -#define DSP_DEFAULT_BUFFER_SIZE 1000000 +#include "block.h" namespace dsp { /** - * Represents a class that can signal to its acting threads to stop. (TODO: Better name) + * Set of buffers and associated metadata. */ - class Signaler { - public: + template + struct BufferSet { /** - * Notify the sending thread that it should stop. + * Sample buffer for each channel. */ - virtual void stopSender() = 0; + T** buffer; /** - * Clear the sender stop flag to allow restarting the sender thread. + * Number of channels, and thus channel buffers. */ - virtual void clearSendStop() = 0; + size_t channels; /** - * Notify the receiving thread that it should stop. + * Maximum number of samples that each buffer can contain. This is assigned by Stream::reserve(). */ - virtual void stopReceiver() = 0; + size_t capacity; /** - * Clear the receiver stop flag to allow restarting the sender thread. + * Number of valid samples in each buffer. This is assigned by Stream::send(). */ - virtual void clearRecvStop() = 0; + size_t samples; }; /** * Streams allow to exchange samples between two threads. * The samples have to be of type (u)int8_t, (u)int16_t, (u)int32_t, (u)int64_t, float, double or Complex. + * This class is thread-safe. */ template - class Stream : public Signaler { + class Stream : public StopNotifier { public: - /** - * Create a stream object. - * @param bufferSize Number of items in the buffers. - */ - Stream(int channels = 1, int bufferSize = DSP_DEFAULT_BUFFER_SIZE); + // Default constructor + Stream(); + // Destructor ~Stream(); /** @@ -51,8 +49,6 @@ namespace dsp { */ void stopSender(); - // TODO: More consistent naming - /** * Clear the sender stop flag to allow restarting the sender thread. */ @@ -69,36 +65,31 @@ namespace dsp { void clearRecvStop(); /** - * Send a buffer of samples. - * @param count Number of samples in the send buffer. - * @return True if the sender thread must exist, false otherwise. + * Obtain a buffer set for sending. + * @param bufferSize Number of samples in each channel buffer. + * @param channels Number of channels. + * @return Buffer set to use for sending. */ - bool send(int count); + const BufferSet& reserve(size_t bufferSize, size_t channels = 1); /** - * Wait for buffer of samples. May also return in case of a signal to exit. ack() must be called as soon as the receive buffer has been entirely processed. - * @return Number of samples or -1 if the worker thread must exit. + * Send a set of sample buffers. + * @param count Number of valid samples in each channel buffer. + * @param channels Number of valid channels channels. + * @return False if the sender thread must exist, true otherwise. */ - int recv(); + bool send(size_t count); /** - * Acknowledge reception and processing of the samples. Allows sender thread to send a new buffer. + * Receive a set of sample buffers. May also return in case of a signal to exit, in which case the number of samples in the set is zero. + * @return Set of sample buffers. */ - void ack(); + const BufferSet& recv(); /** - * Get the sending buffer. - * @param channel ID of the channel to get the buffer for. - * @return Sending buffer. + * Flush the received buffer. Allows sender thread to send a new buffer. */ - T* getSendBuffer(int channel = 0); - - /** - * Get the receiving buffer. - * @param channel ID of the channel to get the buffer for. - * @return Sending buffer. - */ - const T* getRecvBuffer(int channel = 0); + void flush(); private: // Sender variables @@ -106,13 +97,13 @@ namespace dsp { std::mutex sendMtx; bool canSend = true; bool stopSend = false; - T* sendBuf; + BufferSet* sendSet = NULL; // Receiver variables std::condition_variable recvCV; std::mutex recvMtx; - int available = 0; + bool available = false; bool stopRecv = false; - T* recvBuf; + BufferSet* recvSet = NULL; }; } \ No newline at end of file diff --git a/dsp/taps.cpp b/dsp/taps.cpp index fbf615c..82bfdfb 100644 --- a/dsp/taps.cpp +++ b/dsp/taps.cpp @@ -12,7 +12,7 @@ namespace dsp { reallocate(count); // Null out if requested - if (zero) { memset(data, 0, count*sizeof(T)); } + if (zero) { memset(buffer, 0, count*sizeof(T)); } } template @@ -21,7 +21,7 @@ namespace dsp { reallocate(count); // Copy data over - memcpy(data, taps, count*sizeof(T)); + memcpy(buffer, taps, count*sizeof(T)); } template @@ -30,24 +30,24 @@ namespace dsp { reallocate(b.count); // Copy data over - memcpy(data, b.data, b.count*sizeof(T)); + memcpy(buffer, b.buffer, b.count*sizeof(T)); } template Taps::Taps(Taps&& b) { // Copy members - data = b.data; + buffer = b.buffer; count = b.count; // Neutralize old instance - b.data = NULL; + b.buffer = NULL; b.count = 0; } template Taps::~Taps() { // Free the buffer if it is allocated - if (data) { delete[] data; } + if (buffer) { delete[] buffer; } } template @@ -56,7 +56,7 @@ namespace dsp { reallocate(b.count); // Copy data over - memcpy(data, b.data, b.count*sizeof(T)); + memcpy(buffer, b.buffer, b.count*sizeof(T)); // Return self return *this; @@ -64,12 +64,15 @@ namespace dsp { template Taps& Taps::operator=(Taps&& b) { + // Destroy current instance + if (buffer) { delete[] buffer; } + // Copy members - data = b.data; + buffer = b.buffer; count = b.count; // Neutralize old instance - b.data = NULL; + b.buffer = NULL; b.count = 0; // Return self @@ -79,13 +82,14 @@ namespace dsp { template void Taps::reallocate(int count) { // If the new count is no different and the buffer is allocated, no need to realloc - if (data && this->count == count) { return; } + if (buffer && this->count == count) { return; } // Free buffer - if (data) { delete[] data; } + if (buffer) { delete[] buffer; } // Allocate buffer - data = new T[count, sizeof(T)]; + // TODO: Use volk instead + buffer = new T[count, sizeof(T)]; // Update tap count this->count = count; diff --git a/dsp/taps.h b/dsp/taps.h index df47e4d..af4c53a 100644 --- a/dsp/taps.h +++ b/dsp/taps.h @@ -3,10 +3,12 @@ namespace dsp { /** * Filter tap container. + * This class is NOT thread-safe. */ template class Taps { public: + // Default constructor Taps(); /** @@ -23,11 +25,19 @@ namespace dsp { */ Taps(T* taps, int count); + // Copy constructor Taps(const Taps& b); - Taps(Taps&& b); - ~Taps(); + // Move constructor + Taps(Taps&& b); + + // Destructor + virtual ~Taps(); + + // Copy assignment operator Taps& operator=(const Taps& b); + + // Move assignment operator Taps& operator=(Taps&& b); /** @@ -35,21 +45,25 @@ namespace dsp { */ inline int size() const { return count; } - inline operator const T*() const { return data; } + /** + * Get a pointer to the tap buffer. + */ + inline const T* data() const { return buffer; } - inline operator bool() const { return data && count; } + // Cast to bool operator + inline operator bool() const { return count; } /** * Get a tap by index. * @param index Index of the tap * @return Tap at index. */ - inline const T& operator[](int index) const { return data[index]; } + inline const T& operator[](int index) const { return buffer[index]; } protected: void reallocate(int count); int count = 0; - T* data = nullptr; + T* buffer = nullptr; }; } \ No newline at end of file diff --git a/dsp/taps/low_pass.h b/dsp/taps/low_pass.h index 1757f08..453b18b 100644 --- a/dsp/taps/low_pass.h +++ b/dsp/taps/low_pass.h @@ -2,11 +2,15 @@ #include "../taps.h" namespace dsp::taps { + // TODO: Add option to only use an odd number of taps and an option to query the delay + /** * Low-pass filter taps. + * This class is NOT thread-safe. */ class LowPass : public Taps { public: + // Default constructor LowPass(); /** @@ -54,8 +58,8 @@ namespace dsp::taps { private: void generate(); - float cutoff; - float transWidth; - float samplerate; + float cutoff = 0.0f; + float transWidth = 0.0f; + float samplerate = 0.0f; }; } \ No newline at end of file diff --git a/dsp/window.cpp b/dsp/window.cpp new file mode 100644 index 0000000..2237763 --- /dev/null +++ b/dsp/window.cpp @@ -0,0 +1,16 @@ +#include "window.h" + +namespace dsp { + Window::~Window() {} + + void Window::generate(float* data, size_t len) const { + for (size_t i = 0; i < len; i++) { + + } + } + + template + void Window::apply(T* data, size_t len) const { + // TODO + } +} \ No newline at end of file diff --git a/dsp/window.h b/dsp/window.h new file mode 100644 index 0000000..564cfa9 --- /dev/null +++ b/dsp/window.h @@ -0,0 +1,34 @@ +#pragma once +#include + +// TODO: Make something better, this sucks to use + +namespace dsp { + class Window { + public: + // Virtual destructor + virtual ~Window(); + + /** + * Compute the value of the window function. + * @param x Point at which to compute the window at. Must be bounded between 0 and 1. + * @return Value of the window bounded between 0 and 1. + */ + virtual float operator()(float x) const = 0; + + /** + * Generate a window of a given length. + * @param data Samples buffer to write the window to. + * @param len Length of the window. + */ + void generate(float* data, size_t len) const; + + /** + * Apply the window function to samples. + * @param data Sample buffer to apply the window to. + * @param len Length of the sample buffer. + */ + template + void apply(T* data, size_t len) const; + }; +} \ No newline at end of file diff --git a/dsp/window/rectangular.h b/dsp/window/rectangular.h new file mode 100644 index 0000000..ba71766 --- /dev/null +++ b/dsp/window/rectangular.h @@ -0,0 +1,14 @@ +#pragma once +#include "../window.h" + +namespace dsp::window { + class Rectangular : public Window { + public: + /** + * Compute the value of the window function. + * @param x Point at which to compute the window at. Must be bounded between 0 and 1. + * @return Value of the window bounded between 0 and 1. + */ + inline float operator()(float x) const { return 1.0f; } + }; +} \ No newline at end of file diff --git a/license b/license new file mode 100644 index 0000000..cf9d21f --- /dev/null +++ b/license @@ -0,0 +1 @@ +Copyright Alexandre Rouma 2025, All rights reserved. \ No newline at end of file diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..eb5aad2 --- /dev/null +++ b/main.cpp @@ -0,0 +1,48 @@ +#include +// #include "dsp/stream.h" +// #include "dsp/complex.h" +#include "dsp/window.h" +#include "dsp/window/rectangular.h" + +// void sendWorker(dsp::Stream& stream) { +// while (true) { +// // Obtrain a stereo buffer set +// auto set = stream.reserve(240, 2); + +// // Fill the buffer +// for (int i = 0; i < set.capacity; i++) { +// set.buffer[0][i] = rand(); +// set.buffer[1][i] = rand(); +// } + +// // Send the buffer +// if (!stream.send(set.capacity)) { break; } +// } +// } + +// void recvWorker(dsp::Stream& stream) { +// while (true) { +// // Receive a buffer set +// auto set = stream.recv(); +// if (!set.samples) { break; } + +// // Process +// // TODO: Do something + +// // Flush the buffer set +// stream.flush(); +// } +// } + +void test(const dsp::Window& win) { + float test[7]; + win.generate(test, 7); + for (int i = 0; i < 7; i++) { + printf("window[%d] = %f\n", i, test[i]); + } +} + +int main() { + test(dsp::window::Rectangular()); + return 0; +} \ No newline at end of file diff --git a/nuttall.txt b/nuttall.txt new file mode 100644 index 0000000..23bdd40 --- /dev/null +++ b/nuttall.txt @@ -0,0 +1,3 @@ +-1.000000000000000000e+02 -9.700000000000000000e+01 -9.400000000000000000e+01 -9.100000000000000000e+01 -8.800000000000000000e+01 -8.500000000000000000e+01 -8.200000000000000000e+01 -7.900000000000000000e+01 -7.600000000000000000e+01 -7.300000000000000000e+01 -7.000000000000000000e+01 -6.700000000000000000e+01 -6.400000000000000000e+01 -6.100000000000000000e+01 -5.800000000000000000e+01 -5.500000000000000000e+01 -5.200000000000000000e+01 -4.900000000000000000e+01 -4.600000000000000000e+01 -4.300000000000000000e+01 -4.000000000000000000e+01 -3.700000000000000000e+01 -3.400000000000000000e+01 -3.100000000000000000e+01 -2.800000000000000000e+01 -2.500000000000000000e+01 -2.200000000000000000e+01 -1.900000000000000000e+01 -1.600000000000000000e+01 -1.300000000000000000e+01 -1.000000000000000000e+01 +7.498058144869045982e+00 7.429123315241205283e+00 7.355936176181915087e+00 7.277968040549898987e+00 7.194763779640090284e+00 7.105923448779879692e+00 7.011078077859066227e+00 6.909876701203756966e+00 6.801969049096642017e+00 6.687001930717877407e+00 6.564605325502163247e+00 6.434386201983192777e+00 6.295917325052794666e+00 6.148729084702260650e+00 5.992296590643461762e+00 5.826027579346864549e+00 5.649247385514478914e+00 5.461178151598577557e+00 5.260915691572175312e+00 5.047397261226342025e+00 4.819359280697076642e+00 4.575281032595127861e+00 4.313304619145223562e+00 4.031122883301486937e+00 3.725812809806068771e+00 3.393582661555043956e+00 3.029366598614969597e+00 2.626145233288947889e+00 2.173720548942118302e+00 1.656301958095067084e+00 1.047129624943518911e+00 +-1.210822659459354744e+00 -1.175814880761542947e+00 -1.144720301542903274e+00 -1.117665666877913244e+00 -1.094173482616129345e+00 -1.074021147811728305e+00 -1.056953944892579100e+00 -1.042717251608808127e+00 -1.030994512107588079e+00 -1.021508038508901484e+00 -1.013982752724297276e+00 -1.008168659036855708e+00 -1.003789637790280409e+00 -1.000637525267440431e+00 -9.985171718632918081e-01 -9.971756639112729914e-01 -9.965044615951818008e-01 -9.963056341936008531e-01 -9.964478867235371240e-01 -9.968157311727851022e-01 -9.972917017908552451e-01 -9.978346445745125415e-01 -9.983239887393334788e-01 -9.987150255387784448e-01 -9.989609502361606053e-01 -9.990367948406020382e-01 -9.988947527069442778e-01 -9.984655720413724289e-01 -9.977159840282678882e-01 -9.964448203599918230e-01 -9.937137470899454206e-01 diff --git a/src/main.cpp b/src/main.cpp index f113d4a..734dde6 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,20 +1,6 @@ -#include "dsp/taps.h" -#include "dsp/taps/low_pass.h" -#include "dsp/complex.h" -#include "dsp/stream.h" #include -struct TestStruct { - bool a; - int b; -}; - int main() { - dsp::taps::LowPass lp(3000, 100, 15000); - - float test = lp[0]; - - dsp::Stream str; return 0; } \ No newline at end of file diff --git a/todo.txt b/todo.txt new file mode 100644 index 0000000..824563e --- /dev/null +++ b/todo.txt @@ -0,0 +1,9 @@ +Create buffer class that replaces the BufferSet struct + +Move the fixes from the SFAP library's stream.c to this one + +Think a lot about how multichannel streams are to be handled. Doing it all in one thread might not be possible so perhapse a built-in synchronisation system is worth it? + +Think about how blocks handle multi-channel streams. Maybe sometimes you want different behavior for each channel, so duplicating the DSP as is could be stupid + +Using multiple threads for audio may be slower than interleaved channels. It could also add latency since it has to be re-interleaved at the end of the DSP. \ No newline at end of file diff --git a/tools/filt_spec_reg.py b/tools/filt_spec_reg.py new file mode 100644 index 0000000..7eb4266 --- /dev/null +++ b/tools/filt_spec_reg.py @@ -0,0 +1,191 @@ +import numpy as np +import scipy.fft as fft +import scipy.interpolate as spi +import matplotlib.pyplot as plt +import scipy.optimize as opt + +def sincFIR(b, N, win): + lim = b * (N-1) * 0.5 + x = np.linspace(-lim, lim, N) + taps = np.sinc(x) * win(N) + taps *= 1 / np.sum(taps) + return taps + +def response(taps, points = 65536): + # Zero pad + ztaps = np.concatenate([taps, np.zeros(points - len(taps))]) + + # Compute FFT + ff = fft.fft(ztaps) + ff = ff[:len(ff)//2] + + # Compite the frequency of each bin + freqs = np.array(range(points))/points + freqs = freqs[:len(freqs)//2]*2 + + # Return the amplitude and phase + return np.array(freqs), 20*np.log10(np.abs(ff)), np.angle(ff)*180/np.pi + +def getSpecs(taps, att = None): + # Get the response + freqs,amp,_ = response(taps) + + # Find the inflection point + searchIdx = -1 + mid = len(amp)//2 + last = amp[mid-1] + for i in range(mid, len(amp)): + diff = amp[i]-last + last = amp[i] + if diff > 0.0: + searchIdx = i + break + searchFreq = freqs[searchIdx-1] + + # Find the stop band attenuation if none given + if att == None: + maxIdx = np.argmax(amp[searchIdx:]) + searchIdx + att = amp[maxIdx] + + # Find the transition point + interp = spi.interp1d(freqs, amp) + res = opt.root_scalar(lambda x : interp(x)-att, bracket=[0.5, searchFreq]) + return att, res.root - 0.5 + +def bnuttallw(n, N): + return 0.3635819 - 0.4891775*np.cos(2*np.pi*n/N) + 0.13659955*np.cos(4*np.pi*n/N) - 0.0106411*np.cos(6*np.pi*n/N) + +# def lolw(n, N): +# return 0.35881115 - 0.48845026*np.cos(2*np.pi*n/N) + 0.14118885*np.cos(4*np.pi*n/N) - 0.01154974*np.cos(6*np.pi*n/N) + +def nuttallw(n, N): + return 0.355768 - 0.487396*np.cos(2*np.pi*n/N) + 0.144232*np.cos(4*np.pi*n/N) - 0.012604*np.cos(6*np.pi*n/N) + +def lol2w(n, N): + return 0.3560283 - 0.48748693*np.cos(2*np.pi*n/N) + 0.1439717*np.cos(4*np.pi*n/N) - 0.01251307*np.cos(6*np.pi*n/N) + +def nuttall(N): + return bnuttallw(np.array(range(N+1)), N+1)[1:] + +# def lol(N): +# return lolw(np.array(range(N+1)), N+1)[1:] + +def lol2(N): + return lol2w(np.array(range(N+1)), N+1)[1:] + +def nuttallT(N): + return nuttallw(np.array(range(N)), N-1) + +# Previous prediction 1824 +# nuttall = [ 7.74932575 -0.99720463] 1859 +# nuttallT = [ 7.749348 1.00079625] 1861 +# np.blackman = [5.58715521 0.98759599] 1341 +# np.hanning = [3.12375105 1.02086163] 751 +# bnuttal = [ 7.79530085 -2.42566429] 1867 + +def charactWindow(win, minTaps = 32, maxTaps = 1024): + ii = [] + att = [] + tw = [] + for i in range(minTaps, maxTaps+1): + taps = sincFIR(0.5, i, win) + a,t = getSpecs(taps) + ii.append(i) + att.append(a) + tw.append(t) + res = np.polyfit(1/np.array(tw), ii, 1) + return res[0], res[1], np.mean(att) + +# print(charactWindow(lol2)) +# exit(0) + +freqs, amps, phases = response(sincFIR(0.5, 101, lol2)) +plt.plot(freqs, amps) +freqs, amps, phases = response(sincFIR(0.5, 101, nuttall)) +plt.plot(freqs, amps) +plt.show() + +# params = [] +# alphas = [] +# atts = [] +# for i in range(1000): +# p = np.random.uniform(0.0, 1.0, 3) +# p = np.concatenate([ [p[0] - p[1] + p[2]], p ]) + +# def winw(n, N): +# return p[0] - p[1]*np.cos(2*np.pi*n/N) + p[2]*np.cos(4*np.pi*n/N) - p[3]*np.cos(6*np.pi*n/N) + +# def win(N): +# return winw(np.array(range(N+1)), N+1)[1:] + +# try: +# alpha, beta, att = charactWindow(win, 1000, 1024) +# params.append(p) +# alphas.append(alpha) +# atts.append(-att) +# print(p, alpha, beta, att) +# except: +# pass + + +# best = np.argmax(atts) +# print(alphas[best], atts[best], params[best]) + +# plt.scatter(atts, alphas) +# plt.show() + +# (5.587154476193536, 0.9876436878253441, -75.28786662414193) +# (5.587618976822818, 0.9168073377669025, -75.28843748397959) + +# 7.4081533102937565 75.34448596788963 [0.35417698 0.95739698 0.78534742 0.18212743] + +p = [0.35417698, 0.95739698, 0.78534742, 0.18212743] + +def customWin(N, p): + def winw(n, N): + return p[0] - p[1]*np.cos(2*np.pi*n/N) + p[2]*np.cos(4*np.pi*n/N) - p[3]*np.cos(6*np.pi*n/N) + return winw(np.array(range(N+1)), N+1)[1:] + +def optfunc(x): + x = [ np.power(10, x[0]), np.power(10, x[1]), np.power(10, x[2]) ] + p = np.concatenate([ [x[0] - x[1] + x[2]], x ]) + test = customWin(301, p) + test *= 1/np.sum(test) + + p *= 1/np.sum(p) + print(p) + plt.plot(test) + plt.show() + + freqs, amps, _ = response(test) + + # Find the inflection point + searchIdx = -1 + last = amps[0] + for i in range(1, len(amps)): + diff = amps[i]-last + last = amps[i] + if diff > 0.0: + searchIdx = i + break + + # Find stop band attenuation + att = np.max(amps[searchIdx:]) + + # Optimise for attenuation + return att + +optfunc([-0.57685987, -1.10654573, -2.16745906]) + +# best = 0 +# bestCoeffs = None +# for i in range(1000): +# res = opt.minimize(optfunc, np.random.uniform(-10, 0, 3)) +# if res.fun < best: +# best = res.fun +# bestCoeffs = res.x +# print('Result:', i, res.fun, best) +# print('Best:', best, bestCoeffs) + +# # Cool window: [0.38373144 0.49571507 0.11626856 0.00428493] + diff --git a/tools/spec_fit.py b/tools/spec_fit.py new file mode 100644 index 0000000..b1cc88e --- /dev/null +++ b/tools/spec_fit.py @@ -0,0 +1,33 @@ +import numpy as np +import scipy.fft as fft +import scipy.interpolate as spi +import matplotlib.pyplot as plt +import scipy.optimize as opt + +data = np.loadtxt('nuttall.txt') + +atts = data[0] +alphas = data[1] +betas = data[2] + +# def sqerr(xs, ys, f): +# err = 0 +# for i in range(len(xs)): +# x = xs[i] +# dy = ys[i] +# fy = f(x) +# err += (fy - dy)**2 +# return err + +# def compErr(d): +# p = np.polyfit(atts, alphas, d) +# return sqerr(atts, alphas, lambda x : np.polyval(p, x)) + +# p = np.polyfit(atts, alphas, 2) + +# plt.plot(atts, alphas) +# x = np.linspace(atts[0], atts[-1], 1000) +# plt.plot(x, np.polyval(p, x)) +# plt.show() +plt.plot(atts, alphas) +plt.show() \ No newline at end of file