progress

.gitignore

@@ -1,2 +1,3 @@
 .vscode/
 build/
+.old/

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)
+
+
+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 ()

dsp/block.cpp

@@ -1,8 +1,12 @@
 #include "block.h"
+#include <stdexcept>
 
 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<std::mutex> 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<std::mutex> 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<std::mutex> 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<std::mutex> 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<std::mutex> 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());
     }
 }

dsp/block.h

@@ -1,16 +1,46 @@
 #pragma once
 #include <thread>
+#include <mutex>
 #include <vector>
-#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.
+         * Register an input stop notifier.
-         * @param input Input stream to register.
+         * @param input Input stop notifier to register.
         */
-        void registerInput(Signaler* input);
+        void registerInput(StopNotifier* input);
 
         /**
-         * Unregister an input stream.
+         * Unregister an input stop notifier.
-         * @param input Input stream to unregister.
+         * @param input Input stop notifier to unregister.
         */
-        void unregisterInput(Signaler* input);
+        void unregisterInput(StopNotifier* input);
 
         /**
-         * Register an output stream.
+         * Register an output stop notifier.
-         * @param input Output stream to register.
+         * @param input Output stop notifier to register.
         */
-        void registerOutput(Signaler* output);
+        void registerOutput(StopNotifier* output);
 
         /**
-         * Unregister an output stream.
+         * Unregister an output stop notifier.
-         * @param input Output stream to unregister.
+         * @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();
+        virtual bool run() = 0;
-
-        /**
-         * Mutex to be used for block settings.
-        */
-        std::recursive_mutex settingsMtx;
 
     private:
         /**
@@ -71,10 +98,10 @@ namespace dsp {
         void worker();
 
         // Worker variables
-        std::mutex workerMtx;
+        mutable std::mutex workerMtx;
         std::thread workerThread;
-        std::vector<Signaler*> inputs;
+        std::vector<StopNotifier*> inputs;
-        std::vector<Signaler*> outputs;
+        std::vector<StopNotifier*> outputs;
         bool _running = false;
     };
 }

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 };
 }

dsp/demod/fm.cpp

@@ -0,0 +1,94 @@
+#include "fm.h"
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+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<std::recursive_mutex> lck(settingsMtx);
+
+        // Return the deviation
+        return deviation;
+    }
+
+    void FM_s::setDeviation(float deviation) {
+        // Acquire the settings mutex
+        std::lock_guard<std::recursive_mutex> 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<std::recursive_mutex> lck(settingsMtx);
+
+        // Return the samplerate
+        return samplerate;
+    }
+
+    void FM_s::setSamplerate(float samplerate) {
+        // Acquire the settings mutex
+        std::lock_guard<std::recursive_mutex> 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<std::recursive_mutex> 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<Complex>* in, float deviation, float samplerate) :
+        FM_s(deviation, samplerate),
+        ProcessorAsync(this, in)
+    {}
+
+    FM::FM(Stream<Complex>* in, Stream<float>* out, float deviation, float samplerate) :
+        FM_s(deviation, samplerate),
+        ProcessorAsync(this, in, out)
+    {}
+}

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<Complex, float> {
+    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<Complex, float> {
+    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<Complex>* 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<Complex>* in, Stream<float>* out, float deviation, float samplerate);
+    };
+}

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 SAMP_T, class TAPS_T>
+    class FIR_s : public ProcessorSync<SAMP_T, SAMP_T> {
+    public:
+        /**
+         * Create a FIR filter.
+         * @param taps Filter taps.
+        */
+        FIR_s(const Taps<TAPS_T>& taps);
+
+        /**
+         * Get the filter taps.
+         * @return Filter taps.
+        */
+        Taps<TAPS_T> getTaps();
+
+        /**
+         * Set the filter taps.
+         * @param taps Filter taps.
+        */
+        void setTaps(const Taps<TAPS_T>& 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_T> taps;
+    };
+
+    /**
+     * Finite Inpulse Response filter.
+     * This class is thread-safe except for `process()`.
+    */
+    template <class SAMP_T, class TAPS_T>
+    class FIR : public FIR_s<SAMP_T, TAPS_T>, public ProcessorAsync<SAMP_T, SAMP_T> {
+    public:
+        /**
+         * Create a FIR filter.
+         * @param taps Filter taps.
+        */
+        FIR(const Taps<TAPS_T>& taps);
+
+        /**
+         * Create a FIR filter.
+         * @param in Input samples.
+         * @param taps Filter taps.
+        */
+        FIR(Stream<SAMP_T>* in, const Taps<TAPS_T>& taps);
+
+        /**
+         * Create a FIR filter.
+         * @param in Input samples.
+         * @param in Filtered samples.
+         * @param taps Filter taps.
+        */
+        FIR(Stream<SAMP_T>* in, Stream<SAMP_T>* out, const Taps<TAPS_T>& taps);
+    };
+}

dsp/mailbox.h

