Added very basic peak finder - not very good, just an example for now

Author: thomasms

examples/realspectrum.py

@@ -11,22 +11,27 @@
 pkd.io.from_csv(hist_raw, filename)
 energies, counts = hist_raw.X, hist_raw.Y
 
-# smooth it
-smoother = pkd.core.MovingAverageSmoother(3)
-smoother2 = pkd.core.SavitzkyGolaySmoother(3)
-counts1 = smoother.go(counts)
-hist_smoothed1 = pkd.core.SpectrumEnergyBased(energies, counts1)
-counts2 = smoother2.go(counts)
-hist_smoothed2 = pkd.core.SpectrumEnergyBased(energies, counts2)
-
 # estimate and then remove background
 ITERS = range(1, 21)
 # forward windowing
-background_forw = hist_smoothed1.estimateBackground(ITERS)
+background_forw = hist_raw.estimateBackground(ITERS)
 # backward windowing
-background_back = hist_smoothed1.estimateBackground(list(reversed(ITERS)))
-hist_smoothed1.removeBackground(ITERS)
-snippedCounts = hist_smoothed1.Y
+background_back = hist_raw.estimateBackground(list(reversed(ITERS)))
+hist_raw.removeBackground(ITERS)
+snippedCounts = hist_raw.Y
+
+# process manager
+pm = pkd.core.SimpleProcessManager()
+
+# smooth
+sg = pkd.core.SavitzkyGolaySmoother(31)
+pm.append(sg)
+
+# peak finding
+pm.append(pkd.core.MovingAveragePeakFinder(400))
+
+# process
+processed_counts = pm.run(snippedCounts)
 
 # function to plot histogram
 def getplotvalues(x,y):
@@ -47,11 +52,11 @@ def getplotvalues(x,y):
 import matplotlib.pyplot as plt
 
 f = plt.figure(figsize=(10,7))
-plt.semilogy(*getplotvalues(energies, counts1), 'k', linewidth=0.8, alpha=0.5, label="smoothed1")
-plt.semilogy(*getplotvalues(energies, counts2), 'k--', linewidth=0.8, alpha=1.0, label="smoothed2")
+plt.semilogy(*getplotvalues(energies, counts), 'k', linewidth=0.8, alpha=0.5, label="raw")
+plt.semilogy(*getplotvalues(energies, processed_counts), 'r', linewidth=2.0, alpha=0.8, label="processed")
 # plt.semilogy(*getplotvalues(energies, snippedCounts), 'r', linewidth=0.4, alpha=0.5, label="no background")
 plt.semilogy(*getplotvalues(energies, background_forw), 'g', linewidth=1.0, alpha=0.8, label="background - forwards")
-plt.semilogy(*getplotvalues(energies, background_back), 'r', linewidth=1.0, alpha=0.8, label="background - backwards")
+plt.semilogy(*getplotvalues(energies, background_back), 'c', linewidth=1.0, alpha=0.8, label="background - backwards")
 plt.ylim(bottom=0.5, top=1e4)
 plt.xlabel("Energy (eV)", fontsize=24)
 plt.ylabel("Counts", fontsize=24)

include/core.hpp

@@ -18,5 +18,6 @@
 #include "core/process.hpp"
 #include "core/smoothing.hpp"
 #include "core/spectral.hpp"
+#include "core/peaking.hpp"
 
 #endif //CORE_HPP

include/core/peaking.hpp

@@ -0,0 +1,56 @@
+////////////////////////////////////////////////////////////////////
+//                                                                //
+//    Copyright (c) 2019-20, UK Atomic Energy Authority (UKAEA)   //
+//                                                                //
+////////////////////////////////////////////////////////////////////
+
+/*!
+    @file
+    Defines peak finding functions to be applied to numerical array
+
+    @copyright UK Atomic Energy Authority (UKAEA) - 2019-20
+*/
+#ifndef CORE_PEAKING_HPP
+#define CORE_PEAKING_HPP
+
+#include <cmath>
+#include <memory>
+
+#include "common.hpp"
+#include "core/numerical.hpp"
+#include "core/process.hpp"
+#include "core/smoothing.hpp"
+
+PEAKINGDUCK_NAMESPACE_START(peakingduck)
+PEAKINGDUCK_NAMESPACE_START(core)
+
+    /*!
+       @brief Simple moving average peak finding 
+    */
+    template<typename T=DefaultType, int Size=ArrayTypeDynamic>
+    struct MovingAveragePeakFinder : public IProcess<T, Size>
+    {
+
+        explicit MovingAveragePeakFinder(int windowsize)
+        {
+            _movingAverageSmoother = std::make_shared<MovingAverageSmoother<T,Size>>(windowsize);
+        }
+
+        NumericalData<T, Size> 
+        go(const NumericalData<T, Size>& data) const override final
+        {
+            NumericalData<T, Size> smoothed = data - _movingAverageSmoother->go(data);
+            for(int i=0;i<data.size();++i)
+                smoothed[i] = smoothed[i] > 0 ? smoothed[i] : 0.0;
+            return smoothed;
+        };
+
+      private:
+        std::shared_ptr<IProcess<T,Size>> _movingAverageSmoother;
+    };  
+
+
+PEAKINGDUCK_NAMESPACE_END
+PEAKINGDUCK_NAMESPACE_END
+
+#endif // CORE_PEAKING_HPP

py/peakingduck.cpp

@@ -263,6 +263,11 @@ PYBIND11_MODULE(PEAKINGDUCK, m) {
     py::class_<WeightedMovingAverageSmootherPyType, IProcessPyType, std::shared_ptr<WeightedMovingAverageSmootherPyType>>(m_core, "WeightedMovingAverageSmoother")
         .def(py::init<int>());
 
+    // peak finding objects
+    using MovingAveragePeakFinderPyType = core::MovingAveragePeakFinder<NumericalDataCoreType,core::ArrayTypeDynamic>;
+    py::class_<MovingAveragePeakFinderPyType, IProcessPyType, std::shared_ptr<MovingAveragePeakFinderPyType>>(m_core, "MovingAveragePeakFinder")
+        .def(py::init<int>());
+
     // histogram objects
     using HistPyType = core::Histogram<double,double>;
     using HistChannelPyType = core::Histogram<int,double>;