lasp/scripts/lasp_postprocess.py

200 lines
5.1 KiB
Python
Raw Normal View History

#!/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')