# MATLAB CODE TO FIND THE CALIBRATION COEFFICIENTS OF PHASE AND PRESSURE OF MICROPHONES
---

clear
close all
Fs = 8192; % sampling frequency

data_struct = lvm_import('test_4258.lvm'); % importing lvm data
datasp = data_struct.Segment1.data(:,1);
dataAC = data_struct.Segment1.data(:,6);
datamic1 = data_struct.Segment1.data(:,5);
datamic2 = data_struct.Segment1.data(:,3);
datamic3 = data_struct.Segment1.data(:,4);
datamic4 = data_struct.Segment1.data(:,2);
% number of samples = 819200
time = 0:100/819199:100; % time duration 

v_sp = datasp'; % voltage
v_AC = dataAC';
v_datamic1 = datamic1' ;
v_datamic2 = datamic2';
v_datamic3 = datamic3';
v_datamic4 = datamic4';

p_AC = 1000*v_AC; % pressure in absolutely calibrated microphone
p_datamic1 = 1000*v_datamic1; %pressure of mic 1
p_datamic2 = 1000*v_datamic2; %pressure of mic 2
p_datamic3 = 1000*v_datamic3; %pressure of mic 3
p_datamic4 = 1000*v_datamic4; %pressure of mic 4

% fft of absolutely calibrated microphone

j_AC = length(p_AC);
k_AC = 2^nextpow2(j_AC);   % frequency spectrum for AC microphone
freq_AC = Fs*(0:k_AC/2-1)/k_AC; % frequency range
ff_AC = fft(p_AC,k_AC);
z_AC= ff_AC(1:k_AC/2); % obtaining the value of pressure using fft

a_AC = abs(z_AC); % absolute value of pressure amplitude
a1_AC = a_AC-mean(a_AC); 


% pwelch of AC

Nx_AC = length(p_AC);
% length of signal 
nsc_AC = floor(Nx_AC/7);
% dividing the signal into eight sections of equal length
nov_AC = floor(nsc_AC);
% overlap samples
nfft_AC = max(256,2^nextpow2(nsc_AC));
window_AC = hanning(nfft_AC); 
% we use hanning window

kk_AC = 2^nextpow2(Nx_AC/7);
ff_AC = fft(p_AC,kk_AC);
phi_AC = angle(ff_AC)
% calculating the phase angle of absolutely calibrated microphone

freq_0 = Fs*(0:kk_AC/2)/kk_AC; 
% frequency range

pxx = pwelch(p_AC); 
t1_AC = pwelch(p_AC,window_AC,nov_AC,nfft_AC);
% using pwelch command 

kk1_AC = 10*log10(t1_AC); 
% power spectral density of absolute calibrated microphone

% fft of mic1
j_datamic1 = length(p_datamic1);
k_datamic1 = 2^nextpow2(j_datamic1);   
% frequency spectrum for microphone 1
freq_datamic1 = Fs*(0:k_datamic1/2-1)/k_datamic1;
% frequency range
ff_datamic1 = fft(p_datamic1,k_datamic1);
z_datamic1 = ff_datamic1(1:k_datamic1/2); 
% obtaining the value of pressure using fft

a_datamic1 = abs(z_datamic1);
% absolute value of pressure amplitude of mic1
a1_datamic1 = a_datamic1-mean(a_datamic1);


% pwelch of mic1
Nx_mic1 = length(p_datamic1); 
% length of signal 
nsc_mic1 = floor(Nx_mic1/7);
% dividing the signal into eight sections of equal length
nov_mic1 = floor(nsc_mic1); % overlap samples
nfft_mic1 = max(256,2^nextpow2(nsc_mic1));
window_mic1 = hanning(nfft_mic1); % we use hanning window
kk_mic1 = 2^nextpow2(Nx_mic1/7);
ff_mic1 = fft(p_mic1,kk_mic1);
phi_mic1 = = angle(ff_mic1);
% calculating the phase angle of mic1

freq_1 = Fs*(0:kk_mic1/2)/kk_mic1; 
% frequency range

pxx_mic1 = pwelch(p_datamic1); 
t1_mic1 = pwelch(p_datamic1,window_mic1,nov_mic1,nfft_mic1); % using pwelch method
kk1_mic1 = 10*log10(t1_mic1);
% power spectral density of mic1


