# PSEUDO CODE FOR MICROPHONE MEASUREMENTS USING MUTLTI MICROPHONE METHOD

* p-t series data at the microphone locations
- export the lvm data into matlab using lvm import file
<https://www.mathworks.com/matlabcentral/fileexchange/19913-lvm-file-import>
- get the p-t series data at the microphone locations

* p-f data i.e, pressure frequency domain data at microphone locations
- the pressure-time series data is converted to pressure frequency domain using fft (fast fourier transform) command in matlab

* matlab code for extracing the lvm data
- extract the lvm file into matlab and get p-t series, after getting p-t series we need to convert into p-f domain using fft and this code also contains to find the distribution of phase angle along the time
---

    clear
    close all
    Fs = 8192; %sampling frequency
    data_struct = lvm_import('test_4230.lvm'); %test lvm file
   data = data_struct.Segment1.data(:,1); %extracting the first column of the data
   time = 0:100/819199:100; % time
   v = data'; % voltage 
   p = 1000*v; % pressure in pa
   plot(time,p) % plotting of pressure-time series
   xlabel('time (s)');
   ylabel('Pressure');
   title('pressure time series data');

   %FFT
   j = length(p);
   k = 2^nextpow2(j);
   freq = Fs*(0:k/2-1)/k; % frequency range
   ff = fft(p,k);   % doing fft (p-f domain)
   z= ff(1:k/2);
   a = abs(z);
   plot(freq,a) % graph
   xlabel('Frequency (Hz)');
   ylabel('Fourier spectrum pressure ampitude');
   title('FFT presure frequency domain');
   
   % code for phase distribution
   p = angle(ff)*180/pi; % angle in degrees
   figure,plot(k,p); % plotting of phase distribution
   xlabel('frequency');
   ylabel('phase angle')
   

---

cutoff frequency of the tube for the first order modes = 1596Hz
($\phi=140mm$)

* after doing fft of the data we obtained, we check the values of pressure amplitudes which has more amplitude i.e, at a particular frequency (frequency is nearer to the cutoff frequency) and obtain the pressures at locations 1,2,3,4,5 (absolutely calibrated microphone is placed at 5th location) i.e, at locations 1,2,3,4,5 pressures $p_1,p_2,p_3,p_4,p_5$

* using MMM (multi microphone method)

$$ p(f) = p_+(f)exp(-ikx) + p_-(f)exp(ikx) $$


the over determined linear system of equations can be written as

$$Ax=b$$

* where

$$ A=\begin{bmatrix} exp(-ikx_1)&exp(ikx_1) \\ exp(-ikx_2)&exp(ikx_2)\\\vdots&\vdots\\exp(-ikx_5)&exp(ikx_5) \end{bmatrix} $$

$$ x=\begin{bmatrix} p_+(f) \\ p_-(f) \end{bmatrix} $$

$$ b=\begin{bmatrix} p_1(f)\\p_2(f)\\\vdots&\\p_5(f) \end{bmatrix} $$

* as $x_1,x_2,x_3,x_4,x_5$ and $p_1,p_2,p_3,p_4,p_5$ are already known, we can find $p_+(f)$ and $p_-(f)$ i.e, forward and backward wave propogation constants

* using $x=pinv(A)*b$
where $pinv(A)$ = pseudo inverse of matrix A

* after getting $p_+(f)$ and $p_-(f)$ subtituing in the pressure variation equation i.e, $$ p(f) = p_+(f)exp(-ikx) + p_-(f)exp(ikx) $$ and get the pressure variation

* for cutoff frequency = 1596Hz (first order node) we plot variation of pressure, velocity along the length

* matlab code to find variation of pressure and velocity along the length of the tube for time interval ($t=50sec$)

---

 clear all;
clc

x1 = 0.22; % distance of microphones from source
x2 = 0.38;
x3 = 0.51;
x4 = 0.62;
x5 = 0.7;
speed = 346.02753;
rho = 1.2041;
P = zeros(1,1000);
a = zeros(1,1000);
U = zeros(1,1000);
b = zeros(1,1000);
x = zeros(1,1000);
c = zeros(1,1000);
X = zeros(1,1000);
d = zeros(1,1000);
f = 1600; % cutoff frequency of first mode
C = p+(f); % forward propagation constant
D = p_(f); % backward propagation constant
k = 2*pi*f/speed; % wave number
w=2*pi*f; 
nc=1;
t = 0;
while t<=50
for j = 1:1:1000
    x(1,j) =  1*j*0.001;
    a(1,j) =  cos(k*x(1,j));
    b(1,j) =  sin(k*x(1,j));
    c(1,j) =  cos(w*t);
    d(1,j) =  sin(w*t);
    P(1,j) = (2*C*a(1,j)+ 2*D*b(1,j))*c(1,j); % pressure variation
    U(1,j) = ((2*D*a(1,j)- 2*C*b(1,j))*c(1,j))/(rho*speed); % velocity variation

end

figure(nc) % plot of pressure variation along the length fot time interval t = 0to50 sec
plot(x,P,'r-','LineWidth' ,2)
xlabel ('length(l)')
ylabel ('pressure variation')
title ("pressure vs length of the tube at t="+num2str(t)+"seconds");

