Improved bar colors. Ticks on y axis in float on bars. Legend positions in bars. Moved asceefigs repo into lasp.tools.

This commit is contained in:
Anne de Jong 2018-10-01 14:25:09 +02:00 committed by J.A. de Jong - ASCEE
parent 95f297dd46
commit c152bf65f3
10 changed files with 669 additions and 8 deletions

View File

@ -7,7 +7,16 @@ Common definitions used throughout the code.
"""
__all__ = ['P_REF', 'FreqWeighting', 'TimeWeighting', 'getTime', 'calfile',
'W_REF']
'W_REF', 'DEFAULT_FIGSIZE_H', 'DEFAULT_FIGSIZE_W',
'GOLDEN', 'PLOT_COLORS_LIST', 'PLOT_NOCOLORS_LIST']
PLOT_COLORS_LIST = ['b', 'g', 'r', 'c', 'm', 'y', 'k', '#BE6500']
PLOT_NOCOLORS_LIST = ['k', '0.5', 'k', '0.5', '0.5', '0.5', '0.5', 'k']
DEFAULT_FIGSIZE_W = 8
GOLDEN = (np.sqrt(5.)+1)/2
DEFAULT_FIGSIZE_H = DEFAULT_FIGSIZE_W/GOLDEN
# Reference sound pressure level
P_REF = 2e-5

View File

@ -42,6 +42,7 @@ import numpy as np
from .lasp_config import LASP_NUMPY_FLOAT_TYPE
import wave
import os
import numbers
class BlockIter:
@ -106,8 +107,11 @@ def scaleBlockSens(block, sens):
"""
assert sens.ndim == 1
assert sens.size == block.shape[1]
sw = getSampWidth(block.dtype)
fac = 2**(8*sw - 1) - 1
if isinstance(block.dtype, numbers.Integral):
sw = getSampWidth(block.dtype)
fac = 2**(8*sw - 1) - 1
else:
fac = 1,
return block.astype(LASP_NUMPY_FLOAT_TYPE)/fac/sens[np.newaxis, :]

View File

@ -8,7 +8,7 @@ Description:
Class for plotting bars on a QGraphicsScene.
"""
from ..lasp_gui_tools import ASCEEColors
from ..lasp_gui_tools import ASCEEColors, Branding
from PySide.QtGui import (
QGraphicsScene, QPen, QBrush, QGraphicsRectItem,
QGraphicsTextItem, QPainter, QImage, QPrinter
@ -33,7 +33,9 @@ ticklength = 10
# Distance between two bar groups in units of bar thicknesses
dxbars = 2
DEFAULT_COLORS = [ASCEEColors.blue, ASCEEColors.green, Qt.red, Qt.yellow]
DEFAULT_COLORS = [ASCEEColors.blue, ASCEEColors.green, Qt.red, Qt.cyan,
Qt.darkYellow,
Qt.darkMagenta]
class BarScene(QGraphicsScene):
@ -109,7 +111,7 @@ class BarScene(QGraphicsScene):
range_ = ylim[1]-ylim[0]
ytickval = i/(nyticks-1)*range_ + ylim[0]
yticklabel = f'{ytickval:3.3}'
yticklabel = f'{ytickval:.0f}'
txt = QGraphicsTextItem(yticklabel)
txtwidth = txt.boundingRect().width()
txtmaxwidth = max(txtmaxwidth, txtwidth)
@ -172,8 +174,7 @@ class BarScene(QGraphicsScene):
if legend is not None:
maxlegtxtwidth = 0
legposx = 0 if legendpos is None else legendpos[0]
legposy = 0 if legendpos is None else legendpos[1]
legposx, legposy = (0,0) if legendpos is None else legendpos
legpos = (xsize-rightoffset-300+legposx,
ysize-topoffset-30+legposy)

7
lasp/tools/__init__.py Normal file
View File

@ -0,0 +1,7 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from .config import report_quality
from .report_tools import *
from lasp.lasp_common import FreqWeighting, TimeWeighting
__all__ = ['report_quality', 'PSPlot', 'LevelBars', 'Levels', 'close',
'FreqWeighting', 'TimeWeighting']

44
lasp/tools/config.py Normal file
View File

@ -0,0 +1,44 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""!
Author: J.A. de Jong - ASCEE
Description:
"""
__all__ = ['report_quality', 'getReportQuality']
_report_quality = False
def getReportQuality():
global _report_quality
return _report_quality
def report_quality():
import matplotlib as mpl
# mpl.use('Qt5Agg')
global _report_quality
# mpl.use("pgf")
# rc('font',**{'family':'serif','serif':['Libertine']})
# for Palatino and other serif fonts use:
# rc('font',**{'family':'serif','serif':['Palatino']})
# rc('text', usetex=True)
# TeX preamble
preamble = [
r'\usepackage{libertine-type1}'
r'\usepackage[libertine]{newtxmath}'
# r'\usepackage{fontspec}',
# r'\setmainfont{Libertine}',
]
params = {
'font.family': 'serif',
'text.usetex': True,
'text.latex.unicode': True,
'pgf.rcfonts': False,
'pgf.texsystem': 'pdflatex',
'pgf.preamble': preamble,
}
mpl.rcParams.update(params)
_report_quality = True

