From 96edc90c55b355328810a013f1d5178aeaa6849b Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Thu, 25 Feb 2021 12:13:05 +0100 Subject: [PATCH] Implementation of impedance tube code --- lasp/lasp_imptube.py | 148 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 148 insertions(+) create mode 100644 lasp/lasp_imptube.py diff --git a/lasp/lasp_imptube.py b/lasp/lasp_imptube.py new file mode 100644 index 0000000..cda9725 --- /dev/null +++ b/lasp/lasp_imptube.py @@ -0,0 +1,148 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +"""! +Author: J.A. de Jong - ASCEE + +Description: Two-microphone impedance tube methods +""" +__all__ = ['TwoMicImpedanceTube'] +from lrftubes import Air +from .lasp_measurement import Measurement +from numpy import pi, sqrt, exp +import numpy as np +from scipy.interpolate import UnivariateSpline + +class TwoMicImpedanceTube: + def __init__(self, mnormal: Measurement, + mswitched: Measurement, + s: float, + d1: float, + d2: float, + fl: float = None, + fu: float = None, + periodic_method=False, + mat= Air(), + **kwargs): + + """ + Initialize two-microphone impedance tube methods + + Args: + mnormal: Measurement in normal configuration + mswitched: Measurement in normal configuration + s: Microphone distance + fl: Lower evaluation frequency + fu: Lower evaluation frequency + + kwargs: tuple with extra arguments, of which: + N: Period length of periodic excitation *obligatory* + chan0: Measurement channel index of mic 0 + chan1: Measurement channel index of mic 1 + + """ + self.mnormal = mnormal + self.mswitched = mswitched + + self.mat = mat + self.s = s + self.d1 = d1 + self.d2 = d2 + + self.fl = fl + if fl is None: + ksmin = 0.1*pi + kmin = ksmin/s + self.fl = kmin*mat.c0/2/pi + + self.fu = fu + if fu is None: + ksmax = 0.8*pi + kmax = ksmax/s + self.fu = kmax*mat.c0/2/pi + + self.periodic_method = periodic_method + self.channels = [kwargs.pop('chan0', 0), kwargs.pop('chan1', 1)] + # Compute calibration correction + if periodic_method: + self.N = kwargs.pop('N') + freq, C1 = mnormal.periodicCPS(self.N, channels=self.channels) + freq, C2 = mswitched.periodicCPS(self.N, channels=self.channels) + else: + self.nfft = kwargs.pop('nfft', 16000) + freq, C1 = mnormal.CPS(nfft=self.nfft, channels=self.channels) + freq, C2 = mswitched.CPS(nfft=self.nfft, channels=self.channels) + + # Store this, as it is often used for just a single sample. + self.C1 = C1 + self.freq = freq + + self.il = np.where(self.freq<= self.fl)[0][-1] + self.ul = np.where(self.freq > self.fu)[0][0] + + # Calibration correction factor + # self.K = 0*self.freq + 1.0 + K = sqrt(C2[:,0,1]*C1[:,0,0]/(C2[:,1,1]*C1[:,1,0])) + self.K = UnivariateSpline(self.freq, K.real)(self.freq) +\ + 1j*UnivariateSpline(self.freq, K.imag)(self.freq) + + def cut_to_limits(self, ar): + return ar[self.il:self.ul] + + def getFreq(self): + """ + Array of frequencies, cut to limits of validity + """ + return self.cut_to_limits(self.freq) + + def G_AB(self, meas): + if meas is self.mnormal: + C = self.C1 + freq = self.freq + else: + if self.periodic_method: + freq, C = meas.periodicCPS(self.N, self.channels) + else: + freq, C = meas.CPS(nfft=self.nfft, channels=self.channels) + + # Microphone transfer function + G_AB = self.K*C[:,1,0]/C[:,0,0] + return self.getFreq(), self.cut_to_limits(G_AB) + + def R(self, meas): + freq, G_AB = self.G_AB(meas) + s = self.s + k = 2*pi*freq/self.mat.c0 + d1 = self.d1 + RpA = (G_AB - exp(-1j*k*s))/(exp(1j*k*s)-G_AB) + + R = RpA*exp(2*1j*k*(s+d1)) + + return freq, R + + def alpha(self, meas): + """ + Acoustic absorption coefficient + """ + freq, R = self.R(meas) + return freq, 1 - np.abs(R)**2 + + def z(self, meas): + """ + Acoustic impedance at the position of the sample, in front of the sample + """ + freq, R = self.R(meas) + return freq, self.mat.z0*(1+R)/(1-R) + + def zs(self, meas): + """ + Sample impedance jump, assuming a cavity behind the sample with + thickness d2 + """ + freq, R = self.R(meas) + z0 = self.mat.z0 + k = 2*pi*freq/self.mat.c0 + d2 = self.d2 + + zs = 2*z0*(1-R*exp(2*1j*k*d2))/((R-1)*(exp(2*d2*k*1j)-1)) + return freq, zs +