figure(nc+1) % plot of pressure variation along the length fot time interval t = 0to50 sec
plot(x,U,'m-','LineWidth' ,2)
xlabel ('length(l)')
ylabel ('velocity variation')
title ("velocity vs length of the tube at t="+num2str(t)+"seconds")
t=t+10;
nc=nc+2;
end

---
* the reflection coefficient is found using the formula
[^1]
[^1]:Seung-Ho Jang and Jeong-Guon Ih , “On the multiple microphone method for measuring in-duct acoustic properties in the presence of mean flow”, The Journal of the Acoustical Society of America 103, 1520-1526 (1998) https://doi.org/10.1121/1.421289

* simplifying the abovve equation in our terms we get i.e, $\Gamma_+= -k$ , $\Gamma_-= k$ , $\Gamma_+ = \Gamma_+^*$ and $\Gamma_- = \Gamma_-^*$ i.e, wave number is real and the complex wave number is neglected (attenuation), we get ref coeff
($|R|$)

$$R(f)= \frac {(\sum_j exp(-2ikx_j)) (\sum_j p_j exp(ikx_j)) } {(\sum_j exp(2ikx_j)) (\sum_j p_j exp(-ikx_j))} $$

* impedance at $x=0$ (mentioned in the above figure) can be obtained as
$$\frac{z(f)}{\rho c} = \frac{1+R(f)}{1-R(f)}$$

* matlab code to find out reflection, absorption coefficient and impedance along the frequency range
---

clear all;
clc

x1 = 0.22;
x2 = 0.38;
x3 = 0.51;
x4 = 0.62;
x5 = 0.7;
speed = 346.02753;
rho = 1.2041;
f=zeros(1,2251);
k=zeros(1,2251);
a=zeros(1,2251);
b=zeros(1,2251);
c=zeros(1,2251);
d=zeros(1,2251);
e=zeros(1,2251);
A=zeros(1,2251);
B=zeros(1,2251);
C=zeros(1,2251);
D=zeros(1,2251);
E=zeros(1,2251);
F=zeros(1,2251);
Z=zeros(1,2251);
Z1=zeros(1,2251);
R=zeros(1,2251);
z=zeros(1,2251);
z1=zeros(1,2251);

P1 % rms value of pressure in frequency domain of mic1
P2 % rms value of pressure in frequency domain of mic2
P3 % rms value of pressure in frequency domain of mic3
P4 % rms value of pressure in frequency domain of mic4
P5 % rms value of pressure in frequency domain of mic5 (absolutely calibrated microphone)

for j= 1:1:2251
        f(1,j) = 249 + 1*j;         %frequency 250 Hz to 2500 Hz
        k(1,j) = 2*pi*f(1,j)/speed; %calculating wave number
        a(1,j) = 1i*k(1,j)*x1;
        b(1,j) = 1i*k(1,j)*x2;
        c(1,j) = 1i*k(1,j)*x3;
        d(1,j) = 1i*k(1,j)*x4;
        e(1,j) = 1i*k(1,j)*x5;
        A(1,j) = exp(-2*a(1,j))+exp(-2*b(1,j))+exp(-2*c(1,j))+exp(-2*d(1,j))+exp(-2*e(1,j));
        B(1,j) = P1*exp(-a(1,j))+P2*exp(-b(1,j))+P3*exp(-c(1,j))+P4*exp(-d(1,j))+P5*exp(-e(1,j));
        C(1,j) = exp(2*a(1,j))+exp(2*b(1,j))+exp(2*c(1,j))+exp(2*d(1,j))+exp(2*e(1,j));
        D(1,j) = P1*exp(a(1,j))+P2*exp(b(1,j))+P3*exp(c(1,j))+P4*exp(d(1,j))+P5*exp(e(1,j));
        E(1,j) = A(1,j)*B(1,j);
        F(1,j) = C(1,j)*D(1,j);
        R(1,j) = E(1,j)/F(1,j);
        z(1,j) = abs(R(1,j)); % calcluating reflection coefficient 
        E(1,j) = 1+R(1,j);
        F(1,j) = 1-R(1,j);
        Z(1,j) = E(1,j)/F(1,j);
        Z1(1,j) = Z(1,j)*rho*speed; % calculating impedance
        z1(1,j) = 1-(z(1,j)^2); % calculating absorption coefficient
end 
        
        plot(f,z,'r-','LineWidth' ,2) % plot of reflection coefficient along the freq range
        xlabel ('frequency (f)')
        ylabel ('Reflection coefficient')
        title ('Reflection coefficient vs frequency')
        
        plot(f,real(Z1),'r',f,imag(Z1),'b-','LineWidth' ,2)  % plot of real and imaginary part of impedance along the freq range
        xlabel('{\itFrequency} [Hz]')
        ylabel('impedance')
        title('impedance vs frequency')
        legend('real part of impedance','imaginary part of impedance')
        xlim ([250 2500])
        
        plot(f,z2,'k-','LineWidth',2) % plot of absorption coefficient along the freq range
        xlabel ('frequency (f)')
        ylabel ('absorption coefficient')
        title ('absorption coefficient vs frequency')
        xlim ([250 2500])

---
* the values of pressure, velocity along the length of the tube - reflection, absorption coefficent and impedance along the freq range are calculated for seven different cases i.e, with and without foam at speaker end (source) and varying the amplitudes of the input signal