# little test program to check window normalization factor for DFT
""" bitpack_vector.py --
-----------------------------------------------------------------------
C 2009 Marina Bosi & Richard E. Goldberg -- All rights reserved
-----------------------------------------------------------------------
vectorized code for packing and unpacking bits into an array of bytes"""

import numpy as np

# create the spectrum (times 2/N) for a sine-windowed sinusoid of unit amplitude
N= 2048
sw = np.sin(np.pi*(np.arange(N)+0.5)/N)  # sine window
peakamp=1.  # set to 1
freq = int(N/4)  # so as to align w/ a frequency line
signal = peakamp * np.cos(2.0 * np.pi * freq * np.arange(N) / N)  # just a unit amplitude sinusoid
spectrum= np.fft.fft(signal * sw)  # computes the DFT of sine-windowed sinusoid
spectrum = spectrum *2/N  # this takes care of everything but the window effect

# Craig computes normalization based on the single peak line
peak = max(np.abs(spectrum))

# We recommend summing intensity over peak when computing SPL
width = 1  # how many neighbors to include in sum over peak (I typically use 1 on either side of peak)
peakplus= np.abs(spectrum)[(freq-width):(freq+width+1)]  # peak plus those neighbors


# print results
print "Normalizaton Example: (2/N*DFT of sine-windowed unit-amplitude sinusoid)\n"
print "\nSpectral Peak: ",peak, " (Approx. 2/pi = ",2./np.pi," as noted by Craig)"  # the peak
print "\nSpectral Peak + ",width," neighbor(s) on each side: \n", peakplus  # the extended peak

print "\n\nSum over peak: ", peak**2, " (Approx. (2/pi)^2 = ",(2./np.pi)**2, " as noted by Craig)" # what Craig is summing over
print "\nSum over Broadened peak: ", sum(peakplus**2), " (Approx. <w^2> = 1/2 as in book p. 229 for sine window)"