"""
spect.py - Play with getting room frequency response from pink noise excitation.
================================================================================
This reads a record of a microphone's output while listening to pink noise and works
out the frequency response of the room (sort of, with *many* caveats) by finding the
power spectrum of the recorded noise and integrating that over frequency bands which
are fractions of an octave.
Pink noise has equal energy per unit time (power) per octave, so integrating that in
bands that are fixed fractions of an octave should give a constant value for each band
if the room response was flat. Deviations from flat are due to the room.
While this is fine as far as it goes, there are many additional things to take in to
account.
- The room response is really due to reflections of the sound from the walls (etc.)
and the consequent comb filtering as the reflections interact. For a fixed
listening position, placing absorbing material at the reflection points (determined
by "ray tracing" knowing the speaker locations) together with (large amounts of)
absorbing material in corners to soak up low frequency energy (which is not very
directional) is the proper way to "flatten" out the response. However - this correct
way of dealing with the problem is very often not possible in practice. Either the
layout of the room (locations of doors, windows, etc.) and hence far from ideal
speaker locations and listening positions, or the perfectly reasonable
objections of one's "significant other", or both, mean that playing around with the
frequency response by adjusting the signals that go in to the speakers is the only
way forward - regardless of its inadequacy. The room response measured with pink
noise is a reasonable way of determining how to adjust those signals (it is a starting
point, anyway).
- At higher frequencies, the spectrum will be exquisitely sensitive to the position at
which it is measured due to comb filtering. Movements of a few centimeters can convert
troughs to peaks and vice versa.
- In addition to the "constant in time" frequency response, the time it takes specific
frequencies to die away (reverberation) is a major influence on how much the frequencies
are emphasised (or otherwise) to a listener. This "time dependent" response can't be
measured by analysing the response to (constant amplitude) pink noise, of course.
This short bit of code is mainly for my own education. It does make you appreciate just
how complex this all really is. It also shows (again) how much can be done with a small
amount of Python.
Nick Glazzard, May 2014.
------------------------
"""
from numpy import *
from time import sleep
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
import wave
def power_spectrum( nsamples, data, deltat ):
"""Given nsamples of real voltage data spaced deltat seconds apart, find the spectrum of the
data (its frequency components). Square the amplitude for each frequency to convert to power
(intensity). First, window the data to make it pseudo-cyclic with a Hanning window.
Return the number of frequencies, frequency step, maximum frequency and the power spectrum."""
data_freq = fft.rfft(data * hanning( nsamples )) # Real to complex.
nfreqs = data_freq.size
data_freq = data_freq / nfreqs # Normalize.
spectrum = absolute(data_freq) ** 2 # Magnitude squared. (Complex to real).
freq_step = 1.0 / (deltat * 2 * nfreqs) # Frequency distance between spectrum samples.
max_freq = nfreqs * freq_step # Highest frequency for which spectrum is returned.
return( nfreqs, freq_step, max_freq, spectrum )
def logsteps(lo,hi,nperoctave):
"""Work out the log step size and number of steps to take to cover the frequency
range [lo:hi] with nperoctave steps per octave."""
r = hi - lo
step = log(2.0) / nperoctave
nsteps = int(nperoctave * (log(r)/log(2.0))) + 1
return step, nsteps
def logband(lo,step,istep):
"""Find the bounds of a frequency band given a minimum frequency, step size
and the step number for the band."""
l = exp(istep*step) - 1 + lo
h = exp((istep+1)*step) - 1 + lo
return l, h
def sumoctaves( nfreqs, freq_step, spectrum, lo, hi, nperoctave ):
"""Integrate the spectrum over frequency ranges that are fixed fractions of octaves.
Use frequency range lo:hi with nperoctave bands per octave."""
sums = []
mids = []
# Get a table of frequencies for each spectrum point.
freqs = arange( 0, nfreqs*freq_step, freq_step )
# Find the total number of steps and log step size to cover lo:hi with nperoctave
# steps in each octave.
step, nsteps = logsteps(lo,hi,nperoctave)
# Prepare to interpolate spectrum at arbitrary frequencies.
f = interp1d( freqs, spectrum )
# Step over the sub-octave bands
for i in range(0,nsteps):
# Get the frequency range of this band
bandlo, bandhi = logband(lo,step,i)
if bandlo < lo:
bandlo = lo
if bandhi > hi:
bandhi = hi
# Integrate the spectrum over the band.
freqfine = linspace(bandlo,bandhi) # Get 50 evenly spaced values in range bandlo:bandhi
spectrumfine = f(freqfine) # Interpolate the spectrum at those values
bandsum = trapz( spectrumfine, freqfine ) # Integrate the interpolated spectrum chunk
# Store the log of the integral and the mid band frequency.
sums.append(bandsum)
mids.append(0.5*(bandlo+bandhi))
return mids, sums
def readtextdata(n,c=1):
"""Read amplitude data from a text file, n, as written by Audacity with no headers.
Take the data from the left or right channel, according to c (0,1)."""
d = []
# Read the data file.
f = open(n)
for line in f:
w = line.split()
d.append(float(w[c]))
f.close()
return d
def getspectrum(d,s=44100,lo=20.0,hi=20000.0,nsuboct=3):
"""Analyze amplitude data, d.
The sample rate (samples per second) is s. Find the data's power spectrum,
then integrate this over the frequency range lo:hi in nsuboct bands per octave.
Return the mid frequncies and integrated power for each such band."""
# Find the power spectrum.
nsamples = len(d)
deltat = 1.0 / float(s)
nfreqs, freq_step, max_freq, spectrum = power_spectrum( nsamples, d, deltat )
print 'Spectrum: Samples:{0}, Sample width:{1:.3f}Hz, Max freq:{2}Hz.'.format(nfreqs,freq_step,max_freq)
# Integrate over sub-octave bands.
mids, sums = sumoctaves( nfreqs, freq_step, spectrum, lo, hi, nsuboct )
return asarray(mids), asarray(sums)
def readwave(n,c=1):
"""Read an entire WAV format file called n and return one of its channels c (0,1)."""
# Open the WAV file for read. Get some basic parameters and ensure sanity.
print 'Reading WAV file ', n
wr = wave.open(n,'rb')
nchannels = wr.getnchannels()
assert nchannels == 2
sampwidth = wr.getsampwidth()
assert sampwidth == 2
framerate = wr.getframerate()
frames = wr.getnframes()
print 'Channels:',nchannels
print 'Sample width:',sampwidth
print 'Frame rate:',framerate
print 'Frames:',frames
# Read all the data as bytes.
bd = wr.readframes(frames)
# Convert to signed 16 bit integers.
signal = frombuffer(bd,dtype=' n:
return None
else:
end = start + length
if end > n:
end = n
fsig = wave16[start:end].astype(float)
return fsig * (1.0/32768.0)
def getaveragedlogspectrum(n, chunksize, nsuboct, lo, hi):
"""Read the WAV file called n in chunksize blocks and get the integrated power in
nsuboct suboctave bands over the frequency range lo:hi. The chunksize determines
the internal frequency resolution at which the spectrum is computed."""
rate, nsamp, d = readwave(n)
start = 0
end = chunksize
navg = 0
# Iterate over all possible full chunks.
while end < nsamp:
navg += 1
print 'Chunk number:{0}, Sample range:{1} to {2}'.format(navg, start, end)
# Get a chunk of the amplitude data and find the spectral power in
# the sub-octave bands of that chunk.
chunk = getwavefloatchunk(d,start,chunksize)
mids, sums = getspectrum(chunk,rate,lo,hi,nsuboct)
# Sum the spectra for each chunk.
if start == 0:
avgspectrum = sums
else:
avgspectrum += sums
start = end
end += chunksize
# Convert spectra sum to average and convert to log power relative to average power
# of the bands.
if navg > 0:
avgspectrum /= float(navg)
avgval = sum(avgspectrum)/len(avgspectrum)
logavgspectrum = 10.0 * log10(avgspectrum/avgval)
return mids, logavgspectrum, (float(nsamp)/rate)
else:
return None, None, None
def expspectrum(logspect):
"""Convert dB log spectrum back to linear."""
alogspect = asarray(logspect)
return exp(0.23026*alogspect)
def logspectrum(linspect):
"""Convert linear to dB log spectrum."""
alinspect = asarray(linspect)
return 10.0 * log10(alinspect)
if True:
do_before_after = False # Uncorrected and corrected data OR corrected at different locations.
plot_all = False # Plot all curves for different locations OR average of them.
nsuboct = 3 #12 #24 #48 #12 #3 #48 #24 #3 #24
hi = 20000 #20000 #400 #20000.0 #500.0
if do_before_after:
mids1, logavgspectrum1, secs1 = getaveragedlogspectrum('longbypass.wav', 262144, nsuboct, 20.0, hi)
p1, = plt.plot( mids1, logavgspectrum1 )
mids2, logavgspectrum2, secs2 = getaveragedlogspectrum('longx.wav', 262144, nsuboct, 20.0, hi)
p2, = plt.plot( mids2, logavgspectrum2 )
secs = 0.5 * ( secs1 + secs2 )
else:
mids1, logavgspectrum1, secs1 = getaveragedlogspectrum('leftseat.wav', 262144, nsuboct, 20.0, hi)
mids2, logavgspectrum2, secs2 = getaveragedlogspectrum('leftmid.wav', 262144, nsuboct, 20.0, hi)
mids3, logavgspectrum3, secs3 = getaveragedlogspectrum('rightmidr.wav', 262144, nsuboct, 20.0, hi)
mids4, logavgspectrum4, secs4 = getaveragedlogspectrum('right.wav', 262144, nsuboct, 20.0, hi)
mids5, logavgspectrum5, secs5 = getaveragedlogspectrum('longx.wav', 262144, nsuboct, 20.0, hi)
secs = ( secs1 + secs2 + secs3 + secs4 + secs5 ) / 5.0
spaceavgspectrum = expspectrum(logavgspectrum1)
spaceavgspectrum += expspectrum(logavgspectrum2)
spaceavgspectrum += expspectrum(logavgspectrum3)
spaceavgspectrum += expspectrum(logavgspectrum4)
spaceavgspectrum += expspectrum(logavgspectrum5)
spaceavgspectrum = logspectrum(0.2*spaceavgspectrum)
if plot_all:
p1, = plt.plot( mids1, logavgspectrum1 )
p2, = plt.plot( mids2, logavgspectrum2 )
p3, = plt.plot( mids3, logavgspectrum3 )
p4, = plt.plot( mids4, logavgspectrum4 )
p5, = plt.plot( mids5, logavgspectrum5 )
else:
p6, = plt.plot( mids1, spaceavgspectrum )
secstr = ' Over {0:.1f} secs'.format(secs)
if nsuboct == 3:
teth = 'rd'
else:
teth = 'th'
title = 'Pink Noise Room Responses at 1/'+str(nsuboct)+teth+' Octave Bands.'
plt.title(title)
plt.xlabel( 'Freq (Hz)' )
plt.xlim(20.0, hi )
plt.ylim(-20.0,20.0)
plt.ylabel( 'Relative Response (dB)' )
plt.xscale('log')
plt.grid( b=True, which='major', linestyle='-' )
plt.grid( b=True, which='minor', linestyle='--' )
if do_before_after:
plt.legend([p1,p2], ['Before','After'])
else:
if plot_all:
plt.legend([p1,p2,p3,p4,p5], ['Left seat','Left middle','Right middle','Right','Center'])
else:
plt.legend([p6],['Position average'])
extradata = '20:'+str(int(hi))+'Hz.'+secstr
plt.text(25,17,extradata)
if do_before_after:
plt.savefig('before_after.png',bbox_inches='tight')
else:
plt.savefig('various_positions.png',bbox_inches='tight')
plt.show()
if False:
d = readtextdata('bypass2.txt')
mids, sums = getspectrum(d,nsuboct=3)
plt.plot( mids, sums )
d = readtextdata('corrected2.txt')
mids2, sums2 = getspectrum(d,nsuboct=3)
plt.plot( mids2, sums2 )
plt.xlabel( 'Freq (Hz)' )
plt.ylim( -120.0, -70.0 )
plt.xlim(20.0, 20000.0 )
plt.ylabel( 'Response' )
#plt.yscale('log')
plt.xscale('log')
plt.grid( True )
plt.show()
#freq_plot( nfreqs, spectrum, freq_step, max_freq, 'Room', True, 0.0, 20000.0, -90.0, 0.0 )