214
lasp/tools/plot.py Normal file
View File

@ -0,0 +1,214 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""!
Author: J.A. de Jong - ASCEE
Description:
"""
__all__ = ['close', 'Figure', 'Bode', 'PS', 'PSD']
from .config import getReportQuality
import matplotlib.pyplot as plt
import numpy as np
from cycler import cycler
from lasp.lasp_common import (PLOT_COLORS_LIST, PLOT_NOCOLORS_LIST,
DEFAULT_FIGSIZE_H, DEFAULT_FIGSIZE_W)
def close():
plt.close('all')
class Figure:
def __init__(self, **kwargs):
ncols = kwargs.pop('ncols', 1)
nrows = kwargs.pop('nrows', 1)
color = kwargs.pop('color', (PLOT_NOCOLORS_LIST if getReportQuality()
else PLOT_COLORS_LIST))
if isinstance(color, bool):
if color:
color = PLOT_COLORS_LIST
else:
color = PLOT_NOCOLORS_LIST
colors = cycler('color', color)
figsize = kwargs.pop('figsize', (DEFAULT_FIGSIZE_W, DEFAULT_FIGSIZE_H))
self._f = plt.figure(figsize=figsize)
marker = kwargs.pop('marker', False)
if marker:
markers = cycler(marker=['o', 's', 'D', 'X', 'v', '^', '<', '>'])
else:
markers = cycler(marker=[None]*8)
linewidths = cycler(linewidth=[1, 2, 1, 2, 2, 3, 2, 1])
linestyles = cycler(linestyle=['-', '-', '--', ':', '-', '--', ':', '-.', ])
self._ax = []
self._legend = {}
for row in range(nrows):
self._legend[row] = {}
for col in range(ncols):
self._legend[row][col] = []
ax = self._f.add_subplot(100*nrows
+ 10*ncols
+ (row*ncols + col)+1)
ax.set_prop_cycle(
colors+linestyles+markers+linewidths)
self._ax.append(ax)
self._ncols = ncols
self._cur_ax = self._ax[0]
self._cur_col = 0
self._cur_row = 0
self._zorder = -1
def setAx(self, row, col):
self._cur_ax = self._ax[row*self._ncols+col]
@property
def fig(self):
return self._f
def markup(self):
for ax in self._ax:
ax.grid(True, 'both')
self._zorder -= 1
if not getReportQuality():
self._f.show()
def vline(self, x):
self._ax[0].axvline(x)
def plot(self, *args, **kwargs):
line = self._cur_ax.plot(*args, **kwargs, zorder=self._zorder)
self.markup()
return line
def loglog(self, *args, **kwargs):
line = self._cur_ax.loglog(*args, **kwargs, zorder=self._zorder)
self.markup()
return line
def semilogx(self, *args, **kwargs):
line = self._cur_ax.semilogx(*args, **kwargs, zorder=self._zorder)
self.markup()
return line
def xlabel(self, *args, **kwargs):
all_ax = kwargs.pop('all_ax', False)
if all_ax:
for ax in self._ax:
ax.set_xlabel(*args, **kwargs)
else:
self._cur_ax.set_xlabel(*args, **kwargs)
def ylabel(self, *args, **kwargs):
all_ax = kwargs.pop('all_ax', False)
if all_ax:
for ax in self._ax:
ax.set_ylabel(*args, **kwargs)
else:
self._cur_ax.set_ylabel(*args, **kwargs)
def legend(self, leg, *args, **kwargs):
# all_ax = kwargs.pop('all_ax', False)
if isinstance(leg, list) or isinstance(leg, tuple):
self._legend[self._cur_col][self._cur_col] = list(leg)
else:
self._legend[self._cur_col][self._cur_col].append(leg)
self._cur_ax.legend(self._legend[self._cur_col][self._cur_col])
def savefig(self, *args, **kwargs):
self._f.savefig(*args, **kwargs)
def xlim(self, *args, **kwargs):
all_ax = kwargs.pop('all_ax', False)
if all_ax:
for ax in self._ax:
ax.set_xlim(*args, **kwargs)
else:
self._cur_ax.set_xlim(*args, **kwargs)
def ylim(self, *args, **kwargs):
all_ax = kwargs.pop('all_ax', False)
if all_ax:
for ax in self._ax:
ax.set_ylim(*args, **kwargs)
else:
self._cur_ax.set_ylim(*args, **kwargs)
def title(self, *args, **kwargs):
self._cur_ax.set_title(*args, **kwargs)
def xticks(self, ticks):
for ax in self._ax:
ax.set_xticks(ticks)
def close(self):
plt.close(self._f)
def xscale(self, scale):
for ax in self._ax:
ax.set_xscale(scale)
class Bode(Figure):
def __init__(self, *args, **kwargs):
super().__init__(naxes=2, *args, **kwargs)
def add(self, freq, phasor, qtyname='G', **kwargs):
L = 20*np.log10(np.abs(phasor))
phase = np.angle(phasor)*180/np.pi
self.semilogx(freq, L, axno=0, **kwargs)
self.semilogx(freq, phase, axno=1, **kwargs)
self.ylabel('$L$ [%s] [dB]' % qtyname, axno=0)
self.ylabel(fr'$\angle$ {qtyname} [$^\circ$]', axno=1)
self.xlabel('Frequency [Hz]', axno=1)
class PS(Figure):
def __init__(self, ref, *args, **kwargs):
super().__init__(naxes=1, *args, **kwargs)
self.ref = ref
def add(self, fs, freq, ps, qtyname='C', **kwargs):
overall = np.sum(ps)
print(overall)
overall_db = 10*np.log10(overall/self.ref**2)
L = 10*np.log10(np.abs(ps)/self.ref**2)
self.semilogx(freq, L, **kwargs)
# self.plot(freq,L,**kwargs)
self.ylabel('Level [dB re 20$\\mu$Pa]')
self.xlabel('Frequency [Hz]')
self.legend('%s. Overall SPL = %0.1f dB SPL' % (qtyname, overall_db))
class PSD(PS):
def __init__(self, ref, *args, **kwargs):
"""
Initialize a PSD plot
Args:
ref: Reference value for level in dB's
"""
super().__init__(ref, *args, **kwargs)
def add(self, fs, freq, ps, qtyname='C', **kwargs):
df = freq[1]-freq[0]
nfft = fs/df
df = fs/nfft
psd = ps / df
overall = np.sum(np.abs(ps), axis=0)
overall_db = 10*np.log10(overall/self.ref**2)
L = 10*np.log10(abs(psd)/self.ref**2)
self.semilogx(freq, L, **kwargs)
self.ylabel('$L$ [%s] [dB re %0.0e]' % (qtyname, self.ref))
self.xlabel('Frequency [Hz]')
self.legend('%s. Overall SPL = %0.1f dB SPL' % (qtyname, overall_db))

163
lasp/tools/report_tools.py Normal file
View File

@ -0,0 +1,163 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Author: J.A. de Jong - ASCEE
Description: backend tools for easy postprocessing of measurements
"""
from .plot import Figure
from lasp.wrappers import AvPowerSpectra
from lasp.lasp_measurement import Measurement
from lasp.lasp_common import FreqWeighting, TimeWeighting, getFreq, getTime
from lasp.lasp_weighcal import WeighCal
from lasp.lasp_octavefilter import OctaveFilterBank, ThirdOctaveFilterBank
from lasp.lasp_barfigure import BarFigure
from lasp.lasp_figure import Plotable
from lasp.lasp_figure import Plotable, PlotOptions
from lasp.lasp_slm import SLM
import numpy as np
import sys
def close():
import matplotlib.pyplot as plt
plt.close('all')
def PSPlot(fn_list, **kwargs):
"""
Create a power spectral density plot, ASCEE style
Args:
fn_list: list of measurement filenames to plot PSD for
fw:
fs:
nfft:
xscale:
yscale:
"""
fw = kwargs.pop('fw', FreqWeighting.A)
nfft = kwargs.pop('nfft', 2048)
xscale = kwargs.pop('xscale', 'log')
yscale = kwargs.pop('yscale', 'PSD')
ylim = kwargs.pop('ylim', (0, 100))
xlim = kwargs.pop('xlim', (100, 10000))
f = Figure(**kwargs)
print(kwargs)
if xscale == 'log':
pltfun = f.semilogx
else:
pltfun = f.plot
for fn in fn_list:
meas = Measurement(fn)
fs = meas.samplerate
data = meas.praw()
aps = AvPowerSpectra(nfft, 1, 50.)
weighcal = WeighCal(fw, nchannels=1,
fs=fs)
weighted = weighcal.filter_(data)
ps = aps.addTimeData(weighted)
freq = getFreq(fs, nfft)
if yscale == 'PSD':
df = fs/nfft
type_str = '/$\\sqrt{\\mathrm{Hz}}$'
elif yscale == 'PS':
df = 1.
type_str = ''
else:
raise ValueError("'type' should be either 'PS' or 'PSD'")
psd_log = 10*np.log10(ps[:, 0, 0].real/df/2e-5**2)
pltfun(freq, psd_log)
f.xlabel('Frequency [Hz]')
f.ylabel(f'Level [dB({fw[0]}) re (20$\\mu$ Pa){type_str}')
f.ylim(ylim)
f.xlim(xlim)
return f
def Levels(fn_list, **kwargs):
bank = kwargs.pop('bank', 'third')
fw = kwargs.pop('fw', FreqWeighting.A)
tw = kwargs.pop('tw', TimeWeighting.fast)
xmin_txt = kwargs.pop('xmin', '100')
xmax_txt = kwargs.pop('xmax', '16k')
levels = []
leveltype = 'eq'
for fn in fn_list:
meas = Measurement(fn)
fs = meas.samplerate
weighcal = WeighCal(fw, nchannels=1,
fs=fs, calfile=None)
praw = meas.praw()
weighted = weighcal.filter_(praw)
if bank == 'third':
filt = ThirdOctaveFilterBank(fs)
xmin = filt.nominal_txt_tox(xmin_txt)
xmax = filt.nominal_txt_tox(xmax_txt)
elif bank == 'overall':
slm = SLM(meas.samplerate, tw)
slm.addData(weighted)
levels.append(
Plotable(' ', slm.Lmax if leveltype == 'max' else slm.Leq,
name=meas.name))
continue
else:
raise NotImplementedError()
# Octave bands
# filt = OctaveFilterBank(fs)
filtered_out = filt.filter_(weighted)
level = np.empty((xmax - xmin + 1))
xlabels = []
for i, x in enumerate(range(xmin, xmax+1)):
nom = filt.nominal_txt(x)
xlabels.append(nom)
filt_x = filtered_out[nom]['data']
slm = SLM(filt.fs, tw)
slm.addData(filt_x)
leveli = slm.Lmax if leveltype == 'max' else slm.Leq
level[i] = leveli
levels.append(Plotable(xlabels, level, name=meas.name))
return levels
def LevelBars(levels, show=True, **kwargs):
from PySide.QtGui import QApplication, QFont
app = QApplication.instance()
if not app:
app = QApplication(sys.argv)
app.setFont(QFont('Linux Libertine'))
size = kwargs.pop('size', (1200, 600))
opts = PlotOptions.forLevelBars()
opts.ylim = kwargs.pop('ylim', (0, 100))
opts.ylabel = kwargs.pop('ylabel', 'LAeq [dB(A)]')
opts.xlabel = kwargs.pop('xlabel', None)
opts.legend = kwargs.pop('legend', [level.name for level in levels])
opts.legendpos = kwargs.pop('legendpos', None)
opts.title = kwargs.pop('title', None)
def BarPlotter(ptas, pto):
fig = BarFigure(None, pto)
fig.resize(1200, 300)
for pta in ptas:
fig.add(pta)
return fig
fig = BarPlotter(levels, opts)
if show:
fig.show()
fig.resize(*size)
if show:
app.exec_()
return fig