@@ -1,72 +0,0 @@
-#pragma once
-#include "stream.h"
-
-namespace gui {
-    /**
-     * Mailboxes allow to exchange objects between two threads.
-    */
-    template <typename T>
-    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;
-    };
-}

dsp/processor.h

@@ -1,12 +1,164 @@
 #pragma once
+#include <type_traits>
 #include "block.h"
+#include "stream.h"
 
 namespace dsp {
-    template <typename I, typename O>
+    // TODO: Deal with the fact that always locking will be slow in hier blocks...
-    class Processor : public Block {
+
+    /**
+     * 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 IN, class OUT>
+    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<std::recursive_mutex> getLock() const {
+            return std::lock_guard<std::recursive_mutex>(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 IN, class OUT>
+    class ProcessorBlock : public Block {
+    public:
+        // Default constructor
+        ProcessorBlock() {
+            // TODO: Maybe something to do to prevent bad shit if started?
+        }
+
+        /**
+         * TODO
+        */
+        ProcessorBlock(ProcessorSync<IN, OUT>* proc) {
+            // Save the pointer to the processor
+            this->proc = proc;
+
+            // Set the streams
+            setInput(NULL);
+            setOutput(NULL);
+        }
+
+        /**
+         * TODO
+        */
+        ProcessorBlock(ProcessorSync<IN, OUT>* proc, Stream<IN>* in) {
+            // Save the pointer to the processor
+            this->proc = proc;
+
+            // Set the streams
+            setInput(in);
+            setOutput(NULL);
+        }
+
+        /**
+         * TODO
+        */
+        ProcessorBlock(ProcessorSync<IN, OUT>* proc, Stream<IN>* in, Stream<OUT>* 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>* 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>* 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>* in() const {
+            // TODO: Lock
+            return _in;
+        }
+
+        Stream<OUT>* 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<IN, OUT>* proc;
+        Stream<IN>* _in = NULL;
+        Stream<OUT>* _out = NULL;
+        bool ownInput = false;
+        bool ownOutput = false;
     };
 }

dsp/stereo.h

@@ -1,90 +0,0 @@
-#pragma once
-#include <math.h>
-
-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 };
-    }
-}

dsp/stream.cpp

@@ -1,20 +1,29 @@
-#include "Stream.h"
+#include "stream.h"
 #include "./complex.h"
