From 4808e37f2ff4facadf6630ad377d58b676c314b6 Mon Sep 17 00:00:00 2001
From: Leonardo Robol <leo@robol.it>
Date: Mon, 1 Mar 2010 18:17:17 +0100
Subject: [PATCH] Aggiunto qualche script per il filtering e per un primo
 esempio di decomposizione mediante filter bank ortogonale di un segnale
 (precisamente Sunny di George Benson).

---
 Filtering/Filtering.py | 39 +++++++++++++++++++++++
 Filtering/Splitting.py | 86 ++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 125 insertions(+)
 create mode 100644 Filtering/Filtering.py
 create mode 100644 Filtering/Splitting.py

diff --git a/Filtering/Filtering.py b/Filtering/Filtering.py
new file mode 100644
index 0000000..5f336f3
--- /dev/null
+++ b/Filtering/Filtering.py
@@ -0,0 +1,39 @@
+#
+# -*- coding: utf-8 -*-
+#
+
+from numpy import array, convolve, zeros
+
+class AbstractFilter():
+
+    def __init__(self):
+        pass
+
+    def Apply (self, samples):
+        """
+        Apply the filter to samples. This method has to be
+        overloaded from the specific filter class
+        """
+        pass
+               
+
+class FIR(AbstractFilter):
+
+    def __init__(self, coefficients):
+        
+        self.coefficients = coefficients
+    
+    def Apply(self, samples):
+        """Apply the FIR to samples"""
+        return convolve (self.coefficients,
+                         samples, 'same')
+
+def DownSample(samples):
+
+    return samples[::2]
+
+def UpSample(samples):
+    
+    s = zeros(2 * len(samples))
+    s[::2] = samples
+    return s
diff --git a/Filtering/Splitting.py b/Filtering/Splitting.py
new file mode 100644
index 0000000..0875844
--- /dev/null
+++ b/Filtering/Splitting.py
@@ -0,0 +1,86 @@
+#
+# -*- coding: utf-8 -*-
+# 
+
+import Filtering, pylab
+from numpy import array, sqrt, memmap
+
+class SplittingExample():
+
+    def __init__(self):
+
+        # Creo i filtri che userò in maniera arguta
+        # in seguito. Questo filtro è ortogonale, nel
+        # senso che posso ottenere la filter bank inversa
+        # invertendo semplicemente la matrice. 
+        low = 0.125 * array([1+sqrt(3), 3+sqrt(3), 3-sqrt(3), 1-sqrt(3)])
+        high = 0.125 * array([1-sqrt(3), -3+sqrt(3), 3+sqrt(3),-1-sqrt(3)])
+        
+        self.lowpass  = Filtering.FIR (low)
+        self.highpass = Filtering.FIR (high)
+
+        self.LoadSamples ()
+        self.Split ()
+
+        # self.Show ()
+    
+    def LoadSamples(self):
+        """
+        Load the samples from an audio file
+        """
+        filename = "sunny.pcm"
+        samples = memmap (filename,
+                          dtype="<h",
+                          mode="r")
+        self.samples = samples
+
+    def WriteSamples(self, samples, filename):
+        data = memmap (filename,
+                       dtype="<h",
+                       mode="w+",
+                       shape = len(samples))
+        data[:] = samples[:]
+        data.flush ()
+        
+    def Split(self, suffix = "d"):
+        """
+        Split the signal in two signals filtered with the
+        lowpass and the highpass filter"""
+
+        # Primo giro di filtri
+        self.low = self.lowpass.Apply (self.samples)
+        self.high = self.highpass.Apply (self.samples)
+
+        self.low_d = Filtering.DownSample(self.low)
+        self.high_d = Filtering.DownSample(self.high)
+        
+
+
+        self.WriteSamples(self.high_d, "high_" + suffix + ".pcm")
+
+        # Applichiamo ricorsivamente il filtro alle basse frequenze,
+        # per suddividerle ancora meglio fra details and average.
+        self.samples = self.low_d
+
+        # Per il momento usiamo suffix per controllare il punto della
+        # ricorsione a quale ci troviamo. 
+        if len(suffix) < 4:
+            self.Split(suffix + "d")
+        else:
+            # Memorizziamo i low ed usciamo
+            self.WriteSamples(self.low_d,  "low_" + suffix + ".pcm" )
+
+
+    def Show(self):
+        """
+        Shows the result of filtering
+        """
+        pylab.plot (self.low)
+        pylab.plot (self.high)
+
+        pylab.show ()
+        
+        
+if __name__ == "__main__":
+
+    SplittingExample ()
-- 
2.1.4