@@ -1,8 +1,8 @@
 import itertools
 import matplotlib.pyplot as plt
 import numpy as np
-from scipy.io import wavfile
 from scipy.signal import detrend, lfilter, bilinear, spectrogram, filtfilt, resample, fftconvolve
+import soundfile as sf
 import acoustics
 
 from acoustics.standards.iso_tr_25417_2007 import REFERENCE_PRESSURE
@@ -1019,35 +1019,32 @@ def normalize(self, gap=6.0, inplace=False):
         else:
             return self / factor[..., None]
 
-    def to_wav(self, filename, depth=16):
+    def to_wav(self, filename, depth=16, format="WAV"):
         """Save signal as WAV file.
 
         :param filename: Name of file to save to.
         :param depth: If given, convert to integer with specified depth. Else, try to store using the original data type.
+        :param format: It can be either WAV or FLAC
 
-        By default, this function saves a normalized 16-bit version of the signal with at least 6 dB range till clipping occurs.
+        Note that the floating point samples are not yet supported.
 
         """
-        data = self
-        dtype = data.dtype if not depth else 'int' + str(depth)
-        if depth:
-            data = (data * 2**(depth - 1) - 1).astype(dtype)
-        wavfile.write(filename, int(self.fs), data.T)
-        #wavfile.write(filename, int(self.fs), self._data/np.abs(self._data).max() *  0.5)
-        #wavfile.write(filename, int(self.fs), np.int16(self._data/(np.abs(self._data).max()) * 32767) )
+        subtype = "PCM_{}".format(depth)
+        sf.write(filename, self.T, int(self.fs), format=format, subtype=subtype)
 
     @classmethod
-    def from_wav(cls, filename):
+    def from_wav(cls, filename, normalize=True):
         """
         Create an instance of `Signal` from a WAV file.
 
         :param filename: Filename
+        :param normalize: Optionally normalize the input samples
 
         """
-        fs, data = wavfile.read(filename)
-        data = data.astype(np.float32, copy=False).T
-        data /= np.max(np.abs(data))
-        return cls(data, fs=fs)
+        data, fs = sf.read(filename)
+        if normalize:
+            data /= np.max(np.abs(data))
+        return cls(data.T, fs=fs)
 
 
 _PLOTTING_PARAMS = {
 
@@ -23,6 +23,7 @@
 , bootstrapped-pip
 , stdenv
 , python
+, pysoundfile
 , development ? false
 }:
 
@@ -52,7 +53,7 @@ in buildPythonPackage rec {
 
   checkInputs = [ pytest glibcLocales ];
   nativeBuildInputs = lib.optionals development [ sphinx pylint yapf ];
-  propagatedBuildInputs = [ cytoolz numpy scipy matplotlib pandas six tabulate ];
+  propagatedBuildInputs = [ cytoolz numpy scipy matplotlib pandas six tabulate pysoundfile ];
 
   meta = {
     description = "Acoustics module for Python";
 
@@ -4,3 +4,4 @@ matplotlib
 six
 pandas
 tabulate
+pysoundfile
@@ -24,6 +24,7 @@ install_requires =
     six >= 1.4.1
     pandas >= 0.15
     tabulate
+    pysoundfile
 tests_require =
     pytest
 
 
@@ -20,11 +20,14 @@ class test_wav():
     fs = 10025
     samples = int(fs * duration)
     channels = 3
+    values = np.random.randn(channels, samples)
 
-    signal = Signal(np.random.randn(channels, samples), fs)
+    signal = Signal(values, fs)
+    signal.normalize(inplace=True)
 
     with tempfile.TemporaryFile() as file:
         signal.to_wav(file)
+        file.seek(0)
         signal = Signal.from_wav(file)
         assert signal.samples == samples
         assert signal.fs == fs