% fft of mic2
j_datamic2 = length(p_datamic2);
k_datamic2 = 2^nextpow2(j_datamic2);   
% frequency spectrum for microphone 2
freq_datamic2 = Fs*(0:k_datamic2/2-1)/k_datamic2;
% frequency range
ff_datamic2 = fft(p_datamic2,k_datamic2);
z_datamic2 = ff_datamic2(1:k_datamic2/2); 
% obtaining the value of pressure using fft

a_datamic2 = abs(z_datamic2);
% absolute value of pressure amplitude of mic1
a1_datamic2 = a_datamic2-mean(a_datamic2);


% pwelch of mic2
Nx_mic2 = length(p_datamic2); 
% length of signal 
nsc_mic2 = floor(Nx_mic2/7);
% dividing the signal into eight sections of equal length
nov_mic2 = floor(nsc_mic2); % overlap samples
nfft_mic2 = max(256,2^nextpow2(nsc_mic2));
window_mic2 = hanning(nfft_mic2); % we use hanning window
kk_mic2 = 2^nextpow2(Nx_mic2/7);
ff_mic2 = fft(p_mic2,kk_mic2);
phi_mic2 = angle(ff_mic2);
% calculating the phase angle of mic2

freq_2 = Fs*(0:kk_mic2/2)/kk_mic2; % frequency range

pxx_mic2 = pwelch(p_datamic2); 
t1_mic2 = pwelch(p_datamic2,window_mic2,nov_mic2,nfft_mic2); % using pwelch method
kk1_mic1 = 10*log10(t1_mic2);
% power spectral density of mic2

% fft of mic3
j_datamic3 = length(p_datamic3);
k_datamic3 = 2^nextpow2(j_datamic3);   
% frequency spectrum for microphone 3
freq_datamic3 = Fs*(0:k_datamic3/2-1)/k_datamic3;
% frequency range
ff_datamic3 = fft(p_datamic3,k_datamic3);
z_datamic3 = ff_datamic3(1:k_datamic3/2); 
% obtaining the value of pressure using fft

a_datamic3 = abs(z_datamic3);
% absolute value of pressure amplitude of mic1
a1_datamic3 = a_datamic3-mean(a_datamic3);


% pwelch of mic3
Nx_mic3 = length(p_datamic3); 
% length of signal 
nsc_mic3 = floor(Nx_mic3/7);
% dividing the signal into eight sections of equal length
nov_mic3 = floor(nsc_mic3);
% overlap samples
nfft_mic3 = max(256,2^nextpow2(nsc_mic3));
window_mic3 = hanning(nfft_mic3);
% we use hanning window
kk_mic3 = 2^nextpow2(Nx_mic3/7);
ff_mic3 = fft(p_mic3,kk_mic3);
phi_mic3 = angle(ff_mic3);
% calculating the phase angle of mic3

freq_3 = Fs*(0:kk_mic3/2)/kk_mic3;
% frequency range

pxx_mic3 = pwelch(p_datamic3); 
t1_mic3 = pwelch(p_datamic3,window_mic3,nov_mic3,nfft_mic3);
% using pwelch method
kk1_mic3 = 10*log10(t1_mic3);
% power spectral density of mic3

% fft of mic4
j_datamic4 = length(p_datamic4);
k_datamic4 = 2^nextpow2(j_datamic4);  
% frequency spectrum for microphone 4
freq_datamic4 = Fs*(0:k_datamic4/2-1)/k_datamic4;
% frequency range
ff_datamic4 = fft(p_datamic4,k_datamic4);
z_datamic4 = ff_datamic4(1:k_datamic4/2);
% obtaining the value of pressure using fft

a_datamic4 = abs(z_datamic4);
% absolute value of pressure amplitude of mic4
a1_datamic4 = a_datamic4-mean(a_datamic4);


% pwelch of mic4
Nx_mic4 = length(p_datamic4); 
% length of signal 
nsc_mic4 = floor(Nx_mic4/7);
% dividing the signal sections into equal length
nov_mic4 = floor(nsc_mic4); 
% overlap samples
nfft_mic4 = max(256,2^nextpow2(nsc_mic4));
window_mic4 = hanning(nfft_mic4);
% we use hanning window
kk_mic4 = 2^nextpow2(Nx_mic4/7);
ff_mic4 = fft(p_mic4,kk_mic4);
phi_mic4 = angle(ff_mic4);
% calculating the phase angle of mic4