200
scripts/lasp_postprocess.py Executable file
View File

@ -0,0 +1,200 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""!
Author: J.A. de Jong - ASCEE
Description:
"""
import matplotlib as mpl
mpl.use('Qt5Agg')
preamble = [
r'\usepackage{libertine-type1}'
r'\usepackage[libertine]{newtxmath}'
# r'\usepackage{fontspec}',
#r'\setmainfont{Libertine}',
]
params = {
'font.family': 'serif',
'text.usetex': True,
'text.latex.unicode': True,
'pgf.rcfonts': False,
'pgf.texsystem': 'pdflatex',
'pgf.preamble': preamble,
}
mpl.rcParams.update(params)
from asceefigs.config import report_quality
# report_quality()
from lasp.lasp_measurement import Measurement
from lasp.lasp_weighcal import WeighCal, FreqWeighting
from lasp.lasp_slm import SLM, TimeWeighting, P_REF
from lasp.wrappers import AvPowerSpectra
from lasp.lasp_common import getFreq
import os
import sys
import numpy as np
from asceefigs.plot import Figure, PS, close
import matplotlib.pyplot as pl
import time
#close('all')
def select_file():
measfiles = []
for path, name, filenames in os.walk(os.getcwd()):
for filename in filenames:
if '.h5' in filename:
measfiles.append(os.path.join(path, filename))
if len(measfiles) == 0:
raise ValueError('No measurement files in current directory')
for i, mfile in enumerate(measfiles):
print('- %20s: %i' % (mfile, i))
no = input('Please select a file (%i - %i):' % (0, len(measfiles)))
1
no = int(no)
return measfiles[no]
fn = sys.argv[2]
# type_ : spec, sp
type_ = sys.argv[1]
#else:
# fn = select_file()
##
#fn = 'Cees_1'
#fn = 'holland_1'
fs = 48000
meas = Measurement(fn)
ts = meas.time
print('Measurement time: ', time.ctime(ts))
praw = meas.praw()[:, [0]]
N = praw.shape[0]
if type_ == 'spec':
weighcal_slmA = WeighCal(FreqWeighting.A)
weighcal_slmC = WeighCal(FreqWeighting.C)
filtered_A = weighcal_slmA.filter_(praw)
filtered_C = weighcal_slmC.filter_(praw)
# %% Sound level meter
tw = TimeWeighting.fast
slm_A = SLM(fs, tw=tw)
slm_C = SLM(fs, tw=tw)
LAF = slm_A.addData(filtered_A)
LCF = slm_C.addData(filtered_C)
# Strip off filter initialization part
LAF = LAF[int(tw[0]*fs):, :]
LCF = LCF[int(tw[0]*fs):, :]
N = LAF.shape[0]
t = np.linspace(0, N/fs, N, False)
Lfig = pl.figure(figsize=(8, 8))
# from matplotlib import gridspec
#gs = gridspec.GridSpec(2, 2, height_ratios=(3, 1), width_ratios=(9, 1))
#Lax = pl.subplot(gs[0,0])
ax = Lfig.add_subplot(211)
pl.title('')
pl.plot(t, LAF[:, 0])
pl.plot(t, LCF[:, 0])
pl.ylabel('Level [dB]')
pl.legend(['LAF', 'LCF'])
pl.ylim(40, 80)
pl.xlim(0, meas.T)
pl.grid('on', 'both')
print('Maximum A-level:', slm_A.Lmax)
print('Maximum C-level:', slm_C.Lmax)
from scipy.signal import spectrogram
nfft = 8192
noverlap = 4*nfft//5
freq_sp, t_sp, Sxx = spectrogram(filtered_C[:, 0], fs, scaling='spectrum',
window=(
'hann'), nperseg=nfft, noverlap=noverlap)
# Chop off higher frequencies
nfreq = freq_sp.shape[0]
nf_start = 5
nf_end = int(nfreq/40)
freq_sp = freq_sp[nf_start:nf_end]
Sxx = Sxx[nf_start:nf_end, :]
Sxxpow = 10*np.log10(Sxx/P_REF**2)
maxSxx = np.max(Sxxpow)
minSxx = np.min(Sxxpow)
print(f'max Sxx: {maxSxx:.1f} dB')
print(f'min Sxx: {minSxx:.1f} dB')
#SPax = pl.subplot(gs[1,0])
ax1 = Lfig.add_subplot(212)
pl.title('C-weighted spectrogram')
colormesh = pl.pcolormesh(t_sp, freq_sp, Sxxpow,
cmap='rainbow', vmin=30, vmax=60)
pl.xlabel('Time [s]')
pl.ylabel('Frequency [Hz]')
pl.xlim(0, meas.T)
pl.yscale('linear')
# %% Colorbar
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
axins = inset_axes(ax1,
width="5%", # width = 10% of parent_bbox width
height="100%", # height : 50%
loc=6,
bbox_to_anchor=(1.02, 0., 1, 1),
bbox_transform=ax1.transAxes,
borderpad=0,
)
#pl.colorbar(orientation='horizontal')
#cax = pl.subplot(gs[1,1])
cb = pl.colorbar(colormesh, cax=axins)
fn_base = os.path.splitext(fn)[0]
# pl.savefig('%s.eps' %fn,dpi=600)
pl.show(block=False)
#pl.subplot2grid((2,2),(1,1))
#pl.tight_layout()
#l1 = pl.axvline(2,color='k')
#pl.axvline(3,color='k')
#del l1
# %%
else:
# %%
weighcal_A = WeighCal(FreqWeighting.A,
calfile='/home/anne/wip/UMIK-1/cal/7027430_90deg.txt',
sens=0.053690387255872614)
filtered_A = weighcal_A.filter_(praw)
nfft = 8192
aps = AvPowerSpectra(nfft, 1, 50.)
ps = aps.addTimeData(filtered_A)
psplot = PS(P_REF)
freq = getFreq(fs, nfft)
psplot.add(fs, freq, ps[:, 0, 0], 'A-gewogen')
psplot.xlim(10, 10000)
psplot.ylim(20, 80)
psplot.savefig('%s_sp.eps' % fn)
pl.show()
# %%
input('Press any key to close')
pl.savefig(f'{fn_base}_spectrogram.png', dpi=600)
pl.savefig(f'{fn_base}_spectrogram.eps')

13
scripts/lasp_towav Executable file
View File

@ -0,0 +1,13 @@
#!/usr/bin/python
from lasp.lasp_measurement import Measurement
import argparse, os
parser = argparse.ArgumentParser(
description='Playback recorded measurement'
)
parser.add_argument('filename', help='Filename of measurement', type=str)
args = parser.parse_args()
meas = Measurement(args.filename)
fnwav = os.path.splitext(args.filename)[0]
meas.exportAsWave(fnwav, force=True)

6
tools/__init__.py Normal file
View File

@ -0,0 +1,6 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from .config import report_quality
from .report_tools import *
__all__ = ['report_quality']