From 96edc90c55b355328810a013f1d5178aeaa6849b Mon Sep 17 00:00:00 2001
From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" <j.a.dejong@ascee.nl>
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
+