freq_4 = Fs*(0:kk_mic4/2)/kk_mic4;
% frequency range

pxx_mic4 = pwelch(p_datamic4);
% using pwelch method
kk1_mic3 = 10*log10(t1_mic3);
t1_mic4 = pwelch(p_datamic4,window_mic4,nov_mic4,nfft_mic4);
% power spectral density of mic4

kk1_mic4 = 10*log10(t1_mic4);


% calibration coefficient for pressure of microphones

k1_datamic1 = t1_mic1./t_AC;
k1_datamic2 = t1_mic2./t_AC;
k1_datamic3 = t1_mic3./t_AC;
k1_datamic4 = t1_mic4./t_AC;;
k1_AC = t1_AC./t_AC;;

% calibration coefficients for phase of microphones

phi1_AC = phi_AC-phi_AC;
phi1_mic1 = phi_mic1-phi_AC;
phi1_mic2 = phi_mic2-phi_AC;
phi1_mic3 = phi_mic3-phi_AC;
phi1_mic4 = phi_mic4-phi_AC;

              


% plot for calibration coefficient of pressure

k2_datamic1 = [k1_datamic1(4001:32:40001,1)];
freq1_datamic1 = [freq_1(1,4001:32:40001)];
plot(freq1_datamic1,k2_datamic1,'k-','LineWidth',0.5) % graph 
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
ylim([0 2])
% ylim([0 0.1])

hold on

k2_datamic2 = [k1_datamic2(4001:32:40001,1)];
freq1_datamic2 = [freq_2(1,4001:32:40001)];

plot(freq1_datamic2,k2_datamic2,'g-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
% ylim([-2 2])

hold on

k2_datamic3 = [k1_datamic3(4001:32:40001,1)];
freq1_datamic3 = [freq_3(1,4001:32:40001)];

plot(freq1_datamic3,k2_datamic3,'r-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
% ylim([-2 2])

hold on

k2_datamic4 = [k1_datamic4(4001:32:40001,1)];
freq1_datamic4 = [freq_4(1,4001:32:40001)];

plot(freq1_datamic4,k2_datamic4,'m-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
% ylim([0 5])

hold on

k2_AC = [k1_AC(4001:32:40001,1)];
freq1_AC = [freq_0(1,4001:32:40001)];

plot(freq1_AC,k2_AC,'b-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
% ylim([-2 2])

k1 = rms(k2_datamic1);
k2 = rms(k2_datamic2);
k3 = rms(k2_datamic3);
k4 = rms(k2_datamic4);
k5 = rms(k2_AC);

% plot for calibration coefficient of phase

figure(2)

k2_datamic1 = [phi1_mic1(1,1:16:40001)];
freq1_datamic1 = [freq_1(1,1:16:40001)];
plot(freq_1,phi1_mic1,'k-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)')
ylabel('\phi_{ref}-\phi_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
% ylim([-2 2])
ylim([-0.2 0.2])

hold on

k2_datamic2 = [phi1_mic2(1,1:16:40001)];
freq1_datamic2 = [freq_2(1,1:16:40001)];

plot(freq1_datamic2,k2_datamic2,'g-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
% ylim([-2 2])

hold on

k2_datamic3 = [phi1_mic3(1,1:16:40001)];
freq1_datamic3 = [freq_3(1,1:16:40001)];

plot(freq1_datamic3,k2_datamic3,'r-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
% ylim([-2 2])

hold on

k2_datamic4 = [phi1_mic4(1,1:16:40001)];
freq1_datamic4 = [freq_4(1,1:16:40001)];

plot(freq1_datamic4,k2_datamic4,'m-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])
% ylim([0 5])



hold on

k2_AC = [phi1_AC(1,1:16:40001)];
freq1_AC = [freq_0(1,1:16:40001)];

plot(freq1_AC,k2_AC,'b-','LineWidth',0.5) % graph
xlabel('Frequency (Hz)');
ylabel('p_{ref}/p_{j}');
title('FFT presure frequency domain mic 1');
xlim([250 2500])

---