#!/usr/bin/env python # # Programmer: Craig Stuart Sapp # Creation Date: Mon Feb 1 21:31:24 PST 2010 # Last Modified: Wed Feb 17 21:25:28 PST 2010 # Filename: craigwindow.py # Syntax: Python 2.6 # """ Signal windowing functions for Music 422. Implementations of the Sine Window, the Hanning Window for both even and odd window lengths, as well as the Kaiser-Bessel Derived Window. Also contains a test function called ColaMirrorTest() which can be used to check if a particular window will generate a constant-amplitude reconstruction at 50% overlap between adjacent windows. Compatible importing method: from craigwindow import ApplySineWindow as SineWindow from craigwindow import ApplyHanningWindow as HanningWindow from craigwindow import ApplyKBDWindow as KBDWindow Although it would be more efficient to not regenerate the window values for each application of the window: from craigwindow import SineWindow from craigwindow import HanningWindow from craigwindow import KBDWindow And then store the precalculated window in an array to apply to signals with a multiply, for example: import random signal = [random.random() for i in xrange(1024)] window = SineWindow(1024) windowedsignal = signal * window """ import numpy as np def ApplySineWindow(signal): return signal * SineWindow(len(signal)) def ApplyHanningWindow(signal): return signal * HanningWindow(len(signal)) def ApplyKBDWindow(signal, alpha): return signal * KBDWindow(len(signal), alpha) ########################################################################### ############################## ## ## SineWindow -- page 107, equation 2 (Even-length sine window definition). ## Also the odd-length window definition is implemented by this function ## suitable for 50% overlap-add (derived from continuous sine window ## definition in equation 1 on page 107). ## def SineWindow(length): """Sine Window -- Generates a sine window which, when 50% overlapped, generates a constant overlap-add (COLA) amplitude when the window is squared (applied twice to a signal). Shift of 1/2 sample for even lengths because window is sampled from the middle of the window, not from left end of window. Similarly, scaling by N-1 in the odd case such that the continuous window definition is sampled symmetrically around the middle of the window. Argument: length: Length of the output window array. """ N = int(length) assert N > 1 if (N%2 == 0): halfwindow = np.sin((np.arange(N/2)+0.5)*np.pi/N) return np.concatenate([halfwindow, np.flipud(halfwindow)]) else: halfwindow = np.sin(np.arange(N/2+1)*np.pi/(N-1)) # don't repeat the middle sample of an odd window: return np.concatenate([halfwindow, np.flipud(halfwindow[0:-1])]) ############################## ## ## HanningWindow -- page 107, last equation. Even-length window definition. ## Odd-length version also implemented in this function. ## def HanningWindow(length, power=2): """Hanning Window -- Raised Cosine, Cosine**2 or Sine**2 window. Arguments: length: Length of the output window array which must be positive. power: power=2 for Hanning window, power=1 for Sine Window. """ return SineWindow(length)**power ############################## ## ## KBDWindowKVH -- from Katarina Van Heusen's solution, implemented ## using numpy.kaiser function, and then deriving KBD. ## def KBDWindowKVH(length, beta): """Kaiser Bessel Derived Window Arguments: length: Length of output window array. Must be positive & even. beta: Kaiser window parameter alpha times pi. """ N = int(length) assert N%2 == 0 assert N > 0 M = N/2 w = np.kaiser(M+1, beta) wMsum = 1.0 / np.sqrt(w[0:M+1].sum()) wnsuml = np.sqrt(np.cumsum(w[0:M])) wl = wnsuml * wMsum return np.concatenate([wl, np.flipud(wl)]) def KBDWindow(length, alpha): return KBDWindowKVH(length, np.pi * alpha) ############################## ## ## KBDWindowCSS -- Craig's solution, implemented from bessel I_0 function. ## Based on C program implementation: ## https://ccrma.stanford.edu/courses/422/projects/kbd/kbdwindow.cpp ## def KBDWindowCSS(length, alpha): """Kaiser Bessel Derived Window Arguments: length: Length of output window array. Must be positive & even. alpha: Kaiser window parameter. """ N = int(length) assert N % 2 == 0 assert N > 0 halfN = N/2 beta = np.pi * alpha bessels = np.i0(beta * np.sqrt(1.0-(4.0*np.arange(halfN)/N-1.0)**2.0)) halfwin = np.cumsum(bessels) normalize = halfwin[-1] + np.i0(beta*np.sqrt(1.0-(4.0*halfN/N-1.0)**2.0)) halfwin = np.sqrt(halfwin/normalize) return np.concatenate([halfwin, np.flipud(halfwin)]) ############################## ## ## KBDWindowCAR -- from Colin Raffel's solution: Implements Kaiser window, then ## Applys KBD modification. ## def kaiserWindow(length, alpha): """Kaiser-Bessel Window Arguments: length: Length of output window array. Must be positive & even. alpha: Kaiser window parameter. """ N = int(length) n = np.arange(N) w = (n - (N-1)/2.0)/((N-1)/2.0) w = w**2 w = np.sqrt(1-w) w = np.pi * alpha * w w = np.i0(w) w = w/np.i0(np.pi*alpha) return w def KBDWindowCAR(length, alpha): """Kaiser-Bessel Derived Window Arguments: length: Length of output window array. Must be positive & even. alpha: Kaiser window parameter. """ N = int(length) assert N%2 == 0 assert N > 0 wKai = kaiserWindow(N/2+1, alpha) wHalf = np.sqrt(np.cumsum(wKai[0:-1])/np.sum(wKai)) return np.concatenate([wHalf, np.flipud(wHalf)]) ############################## ## ## ColaMirrorTest -- ## def ColaMirrorTest(window, power = 2): """ColaMirrorTest -- Test for Constant OverLap-Add requirement in TDAC process which is that the square of the window values plus the mirror image of each half window squared sums to one (or at least a constant). The central value of an odd-length window will be checked agains both ends of the window (resulting in an additional middle value in the output array). Check the max an min values of the output array to verify that they are equal (or extremely close). Arguments: window: array containing a window. power: power value to apply to window values. (1 for hanning, 2 for sine or kbd) """ if (type(window) == list): window = np.array(window) N = window.size power = float(power) half = int(N/2.0) output = np.zeros(N, dtype=float) powed = window**power if (N % 2 == 0): firsthalf = powed[0:N/2] secondhalf = powed[N/2:N] return np.concatenate([firsthalf + np.flipud(firsthalf),\ secondhalf + np.flipud(secondhalf)]) else: # returns length N+1, since checking middle sample twice firsthalf = powed[0:N/2+1] secondhalf = powed[N/2:N] return np.concatenate([firsthalf + np.flipud(firsthalf),\ secondhalf + np.flipud(secondhalf)]) ########################################################################### # # Testing code: # if __name__ == "__main__": """ # Plot various windows with 1024 length ############################## #"" import pylab as pl N = 1024 n = np.arange(N)/(N-1.0) # normalized range lines = pl.plot(n, SineWindow(N), \ n, KBDWindow(N, 1), \ n, HanningWindow(N), \ n, KBDWindow(N, 4)) pl.legend(lines, ("Sine Window", \ "KBD Win alpha=1", \ "Hanning Window", \ "KBD Win alpha=4"), \ "lower center") pl.axis('tight') pl.xlabel('normalized windowlength') pl.ylabel('amplitude') pl.suptitle('Four examples of windows') pl.show() #"" # Test windows to make sure they are COLA: ############################ #"" sinewin = SineWindow(1024) sinecola = ColaMirrorTest(sinewin) print "Min value:", np.min(sinecola) print "Max value:", np.max(sinecola) ### Min value: 1.0 ### Max value: 1.0 hanwin = HanningWindow(1024) hancola = ColaMirrorTest(hanwin, 1.0) print "Min value:", np.min(hancola) print "Max value:", np.max(hancola) ### Min value: 1.0 ### Max value: 1.0 kai100win = KBDWindowCSS(1024, 1.00) kai100cola = ColaMirrorTest(kai100win) print "Min value:", np.min(kai100cola) print "Max value:", np.max(kai100cola) ### Min value: 1.0 ### Max value: 1.0 kai200win = KBDWindowKVH(1024, 2.00) kai200cola = ColaMirrorTest(kai200win) print "Min value:", np.min(kai200cola) print "Max value:", np.max(kai200cola) ### Min value: 1.0 ### Max value: 1.0 kai399win = KBDWindowCAR(1024, 3.99) kai399cola = ColaMirrorTest(kai399win) print "Min value:", np.min(kai399cola) print "Max value:", np.max(kai399cola) ### Min value: 1.0 ### Max value: 1.0 # Make a transitional window from Sine into KBD: N = 1024 hywin = np.concatenate([SineWindow(N)[0:N/2], KBDWindow(N,4)[N/2:N]]) hywincola = ColaMirrorTest(hywin) print "Min value:", np.min(hywincola) print "Max value:", np.max(hywincola) ### Min value: 1.0 ### Max value: 1.0 #"" # Comparison test between KBDWindowKVH and KBDWindowCSS and KBDWindowCAR: #"" KBDWindowCSS(2048, 1.0)[0:10] ### array([0.01629672, 0.02310258, 0.02836276, 0.03282897, 0.0367916, ### 0.04039919, 0.04373976, 0.04687052, 0.04983111, 0.05265029]) max(np.abs(KBDWindowKVH(2048, np.pi) - KBDWindowCSS(2048, 1.0))) ### 6.6613381477509392e-16 max(np.abs(KBDWindowCAR(2048, np.pi) - KBDWindowCSS(2048, 1.0))) ### 6.6613381477509392e-16 max(np.abs(KBDWindowCAR(2048, 1.0) - KBDWindowKVH(2048, np.pi))) ### 2.2204460492503131e-16 import time t1 = time.clock() csswin = KBDWindowCSS(1000000, 1.0) t2 = time.clock() print "CSS TIME:", t2 - t1, "seconds" ### CSS TIME: 0.90310835761 seconds t1 = time.clock() kvhwin = KBDWindowKVH(1000000, 1.0) t2 = time.clock() print "KVH TIME:", t2 - t1, "seconds" ### KVH TIME: 0.900658371609 seconds t1 = time.clock() carwin = KBDWindowCAR(1000000, 1.0) t2 = time.clock() print "CAR TIME:", t2 - t1, "seconds" ### CAR TIME: 0.905957039338 seconds #"" # Test window-application functions which internally generate a window. #"" import pylab as pl N = 1024 signal = np.sin(2.0 * np.pi * 3000.0 * np.arange(N) / 44100.0) n = np.arange(N)/(N-1.0) # normalized range pl.clf() lines = pl.plot(n, ApplySineWindow(signal), \ n+1, ApplyHanningWindow(signal), \ n+2, ApplyKBDWindow(signal, 1), \ n+3, ApplyKBDWindow(signal, 4)) pl.axis('tight') pl.xlabel('normalized windowlength') pl.ylabel('amplitude') pl.suptitle('Various windows applied to a sinewave') pl.text(0.5, 0.0, 'sine\nwindowed', horizontalalignment='center', verticalalignment='center') pl.text(1.5, 0.0, 'hanning\nwindowed', horizontalalignment='center', verticalalignment='center') pl.text(2.5, 0.0, 'KBD windowed\nalpha=1', horizontalalignment='center', verticalalignment='center') pl.text(3.5, 0.0, 'KBD windowed\nalpha=4', horizontalalignment='center', verticalalignment='center') pl.show() #"" # Visualize constant overlapp-add of windows #"" import pylab as pl pl.clf() N = 1024 window = KBDWindow(1024, 4.0) winl = window[0:N/2] winr = window[N/2:N] overlapsum = winl + winr win2l = winl**2 win2r = winr**2 overlap2sum = win2l + win2r n = np.arange(N)/(N-1.0) # normalized range nl = n[0:N/2] nr = n[N/2:N] pl.plot(nl, winl, 'b-') pl.plot(nl, winr, 'b-') pl.plot(nl, overlapsum, 'r--') pl.plot(nl+1, win2l, 'b-') pl.plot(nl+1, win2r, 'b-') pl.plot(nl+1, overlap2sum, 'r--') pl.text(0.25, 1.1, 'not COLA$^1$', \ horizontalalignment='center', verticalalignment='center') pl.text(1.25, 1.1, 'COLA$^2$', \ horizontalalignment='center', verticalalignment='center') pl.text(0.5, 0.5, 'KBD Window, $\\alpha=4$\n50% overlap', \ horizontalalignment='center', verticalalignment='center') pl.text(1.0, 0.5, 'KBD Window$^2$, $\\alpha=4$\n50% overlap', \ horizontalalignment='center', verticalalignment='center') pl.show() """