diff --git a/.gitignore b/.gitignore index c937b44..3021a61 100644 --- a/.gitignore +++ b/.gitignore @@ -18,4 +18,4 @@ test/test_workers doc LASP.egg-info lasp_octave_fir.* -lasp/aps_ui.py +lasp/ui_* diff --git a/lasp/__init__.py b/lasp/__init__.py index e96bd63..54bf4e6 100644 --- a/lasp/__init__.py +++ b/lasp/__init__.py @@ -1,3 +1,3 @@ from .wrappers import * -from .tools import * -P_REF = 2e-5 +from .lasp_tools import * +from .lasp_common import * \ No newline at end of file diff --git a/lasp/c/lasp_alg.h b/lasp/c/lasp_alg.h index 1443576..b26daf0 100644 --- a/lasp/c/lasp_alg.h +++ b/lasp/c/lasp_alg.h @@ -25,10 +25,10 @@ static inline d vd_dot(const vd * a,const vd* b) { return d_dot(getvdval(a,0),getvdval(b,0),a->n_rows); } -/** +/** * y = fac * y * - * @param y + * @param y * @param fac scale factor */ static inline void dmat_scale(dmat* y,const c fac){ @@ -38,10 +38,10 @@ static inline void dmat_scale(dmat* y,const c fac){ } } -/** +/** * y = fac * y * - * @param y + * @param y * @param fac scale factor */ static inline void cmat_scale(cmat* y,const c fac){ diff --git a/lasp/c/lasp_filterbank.c b/lasp/c/lasp_filterbank.c index 34b2c73..8ca4512 100644 --- a/lasp/c/lasp_filterbank.c +++ b/lasp/c/lasp_filterbank.c @@ -98,11 +98,13 @@ dmat FilterBank_filter(FilterBank* fb, dmat input_block = dmat_alloc(nfft,1); - /* Output is ready to be multiplied with FFt'th filter */ + /* FFT'th filter coefficients */ cmat input_fft = cmat_alloc(nfft/2+1,1); - /* Output is ready to be multiplied with FFt'th filter */ + /* Output of the fast convolution */ cmat output_fft = cmat_alloc(nfft/2+1,nfilters); + + /* Inverse FFT'th output */ dmat output_block = dmat_alloc(nfft,nfilters); while (dFifo_size(input_fifo) >= nfft) { @@ -111,6 +113,7 @@ dmat FilterBank_filter(FilterBank* fb, &input_block, fb->P_m_1 /* save P-1 samples */ ); + dbgassert(nsamples == nfft,"BUG in dFifo"); Fft_fft(fb->fft,&input_block,&input_fft); @@ -127,7 +130,9 @@ dmat FilterBank_filter(FilterBank* fb, vc_free(&output_fft_col); vc_free(&filter_col); + } + vc_free(&input_fft_col); Fft_ifft(fb->fft,&output_fft,&output_block); @@ -136,6 +141,8 @@ dmat FilterBank_filter(FilterBank* fb, fb->P_m_1,0, /* startrow, startcol */ nfft-fb->P_m_1, /* Number of rows */ output_block.n_cols); + + /* Push the valid samples to the output FIFO */ dFifo_push(fb->output_fifo,&valid_samples); dmat_free(&valid_samples); diff --git a/lasp/c/lasp_mat.h b/lasp/c/lasp_mat.h index f4e31c9..e1b8421 100644 --- a/lasp/c/lasp_mat.h +++ b/lasp/c/lasp_mat.h @@ -357,9 +357,11 @@ static inline void dmat_copy_rows(dmat* to,const dmat* from, us startrow_from, us nrows) { us col,ncols = to->n_cols; - + dbgassert(to && from,NULLPTRDEREF); + dbgassert(to->n_cols == from->n_cols,SIZEINEQUAL); dbgassert(startrow_from+nrows <= from->n_rows,OUTOFBOUNDSMATR); dbgassert(startrow_to+nrows <= to->n_rows,OUTOFBOUNDSMATR); + for(col=0;col self._Lmax: + self._Lmax = curmax + + return Level diff --git a/lasp/lasp_tools.py b/lasp/lasp_tools.py new file mode 100644 index 0000000..0113368 --- /dev/null +++ b/lasp/lasp_tools.py @@ -0,0 +1,8 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +"""! +Author: J.A. de Jong - ASCEE + +Description: +""" +import numpy as np diff --git a/lasp/lasp_weighcal.py b/lasp/lasp_weighcal.py new file mode 100644 index 0000000..b174dda --- /dev/null +++ b/lasp/lasp_weighcal.py @@ -0,0 +1,136 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Weighting and calibration filter in one +@author: J.A. de Jong - ASCEE +""" +from .filter_design.freqweighting_fir import A,C +from .lasp_common import FreqWeighting +from .filter_design.fir_design import (arbitrary_fir_design, + freqResponse as frp) +from lasp.lasp_config import ones, empty +from .wrappers import FilterBank +import numpy as np + +__all__ = ['WeighCal'] + + +class WeighCal: + """ + Time weighting and calibration FIR filter + """ + def __init__(self, fw = FreqWeighting.default, + nchannels = 1, + fs = 48000., + calfile = None, + sens=1.0): + """ + Initialize the frequency weighting and calibration FIR filters. + + Args: + fw: Frequency weighting to apply + nchannels: Number of channels for the input data + fs: Sampling frequency [Hz] + calfile: Calibration file to load. + sens: Sensitivity in units [\f$ Pa^{-1} \f$] + """ + + self.nchannels = nchannels + self.fs = fs + self.fw = fw + self.sens = sens + self.calfile = calfile + + # Frequencies used for the filter design + freq_design = np.linspace(0,17e3,3000) + freq_design[-1] = fs/2 + + # Objective function for the frequency response + frp_obj = self.frpObj(freq_design) + + self._firs = [] + self._fbs = [] + for chan in range(self.nchannels): + fir = arbitrary_fir_design(fs,2048,freq_design, + frp_obj[:,chan], + window='rectangular') + self._firs.append(fir) + + self._fbs.append(FilterBank(fir[:,np.newaxis],4096)) + + self._freq_design = freq_design + + def filter_(self,data): + """ + Filter data using the calibration and frequency weighting filter. + """ + nchan = self.nchannels + assert data.shape[1] == nchan + assert data.shape[0] > 0 + + filtered = [] + for chan in range(nchan): + filtered.append(self._fbs[chan].filter_(data[:,chan])[:,0]) + filtered = np.asarray(filtered).transpose()/self.sens + if filtered.ndim == 1: + filtered = filtered[:,np.newaxis] + return filtered + + + def frpCalObj(self, freq_design): + """ + Computes the objective frequency response of the calibration filter + """ + calfile = self.calfile + if calfile is not None: + cal = np.loadtxt(calfile,skiprows=2) + freq = cal[:,0] + cal = cal[:,1:] + if cal.shape[1] != self.nchannels: + raise ValueError('Number of channels in calibration file does' + ' not equal to given number of channels') + calfac = 10**(-cal/20) + filter_calfac = empty((freq_design.shape[0],self.nchannels)) + + for chan in range(self.nchannels): + filter_calfac[:,chan] = np.interp(freq_design,freq,calfac[:,chan]) + + else: + filter_calfac = ones((freq_design.shape[0],self.nchannels,)) + + return filter_calfac + + def frpWeightingObj(self, freq_design): + """ + Computes the objective frequency response of the frequency weighting + filter. + """ + fw = self.fw + if fw == FreqWeighting.A: + return A(freq_design) + elif fw == FreqWeighting.C: + return C(freq_design) + elif fw == FreqWeighting.Z: + return ones(freq_design.shape[0]) + else: + raise ValueError('Invalid fw parameter') + + def frpObj(self, freq_design): + """ + Combines the frequency weighting and the calibration filter into + one frequency response objective function. + """ + # Objective function for the frequency response + frp_objective = self.frpCalObj(freq_design) * \ + self.frpWeightingObj(freq_design)[:,np.newaxis] + frp_objective[-1] = 0. + + return frp_objective + + def freqResponse(self, chan=0, freq=None): + """ + Returns the frequency response of the designed FIR filter + """ + if freq is None: + freq = np.logspace(1,np.log10(self.fs/2),500) + return freq,frp(self.fs,freq,self._firs[chan]),self.frpObj(freq)[:,chan] diff --git a/lasp/tools.py b/lasp/tools.py deleted file mode 100644 index fbe8729..0000000 --- a/lasp/tools.py +++ /dev/null @@ -1,13 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -"""! -Author: J.A. de Jong - ASCEE - -Description: -""" -import numpy as np - -def getFreq(fs, nfft): - df = fs/nfft # frequency resolution - K = nfft//2+1 # number of frequency bins - return np.linspace(0, (K-1)*df, K) diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx index f59c68f..1ba6aca 100644 --- a/lasp/wrappers.pyx +++ b/lasp/wrappers.pyx @@ -98,8 +98,6 @@ cdef class Fft: return timedata - - cdef extern from "lasp_window.h": ctypedef enum WindowType: Hann @@ -214,7 +212,7 @@ cdef class AvPowerSpectra: cdef vd weighting_vd cdef vd* weighting_ptr = NULL if(weighting.size != 0): - weighting_vd = dmat_foreign_data(weighting.n_rows,1, + weighting_vd = dmat_foreign_data(weighting.size,1, &weighting[0],False) weighting_ptr = &weighting_vd @@ -295,7 +293,7 @@ cdef class FilterBank: cdef dmat hmat = dmat_foreign_data(h.shape[0], h.shape[1], &h[0,0], - True) + False) self.fb = FilterBank_create(&hmat,nfft) dmat_free(&hmat) @@ -376,6 +374,9 @@ cdef class SPLowpass: SPLowpass_free(self.lp) def filter_(self,d[:] input_): + if input_.shape[0] == 0: + return np.array([],dtype=NUMPY_FLOAT_TYPE) + cdef vd input_vd = dmat_foreign_data(input_.shape[0],1, &input_[0],False)