numpy binary file output

README.md

@@ -31,27 +31,24 @@ Mixed-Phase Compensation references:
 Required Python dependencies:
 
     1) Python 2.7
-    2) Scientific Python: SciPy (v0.11+), Numpy, & Matplotlib
-    3) Scikits.audiolab, setuptools, libsndfile
+    2) Scientific Python: SciPy, Numpy, & Matplotlib
 
-This is a command line tool. I will be working to minimuze the number of dependencies in 
-future releases. One may successfully install PORC on Windows, although some extra legwork is 
-required to get all the dependencies and environmental variables set (email me for help).
+The easiest install method on Windows is simply to install the continuum.io Anaconda package.
 
 Measurement
 ===========
 
 One needs to measure the log-frequency impulse response of your speakers with a 
 calibrated Electret Measurement Microphone, e.g. Dayton Audio EMM-6. Software 
-such as Room EQ Wizard (REQ) may be used for this purpose:
+such as Room EQ Wizard (REQ), Holm Impulse, or Arta may be used for this purpose:
 http://www.hometheatershack.com/roomeq/
 
 Usage
 =====
 
-porc.py [-h] [--mixed] [-t FILE] [-n NTAPS] input_file output_file
+porc.py [-h] [--mixed] [-t FILE] [-n NTAPS] [-o OPFORMAT] input_file output_file
 
-    python porc.py -t tact30f.txt -n 6144 l48.wav leq48.wav
+    python porc.py -t tact30f.txt -n -o bin 6144 l48.wav leq48.bin
 
 Use the -h flag for help!
 
@@ -81,7 +78,7 @@ easily do this with Audacity or REQ.
 
 Example:
 
-	python porc.py --mixed -t tact30f.txt -n 6144 l48.wav leq48.wav
+	python porc.py --mixed -t tact30f.txt -n -o bin 6144 l48.wav leq48.wav
 
 Have some patience with this method. The convolution takes a few CPU cycles.
 
@@ -101,17 +98,15 @@ You may need to merge left and right channels into a single stereo .wav
 OpenDRC Convolution
 ===================
 
-Use sox to convert output .wav to raw 32 bit IEEE floating point mono for the left & right channels
+Use -o bin flag to set binary 32bit IEEE floating point mono file format output for OpenDRC.
 
-    sox leq48.wav -t f32 leq48.bin
-    sox req48.wav -t f32 req48.bin
 
 TODO
 ====
 
-	Implement algo to automatically remove leading silence from RIR.
-	Port this code to C 
-    Update this page with better documentation!
+	Implement algo to automatically remove leading silence (zeros) from RIR.
+	Add a GUI Frontend (pretty interactive graphs, drawing target curve, etc...)
+	Update this page with better documentation!
 
 Contact
 =======

porc.py

@@ -42,6 +42,9 @@
 # Python libs
 import sys
 import textwrap
+import wave
+from contextlib import closing
+import struct
 
 # Scipy, Numpy, and matplotlibs
 import numpy as np
@@ -55,9 +58,6 @@
 from scipy.stats import kurtosis, nanstd
 from scipy.stats import norm as Gaussian
 import matplotlib.pyplot as plt
-from scikits.audiolab import Format, Sndfile
-
-import sys
 
 # PORC source files
 from parfiltid import parfiltid
@@ -98,7 +98,7 @@ def mad(a, c=Gaussian.ppf(3/4.), axis=0):  # c \approx .6745
     a = np.asarray(a)
     return np.median((np.fabs(a))/c, axis=axis)
 
-def roomcomp(impresp, filter, target, ntaps, mixed_phase):
+def roomcomp(impresp, filter, target, ntaps, mixed_phase, opformat):
 
   # Read impulse response
   Fs, data = wavfile.read(impresp)
@@ -218,25 +218,23 @@ def roomcomp(impresp, filter, target, ntaps, mixed_phase):
       #eqresp = np.real(conv(equalizer, data))
     else:
       print "zero taps; skipping mixed-phase computation"
-
-  # TODO: Fix the scipi.io wavfile.write method?
-  # For whatver reason, I can't get Scipy's wav io to work with
-  # sox (e.g. sox filter.wav -t f32 filter.bin) and OpenDRC's file import.
-  # Audiolab works well below, but it's an extra set of dependencies
-  #wavfile.write(filter, Fs, equalizer.astype(np.float16))
-
-  # Write data, convert to 16 bits
-  format = Format('wav')
-  f = Sndfile(filter, 'w', format, 1, Fs)
-  f.write_frames(norm(np.real(equalizer)))
-  f.close
-
-  #f = Sndfile('eqresp.wav', 'w', format, 1, Fs)
-  #f.write_frames(norm(eqresp))
-  #f.close
-
-  print '\nOutput filter length =', len(equalizer), 'taps'
-  print 'Output filter written to ' + filter
+  if opformat == 'wav':      
+  # Write data
+    wavwrite_24(filter, Fs, equalizer)
+
+    print '\nOutput filter length =', len(equalizer), 'taps'
+    print 'Output filter written to ' + filter
+
+  elif opformat == 'bin':
+    # direct output to bin avoids float64->pcm16->float32 conversion by going direct 
+    #float64->float32
+    f = open(filter, 'w')
+    norm(np.real(equalizer)).astype('float32').tofile(f)
+    f.close()
+    print '\nOutput filter length =', len(equalizer), 'taps'
+    print 'Output filter written to ' + filter
+  else:
+    print 'Output format not recognized, no file generated.'
 
   print "\nUse sox to convert output .wav to raw 32 bit IEEE floating point if necessary,"
   print "or to merge left and right channels into a stereo .wav"
@@ -281,6 +279,15 @@ def roomcomp(impresp, filter, target, ntaps, mixed_phase):
   plt.legend()
   plt.show()
 
+def wavwrite_24(fname, fs, data):
+    data_as_bytes = (struct.pack('<i', int(samp*(2**23-1))) for samp in data)
+    with closing(wave.open(fname, 'wb')) as wavwriter:
+        wavwriter.setnchannels(1)
+        wavwriter.setsampwidth(3)
+        wavwriter.setframerate(fs)
+        for data_bytes in data_as_bytes:
+            wavwriter.writeframes(data_bytes[0:3])
+
 def main():
 
 	print
@@ -313,13 +320,13 @@ def main():
 	parser.add_argument("-n", dest="ntaps", default = 6144,
 					  help="filter length, in taps. Default = len(input)", type=int)
 	parser.add_argument('--mixed', action='store_true', default = False,
-					  help="implement mixed-phase compensation. see README for details")                 
+					  help="implement mixed-phase compensation. see README for details") 
+	parser.add_argument("-o", dest="opformat", default = 'bin',
+					  help="Output file type, default bin optional wav", type=str)                     
 
 	args = parser.parse_args()
 
-	roomcomp(args.impresp, args.filter, args.target, args.ntaps, args.mixed)
+	roomcomp(args.impresp, args.filter, args.target, args.ntaps, args.mixed, args.opformat)
 
 if __name__=="__main__":
     main()  
-
-