-#include "stereo.h"
 
 namespace dsp {
     template <typename T>
-    Stream<T>::Stream(int channels, int bufferSize) {
+    Stream<T>::Stream() {
-        // Allocate both send and receive buffers aligned by the size of the type
+        // Allocate both send and receive buffer sets
-        sendBuf = new T[bufferSize, sizeof(T)];
+        sendSet = new BufferSet<T>;
-        recvBuf = new T[bufferSize, sizeof(T)];
+        recvSet = new BufferSet<T>;
+
+        // 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 <typename T>
     Stream<T>::~Stream() {
-        // Free both send and receive buffers
+        // Free both send and receive buffer sets
-        delete[] sendBuf;
+        delete sendSet;
-        delete[] recvBuf;
+        delete recvSet;
     }
 
     template <typename T>
@@ -66,15 +75,46 @@ namespace dsp {
     }
 
     template <typename T>
-    bool Stream<T>::send(int count) {
+    const BufferSet<T>& Stream<T>::reserve(size_t bufferSize, size_t channels) {
         // Acquire the sender variables
         std::unique_lock<std::mutex> 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 <typename T>
+    bool Stream<T>::send(size_t count) {
+        // Acquire the sender variables
+        std::unique_lock<std::mutex> 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<std::mutex> rlck(recvMtx);
 
-        // Swap buffers
+        // Swap the buffers sets
-        T* tmp = sendBuf;
+        BufferSet<T>* tmp = sendSet;
-        sendBuf = recvBuf;
+        sendSet = recvSet;
-        recvBuf = tmp;
+        recvSet = tmp;
 
         // Release the sender variables
         slck.unlock();
 
-        // Set the number of items that are readable
+        // Set the available flag
-        available = count;
+        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 <typename T>
-    int Stream<T>::recv() {
+    const BufferSet<T>& Stream<T>::recv() {
         // Acquire the receiver variables
         std::unique_lock<std::mutex> 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
+        // Reset the available flag
-        return  stopRecv ? -1 : available;
+        available = false;
+
+        // Zero out the number of samples if asked to stop
+        if (stopRecv) { recvSet->samples = 0; }
+
+        // Return the buffer set
+        return  *recvSet;
     }
 
     template <typename T>
-    void Stream<T>::ack() {
+    void Stream<T>::flush() {
         // Acquire the sender variables
         std::unique_lock<std::mutex> slck(sendMtx);
 
@@ -139,5 +188,4 @@ namespace dsp {
     template class Stream<float>;
     template class Stream<double>;
     template class Stream<Complex>;
-    template class Stream<Stereo>;
 }

dsp/stream.h

@@ -1,49 +1,47 @@
 #pragma once
 #include <mutex>
 #include <condition_variable>
-
+#include "block.h"
-#define DSP_DEFAULT_BUFFER_SIZE 1000000
 
 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 {
+    template <typename T>
-    public:
+    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<T>::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<T>::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 <typename T>
-    class Stream : public Signaler {
+    class Stream : public StopNotifier {
     public:
-        /**
+        // Default constructor
-         * Create a stream object.
+        Stream();
-         * @param bufferSize Number of items in the buffers.
-        */
-        Stream(int channels = 1, int bufferSize = DSP_DEFAULT_BUFFER_SIZE);
 
+        // 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.
+         * Obtain a buffer set for sending.
-         * @param count Number of samples in the send buffer.
+         * @param bufferSize Number of samples in each channel buffer.
-         * @return True if the sender thread must exist, false otherwise.
+         * @param channels Number of channels.
+         * @return Buffer set to use for sending.
         */
-        bool send(int count);
+        const BufferSet<T>& 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.
+         * Send a set of sample buffers.
-         * @return Number of samples or -1 if the worker thread must exit.
+         * @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<T>& recv();
 
         /**
-         * Get the sending buffer.
+         * Flush the received buffer. Allows sender thread to send a new buffer.
-         * @param channel ID of the channel to get the buffer for.
-         * @return Sending buffer.
         */
-        T* getSendBuffer(int channel = 0);
+        void flush();
-
-        /**
-         * Get the receiving buffer.
-         * @param channel ID of the channel to get the buffer for.
-         * @return Sending buffer.
-        */
-        const T* getRecvBuffer(int channel = 0);
 
     private:
         // Sender variables
@@ -106,13 +97,13 @@ namespace dsp {
         std::mutex sendMtx;
         bool canSend = true;
         bool stopSend = false;
-        T* sendBuf;
+        BufferSet<T>* sendSet = NULL;
 
         // Receiver variables
         std::condition_variable recvCV;
         std::mutex recvMtx;
-        int available = 0;
+        bool available = false;
         bool stopRecv = false;
-        T* recvBuf;
+        BufferSet<T>* recvSet = NULL;
     };
 }

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 <class T>
@@ -21,7 +21,7 @@ namespace dsp {
         reallocate(count);
 
         // Copy data over
-        memcpy(data, taps, count*sizeof(T));
+        memcpy(buffer, taps, count*sizeof(T));
     }
 
     template <class T>
@@ -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 <class T>
     Taps<T>::Taps(Taps<T>&& 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 <class T>
     Taps<T>::~Taps() {
         // Free the buffer if it is allocated
-        if (data) { delete[] data; }
+        if (buffer) { delete[] buffer; }
     }
 
     template <class T>
@@ -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 <class T>
     Taps<T>& Taps<T>::operator=(Taps<T>&& 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 <class T>
     void Taps<T>::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;

dsp/taps.h

@@ -3,10 +3,12 @@
 namespace dsp {
     /**
      * Filter tap container.
+     * This class is NOT thread-safe.
     */
     template <class T>
     class Taps {
     public:
+        // Default constructor
         Taps();
         
         /**
@@ -23,11 +25,19 @@ namespace dsp {
         */
         Taps(T* taps, int count);
 
+        // Copy constructor
         Taps(const Taps<T>& b);
-        Taps(Taps<T>&& b);
-        ~Taps();
 
+        // Move constructor
+        Taps(Taps<T>&& b);
+
+        // Destructor
+        virtual ~Taps();
+
+        // Copy assignment operator
         Taps<T>& operator=(const Taps<T>& b);
+
+        // Move assignment operator
         Taps<T>& operator=(Taps<T>&& 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;
     };
 }

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<float> {
     public:
+        // Default constructor
         LowPass();
 
         /**
@@ -54,8 +58,8 @@ namespace dsp::taps {
     private:
         void generate();
 
-        float cutoff;
+        float cutoff = 0.0f;
-        float transWidth;
+        float transWidth = 0.0f;
-        float samplerate;
+        float samplerate = 0.0f;
     };
 }

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 <class T>
+    void Window::apply(T* data, size_t len) const {
+        // TODO
+    }
+}

dsp/window.h

@@ -0,0 +1,34 @@
+#pragma once
+#include <stdint.h>
+
+// 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 <class T>
+        void apply(T* data, size_t len) const;
+    };
+}

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; }
+    };
+}

license

@@ -0,0 +1 @@
+Copyright Alexandre Rouma 2025, All rights reserved.

main.cpp

@@ -0,0 +1,48 @@
+#include <stdio.h>
+// #include "dsp/stream.h"
+// #include "dsp/complex.h"
+#include "dsp/window.h"
+#include "dsp/window/rectangular.h"
+
+// void sendWorker(dsp::Stream<float>& 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<float>& 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;
+}

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

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 <stdio.h>
 
-struct TestStruct {
-    bool a;
-    int b;
-};
-
 int main() {
-    dsp::taps::LowPass lp(3000, 100, 15000);
-
-    float test = lp[0];
-
-    dsp::Stream<float> str;
 
     return 0;
 }

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.

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]
+

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()