Numpy: calculate Power Spectral Density for two test files

#!/usr/bin/python
import numpy as np
import scipy.io.wavfile
import scipy.signal
import matplotlib.pyplot as plt

pwr_to_dB_fct = 10 / np.log(10)
pwr_to_dB = lambda x: np.log(x)*pwr_to_dB_fct


fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_title('Difference 48kHz - 96kHz')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Power Density [1/Hz]')
ax.set_xlim(400, 40e3)
ax.grid()

for fn in ['dynamic-test-ospac.wav', 'dynamic-highpass.wav'] :
    fS, wave = scipy.io.wavfile.read(fn)
    print('%s: %d samples at %f Hz'%(fn, wave.shape[0], fS))
    freq, Plin = scipy.signal.welch(wave, fS)
    PdB = pwr_to_dB(Plin)

    # "highpass" file is artificially amplified, we correct by
    # substracting the difference between the highest two frequencies
    # in test-ospac and highpass
    if fn == 'dynamic-highpass.wav' :
        PdB -= 46.43541

    print(freq[-5:], PdB[-5])
    ax.semilogx(freq, PdB, label=fn[:-4])

ax.legend()
fig.savefig('plot.png')