diff --git a/src/config.rs b/src/config.rs
index fee5277..cf6b0d9 100644
--- a/src/config.rs
+++ b/src/config.rs
@@ -11,7 +11,7 @@ cfg_if::cfg_if! {
         /// Ratio between circumference and diameter of a circle
         pub const pi: Flt = std::f64::consts::PI;
 
-        
+
     }
     else if #[cfg(feature="f32")] {
         /// Floating-point value, compile time option to make it either f32, or f64
@@ -26,7 +26,8 @@ cfg_if::cfg_if! {
 
 cfg_if::cfg_if! {
 if #[cfg(feature = "python-bindings")] {
-    pub use numpy::{IntoPyArray,PyArray, PyArray1, PyArrayDyn, PyArrayLike1, PyReadonlyArrayDyn};
+    pub use numpy::{IntoPyArray,PyArray, PyArray1, PyArray2, PyArray3, PyArrayDyn, PyArrayLike1,
+        PyArrayLike2,PyArrayLike3,PyReadonlyArrayDyn, convert::ToPyArray};
     pub use pyo3::prelude::*;
     pub use pyo3::exceptions::PyValueError;
     pub use pyo3::{pymodule, types::PyModule, PyResult};
diff --git a/src/lib.rs b/src/lib.rs
index 6cd2602..8430507 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -64,6 +64,9 @@ fn lasprs(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_class::<slm::SLM>()?;
     m.add_class::<ps::WindowType>()?;
     m.add_class::<ps::Overlap>()?;
+    m.add_class::<ps::ApsMode>()?;
+    m.add_class::<ps::ApsSettings>()?;
+    m.add_class::<ps::AvPowerSpectra>()?;
 
     Ok(())
 }
diff --git a/src/ps/aps.rs b/src/ps/aps.rs
index 467d9f4..c32af7f 100644
--- a/src/ps/aps.rs
+++ b/src/ps/aps.rs
@@ -1,198 +1,25 @@
-use super::timebuffer::TimeBuffer;
-use super::CrossPowerSpecra;
 use super::*;
+use super::{timebuffer::TimeBuffer, CrossPowerSpecra};
 use crate::{config::*, TransferFunction, ZPKModel};
 use anyhow::{bail, Error, Result};
-use derive_builder::Builder;
 use freqweighting::FreqWeighting;
 
-/// All settings used for computing averaged power spectra using Welch' method.
-#[derive(Builder, Clone)]
-#[builder(build_fn(validate = "Self::validate", error = "Error"))]
-pub struct ApsSettings {
-    /// Mode of computation, see [ApsMode].
-    #[builder(default)]
-    mode: ApsMode,
-    /// Overlap in time segments. See [Overlap].
-    #[builder(default)]
-    overlap: Overlap,
-    /// Window applied to time segments. See [WindowType].
-    #[builder(default)]
-    windowType: WindowType,
-    /// Kind of freqency weighting. Defaults to Z
-    #[builder(default)]
-    freqWeightingType: FreqWeighting,
-    /// FFT Length
-    nfft: usize,
-    /// Sampling frequency
-    fs: Flt,
-}
-
-impl ApsSettingsBuilder {
-    fn validate(&self) -> Result<()> {
-        if !self.fs.is_some() {
-            bail!("Sampling frequency not given");
-        }
-        let fs = self.fs.unwrap();
-
-        if !fs.is_normal() {
-            bail!("Sampling frequency not a normal number")
-        }
-        if fs <= 0.0 {
-            bail!("Invalid sampling frequency given as parameter");
-        }
-
-        if self.nfft.is_none() {
-            bail!("nfft not specified")
-        };
-        let nfft = self.nfft.unwrap();
-        if nfft % 2 != 0 {
-            bail!("NFFT should be even")
-        }
-        if nfft == 0 {
-            bail!("Invalid NFFT, should be > 0.")
-        }
-        // Perform some checks on ApsMode
-        if let Some(ApsMode::ExponentialWeighting { tau }) = self.mode {
-            if tau <= 0.0 {
-                bail!("Invalid time weighting constant [s]. Should be > 0 if given.");
-            }
-        }
-
-        Ok(())
-    }
-}
-impl ApsSettings {
-    /// Returns nfft
-    pub fn nfft(&self) -> usize {
-        self.nfft
-    }
-    fn get_overlap_keep(&self) -> usize {
-        self.validate_get_overlap_keep().unwrap()
-    }
-    /// Returns the amount of samples to `keep` in the time buffer when
-    /// overlapping time segments using [TimeBuffer].
-    fn validate_get_overlap_keep(&self) -> Result<usize> {
-        let nfft = self.nfft;
-        let overlap_keep = match self.overlap {
-            Overlap::Number { N } if N >= nfft => {
-                bail!("Invalid overlap number of samples. Should be < nfft, which is {nfft}.")
-            }
-            // Keep 1 sample, if overlap is 1 sample etc.
-            Overlap::Number { N } if N < nfft => N,
-
-            // If overlap percentage is >= 100, or < 0.0 its an error
-            Overlap::Percentage { pct } if !(0.0..100.).contains(&pct) => {
-                bail!("Invalid overlap percentage. Should be >= 0. And < 100.")
-            }
-            // If overlap percentage is 0, this gives
-            Overlap::Percentage { pct } => ((pct * nfft as Flt) / 100.) as usize,
-            Overlap::NoOverlap {} => 0,
-            _ => unreachable!(),
-        };
-        if overlap_keep >= nfft {
-            bail!("Computed overlap results in invalid number of overlap samples. Please make sure the FFT length is large enough, when high overlap percentages are required.");
-        }
-        Ok(overlap_keep)
-    }
-
-    /// Return a reasonable acoustic default with a frequency resolution around
-    /// ~ 10 Hz, where nfft is still an integer power of 2.
-    ///
-    /// # Errors
-    ///
-    /// If `fs` is something odd, i.e. < 1 kHz, or higher than 1 MHz.
-    ///
-    pub fn reasonableAcousticDefault(fs: Flt, mode: ApsMode) -> Result<ApsSettings> {
-        if fs < 1e3 || fs > 1e6 {
-            bail!("Sampling frequency for reasonable acoustic data is >= 1 kHz and <= 1 MHz.");
-        }
-        let fs_div_10_rounded = (fs / 10.) as u32;
-
-        // 2^30 is about 1 million. We search for a two-power of an nfft that is
-        // the closest to fs/10. The frequency resolution is about fs/nfft.
-        let nfft = (0..30).map(|i| 2u32.pow(i) - fs_div_10_rounded).fold(
-            // Start wth a value that is always too large
-            fs as u32 * 10,
-            |cur, new| cur.min(new),
-        ) as usize;
-
-        Ok(ApsSettings {
-            mode,
-            fs,
-            nfft,
-            windowType: WindowType::default(),
-            overlap: Overlap::default(),
-            freqWeightingType: FreqWeighting::default(),
-        })
-    }
-
-    /// Return sampling frequency
-    pub fn fs(&self) -> Flt {
-        self.fs
-    }
-
-    /// Return Nyquist frequency
-    pub fn fnyq(&self) -> Flt {
-        self.fs() / 2.
-    }
-
-    /// Returns a single-sided frequency array corresponding to points in Power
-    /// spectra computation.
-    pub fn getFreq(&self) -> Array1<Flt> {
-        let df = self.fs / self.nfft as Flt;
-        let K = self.nfft / 2 + 1;
-        Array1::linspace(0., (K - 1) as Flt * df, K)
-    }
-}
-
-
-/// The 'mode' used in computing averaged power spectra. When providing data in
-/// blocks to the [AvPowerSpectra] the resulting 'current estimate' responds
-/// differently, depending on the model.
-#[derive(Copy, Clone, PartialEq)]
-#[cfg_attr(feature = "python-bindings", pyclass)]
-pub enum ApsMode {
-    /// Averaged over all data provided. New averages can be created by calling
-    /// `AvPowerSpectra::reset()`
-    AllAveraging {},
-    /// In this mode, the `AvPowerSpectra` works a bit like a sound level meter,
-    /// where new data is weighted with old data, and old data exponentially
-    /// backs off. This mode only makes sense when `tau >> nfft/fs`
-    ExponentialWeighting {
-        /// Time weighting constant, follows convention of Sound Level Meters.
-        /// Means the data is approximately low-pass filtered with a cut-off
-        /// frequency f_c of s/tau ≅ 1 → f_c = (2 * pi * tau)^-1.
-        tau: Flt,
-    },
-    /// Spectrogram mode. Only returns the latest estimate(s).
-    Spectrogram {},
-}
-impl Default for ApsMode {
-    fn default() -> Self {
-        ApsMode::AllAveraging {}
-    }
-}
-#[cfg(feature = "python-bindings")]
-#[cfg_attr(feature = "python-bindings", pymethods)]
-impl ApsMode {
-    #[inline]
-    fn __eq__(&self, other: &Self) -> bool {
-        self == other
-    }
-}
-
 /// Averaged power spectra computing engine
 /// Used to compute power spectra estimations on
 /// long datasets, where nfft << length of data. This way, the variance of a
 /// single periodogram is suppressed with increasing number of averages.
 ///
 /// For more information, see the book on numerical recipes.
-///
+
+#[cfg(feature = "python-bindings")]
+#[cfg_attr(feature = "python-bindings", pyclass)]
+#[derive(Debug)]
+
 pub struct AvPowerSpectra {
     // Power spectra estimator for single block
     ps: PowerSpectra,
 
+    // Settings for computing power spectra, see [ApsSettings]
     settings: ApsSettings,
 
     // The number of samples to keep in the time buffer when overlapping time
@@ -454,6 +281,28 @@ impl AvPowerSpectra {
     }
 }
 
+#[cfg(feature = "python-bindings")]
+#[cfg_attr(feature = "python-bindings", pymethods)]
+impl AvPowerSpectra {
+    #[new]
+    fn new_py(s: ApsSettings) -> AvPowerSpectra {
+        AvPowerSpectra::new(s)
+    }
+
+    #[pyo3(name = "compute")]
+    fn compute_py<'py>(
+        &mut self,
+        py: Python<'py>,
+        dat: PyArrayLike2<Flt>,
+    ) -> Bound<'py, PyArray3<Cflt>> {
+        let dat = dat.as_array();
+        if let Some(res) = self.compute_last(dat) {
+            let res = res.clone();
+            return res.to_pyarray_bound(py);
+        }
+        panic!("No data!");
+    }
+}
 #[cfg(test)]
 mod test {
     use approx::assert_abs_diff_eq;
diff --git a/src/ps/apsmode.rs b/src/ps/apsmode.rs
new file mode 100644
index 0000000..508a5e6
--- /dev/null
+++ b/src/ps/apsmode.rs
@@ -0,0 +1,35 @@
+use crate::config::*;
+/// The 'mode' used in computing averaged power spectra. When providing data in
+/// blocks to the [AvPowerSpectra] the resulting 'current estimate' responds
+/// differently, depending on the model.
+#[derive(Copy, Clone, PartialEq, Debug)]
+#[cfg_attr(feature = "python-bindings", pyclass)]
+pub enum ApsMode {
+    /// Averaged over all data provided. New averages can be created by calling
+    /// `AvPowerSpectra::reset()`
+    AllAveraging {},
+    /// In this mode, the `AvPowerSpectra` works a bit like a sound level meter,
+    /// where new data is weighted with old data, and old data exponentially
+    /// backs off. This mode only makes sense when `tau >> nfft/fs`
+    ExponentialWeighting {
+        /// Time weighting constant, follows convention of Sound Level Meters.
+        /// Means the data is approximately low-pass filtered with a cut-off
+        /// frequency f_c of s/tau ≅ 1 → f_c = (2 * pi * tau)^-1.
+        tau: Flt,
+    },
+    /// Spectrogram mode. Only returns the latest estimate(s).
+    Spectrogram {},
+}
+impl Default for ApsMode {
+    fn default() -> Self {
+        ApsMode::AllAveraging {}
+    }
+}
+#[cfg(feature = "python-bindings")]
+#[cfg_attr(feature = "python-bindings", pymethods)]
+impl ApsMode {
+    #[inline]
+    fn __eq__(&self, other: &Self) -> bool {
+        self == other
+    }
+}
diff --git a/src/ps/apssettings.rs b/src/ps/apssettings.rs
new file mode 100644
index 0000000..2045dcc
--- /dev/null
+++ b/src/ps/apssettings.rs
@@ -0,0 +1,168 @@
+use super::*;
+use crate::config::*;
+use anyhow::{bail, Error, Result};
+use derive_builder::Builder;
+
+/// All settings used for computing averaged power spectra using Welch' method.
+#[derive(Builder, Clone, Debug)]
+#[cfg_attr(feature = "python-bindings", pyclass)]
+#[builder(build_fn(validate = "Self::validate", error = "Error"))]
+pub struct ApsSettings {
+    /// Mode of computation, see [ApsMode].
+    #[builder(default)]
+    pub mode: ApsMode,
+    /// Overlap in time segments. See [Overlap].
+    #[builder(default)]
+    pub overlap: Overlap,
+    /// Window applied to time segments. See [WindowType].
+    #[builder(default)]
+    pub windowType: WindowType,
+    /// Kind of freqency weighting. Defaults to Z
+    #[builder(default)]
+    pub freqWeightingType: FreqWeighting,
+    /// FFT Length
+    pub nfft: usize,
+    /// Sampling frequency
+    pub fs: Flt,
+}
+
+impl ApsSettingsBuilder {
+    fn validate(&self) -> Result<()> {
+        if !self.fs.is_some() {
+            bail!("Sampling frequency not given");
+        }
+        let fs = self.fs.unwrap();
+
+        if !fs.is_normal() {
+            bail!("Sampling frequency not a normal number")
+        }
+        if fs <= 0.0 {
+            bail!("Invalid sampling frequency given as parameter");
+        }
+
+        if self.nfft.is_none() {
+            bail!("nfft not specified")
+        };
+        let nfft = self.nfft.unwrap();
+        if nfft % 2 != 0 {
+            bail!("NFFT should be even")
+        }
+        if nfft == 0 {
+            bail!("Invalid NFFT, should be > 0.")
+        }
+        // Perform some checks on ApsMode
+        if let Some(ApsMode::ExponentialWeighting { tau }) = self.mode {
+            if tau <= 0.0 {
+                bail!("Invalid time weighting constant [s]. Should be > 0 if given.");
+            }
+        }
+
+        Ok(())
+    }
+}
+
+#[cfg(feature = "python-bindings")]
+#[cfg_attr(feature = "python-bindings", pymethods)]
+impl ApsSettings {
+    #[new]
+    fn new(
+        mode: ApsMode,
+        overlap: Overlap,
+        windowType: WindowType,
+        freqWeightingType: FreqWeighting,
+        nfft: usize,
+        fs: Flt,
+    ) -> ApsSettings {
+        ApsSettings {
+            mode,
+            overlap,
+            windowType,
+            freqWeightingType,
+            nfft,
+            fs,
+        }
+    }
+}
+
+impl ApsSettings {
+    /// Returns the amount of samples to keep in overlapping blocks of power
+    /// spectra.
+    pub fn get_overlap_keep(&self) -> usize {
+        self.validate_get_overlap_keep().unwrap()
+    }
+
+    /// Returns the amount of samples to `keep` in the time buffer when
+    /// overlapping time segments using [TimeBuffer].
+    fn validate_get_overlap_keep(&self) -> Result<usize> {
+        let nfft = self.nfft;
+        let overlap_keep = match self.overlap {
+            Overlap::Number { N } if N >= nfft => {
+                bail!("Invalid overlap number of samples. Should be < nfft, which is {nfft}.")
+            }
+            // Keep 1 sample, if overlap is 1 sample etc.
+            Overlap::Number { N } if N < nfft => N,
+
+            // If overlap percentage is >= 100, or < 0.0 its an error
+            Overlap::Percentage { pct } if !(0.0..100.).contains(&pct) => {
+                bail!("Invalid overlap percentage. Should be >= 0. And < 100.")
+            }
+            // If overlap percentage is 0, this gives
+            Overlap::Percentage { pct } => ((pct * nfft as Flt) / 100.) as usize,
+            Overlap::NoOverlap {} => 0,
+            _ => unreachable!(),
+        };
+        if overlap_keep >= nfft {
+            bail!("Computed overlap results in invalid number of overlap samples. Please make sure the FFT length is large enough, when high overlap percentages are required.");
+        }
+        Ok(overlap_keep)
+    }
+
+    /// Return a reasonable acoustic default with a frequency resolution around
+    /// ~ 10 Hz, where nfft is still an integer power of 2.
+    ///
+    /// # Errors
+    ///
+    /// If `fs` is something odd, i.e. < 1 kHz, or higher than 1 MHz.
+    ///
+    pub fn reasonableAcousticDefault(fs: Flt, mode: ApsMode) -> Result<ApsSettings> {
+        if fs < 1e3 || fs > 1e6 {
+            bail!("Sampling frequency for reasonable acoustic data is >= 1 kHz and <= 1 MHz.");
+        }
+        let fs_div_10_rounded = (fs / 10.) as u32;
+
+        // 2^30 is about 1 million. We search for a two-power of an nfft that is
+        // the closest to fs/10. The frequency resolution is about fs/nfft.
+        let nfft = (0..30).map(|i| 2u32.pow(i) - fs_div_10_rounded).fold(
+            // Start wth a value that is always too large
+            fs as u32 * 10,
+            |cur, new| cur.min(new),
+        ) as usize;
+
+        Ok(ApsSettings {
+            mode,
+            fs,
+            nfft,
+            windowType: WindowType::default(),
+            overlap: Overlap::default(),
+            freqWeightingType: FreqWeighting::default(),
+        })
+    }
+
+    /// Return sampling frequency
+    pub fn fs(&self) -> Flt {
+        self.fs
+    }
+
+    /// Return Nyquist frequency
+    pub fn fnyq(&self) -> Flt {
+        self.fs / 2.
+    }
+
+    /// Returns a single-sided frequency array corresponding to points in Power
+    /// spectra computation.
+    pub fn getFreq(&self) -> Array1<Flt> {
+        let df = self.fs / self.nfft as Flt;
+        let K = self.nfft / 2 + 1;
+        Array1::linspace(0., (K - 1) as Flt * df, K)
+    }
+}
diff --git a/src/ps/fft.rs b/src/ps/fft.rs
index a7b708b..4433095 100644
--- a/src/ps/fft.rs
+++ b/src/ps/fft.rs
@@ -14,6 +14,11 @@ pub struct FFT {
     // nfft stored as float, this is how it is required most often
     nfftF: Flt,
 }
+impl Debug for FFT {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.debug_struct("Forward FFT engine, lenfth: {self.nfftF}").finish()
+    }
+}
 
 impl FFT {
     /// Create new FFT from given nfft
diff --git a/src/ps/freqweighting.rs b/src/ps/freqweighting.rs
index 869ecea..1d32cc0 100644
--- a/src/ps/freqweighting.rs
+++ b/src/ps/freqweighting.rs
@@ -25,6 +25,9 @@ impl FreqWeighting {
     fn all() -> Vec<Self> {
         Self::iter().collect()
     }
+    fn __str__(&self) -> String {
+        format!("{self}-weighting")
+    }
     #[staticmethod]
     #[pyo3(name = "default")]
     fn default_py() -> Self {
diff --git a/src/ps/mod.rs b/src/ps/mod.rs
index 9fa52e2..3d05bff 100644
--- a/src/ps/mod.rs
+++ b/src/ps/mod.rs
@@ -10,12 +10,16 @@ mod ps;
 mod timebuffer;
 mod window;
 mod freqweighting;
+mod apssettings;
 mod overlap;
+mod apsmode;
 use crate::config::*;
 
 
 pub use freqweighting::FreqWeighting;
 pub use overlap::Overlap;
-pub use aps::{ApsSettings, ApsSettingsBuilder,ApsMode, AvPowerSpectra};
+pub use apssettings::{ApsSettings, ApsSettingsBuilder};
+pub use apsmode::ApsMode;
+pub use aps::AvPowerSpectra;
 pub use ps::{CrossPowerSpecra, PowerSpectra, CPSResult};
 pub use window::{Window, WindowType};
diff --git a/src/ps/ps.rs b/src/ps/ps.rs
index 959a18d..8ae54ab 100644
--- a/src/ps/ps.rs
+++ b/src/ps/ps.rs
@@ -85,6 +85,7 @@ impl CrossPowerSpecra for CPSResult {
 /// example the computations of spectrograms, or Welch' method of spectral
 /// estimation.
 ///
+#[derive(Debug)]
 pub struct PowerSpectra {
     /// Window used in estimator. The actual Window in here is normalized with
     /// the square root of the Window power. This safes one division when
diff --git a/src/ps/timebuffer.rs b/src/ps/timebuffer.rs
index ec1e971..6e620f2 100644
--- a/src/ps/timebuffer.rs
+++ b/src/ps/timebuffer.rs
@@ -8,7 +8,7 @@ use std::collections::VecDeque;
 /// TimeBuffer, storage to add blocks of data in a ring buffer, that can be
 /// extracted by blocks of other size. Also, we can keep samples in a buffer to
 /// create, for example, overlapping windows of time data.
-#[derive(Default)]
+#[derive(Default, Debug)]
 pub struct TimeBuffer {
     data: Vec<VecDeque<Flt>>,
 }
diff --git a/src/ps/window.rs b/src/ps/window.rs
index 3384c3e..c3994f3 100644
--- a/src/ps/window.rs
+++ b/src/ps/window.rs
@@ -108,7 +108,7 @@ impl WindowType {
 
 
 /// Window (taper) computed from specified window type.
-#[derive(Clone)]
+#[derive(Clone, Debug)]
 pub struct Window {
     /// The enum from which it is generated
     pub w: WindowType,