close all;
clear all;

% read wave samples
[xn1 fs] = wavread('micI.wav');    
[yn1 fs] = wavread('piezoI.wav');
[xn2 fs] = wavread('micAh.wav');    
[yn2 fs] = wavread('piezoAh.wav');
[xn3 fs] = wavread('micMi.wav');    
[yn3 fs] = wavread('piezoMi.wav');
[xn4 fs] = wavread('micOh.wav');    
[yn4 fs] = wavread('piezoOh.wav');
[xn5 fs] = wavread('micMe.wav');    
[yn5 fs] = wavread('piezoMe.wav');

N = 8192;         
w = hann(N);
fni = [0:fs/N:fs-fs/N];

% input 'xn' is mic-recorded sound
% output 'yn' is piezo-recorded sound

Xn1 = abs(fft(xn1,N));
Yn1 = abs(fft(yn1,N));
Xn1_N = Xn1/max(Xn1);
Yn1_N = Yn1/max(Yn1);
Hn1 = Yn1_N./Xn1_N;
Xn2 = abs(fft(xn2,N));
Yn2 = abs(fft(yn2,N));
Xn2_N = Xn2/max(Xn2);
Yn2_N = Yn2/max(Yn2);
Hn2 = Yn2_N./Xn2_N;
Xn3 = abs(fft(xn3,N));
Yn3 = abs(fft(yn3,N));
Xn3_N = Xn3/max(Xn3);
Yn3_N = Yn3/max(Yn3);
Hn3 = Yn3_N./Xn3_N;
Xn4 = abs(fft(xn4,N));
Yn4 = abs(fft(yn4,N));
Xn4_N = Xn4/max(Xn4);
Yn4_N = Yn4/max(Yn4);
Hn4 = Yn4_N./Xn4_N;
Xn5 = abs(fft(xn5,N));
Yn5 = abs(fft(yn5,N));
Xn5_N = Xn5/max(Xn5);
Yn5_N = Yn5/max(Yn5);
Hn5 = Yn5_N./Xn5_N;

% for dB plot
Xn1_dB = 20*log10(Xn1);
Yn1_dB = 20*log10(Yn1);
Xn2_dB = 20*log10(Xn2);
Yn2_dB = 20*log10(Yn2);
Xn3_dB = 20*log10(Xn3);
Yn3_dB = 20*log10(Yn3);
Xn4_dB = 20*log10(Xn4);
Yn4_dB = 20*log10(Yn4);
Xn5_dB = 20*log10(Xn5);
Yn5_dB = 20*log10(Yn5);

Xn1_dBN = Xn1_dB - max(Xn1_dB);
Yn1_dBN = Yn1_dB - max(Yn1_dB);
Xn2_dBN = Xn2_dB - max(Xn2_dB);
Yn2_dBN = Yn2_dB - max(Yn2_dB);
Xn3_dBN = Xn3_dB - max(Xn3_dB);
Yn3_dBN = Yn3_dB - max(Yn3_dB);
Xn4_dBN = Xn4_dB - max(Xn4_dB);
Yn4_dBN = Yn4_dB - max(Yn4_dB);
Xn5_dBN = Xn5_dB - max(Xn5_dB);
Yn5_dBN = Yn5_dB - max(Yn5_dB);

W = w(1:N/2)';
H_av = (Hn1+Hn2+Hn3+Hn4+Hn5)/5;

H1 = Hn1(1:N/2);
H2 = Hn2(1:N/2);
H3 = Hn3(1:N/2);
H4 = Hn4(1:N/2);
H5 = Hn5(1:N/2);

% to get filter coeffiecnts of each transfer function
[B1 A1] = invfreqz(H1, W, 10 ,10);
[B2 A2] = invfreqz(H2, W, 10 ,10);
[B3 A3] = invfreqz(H3, W, 10 ,10);
[B4 A4] = invfreqz(H4, W, 10 ,10);
[B5 A5] = invfreqz(H5, W, 10 ,10);

y1_synth = filter(B1,A1,xn1);
y2_synth = filter(B2,A2,xn2);
y3_synth = filter(B3,A3,xn3);
y4_synth = filter(B4,A4,xn4);
y5_synth = filter(B5,A5,xn5);

y1_synthN = y1_synth/max(y1_synth);
y2_synthN = y2_synth/max(y2_synth);
y3_synthN = y3_synth/max(y3_synth);
y4_synthN = y4_synth/max(y4_synth);
y5_synthN = y5_synth/max(y5_synth);

wavwrite(y1_synth,fs,'synthI.wav');
wavwrite(y2_synth,fs,'synthAh.wav');
wavwrite(y3_synth,fs,'synthMi.wav');
wavwrite(y4_synth,fs,'synthOh.wav');
wavwrite(y5_synth,fs,'synthMe.wav');

Y1_synth = abs(fft(y1_synthN,N));
Y2_synth = abs(fft(y2_synthN,N));
Y3_synth = abs(fft(y3_synthN,N));
Y4_synth = abs(fft(y4_synthN,N));
Y5_synth = abs(fft(y5_synthN,N));

Y1_synthN = Y1_synth/max(Y1_synth);
Y2_synthN = Y2_synth/max(Y2_synth);
Y3_synthN = Y3_synth/max(Y3_synth);
Y4_synthN = Y4_synth/max(Y4_synth);
Y5_synthN = Y5_synth/max(Y5_synth);

fni = 0:fs/N:fs-fs/N; 
xbound = 4000;
ybound = 40;
figure(1)
subplot(511)
plot(fni,Xn1_N); xlim([0 xbound]);
title('air conducted sound by microphone(N=8192), singing "I"');
ylabel('Magnitude');
subplot(512)
plot(fni,Xn2_N); xlim([0 xbound]);
title('Ah');
ylabel('Magnitude');
subplot(513)
plot(fni,Xn3_N); xlim([0 xbound]);
title('Mi');
ylabel('Magnitude');
subplot(514)
plot(fni,Xn4_N); xlim([0 xbound]);
title('Oh');
ylabel('Magnitude');
subplot(515)
plot(fni,Xn5_N); xlim([0 xbound]);
title('Me');
xlabel('Frequency (Hz)');
ylabel('Magnitude');

figure(2)
subplot(511)
plot(fni,Yn1_N); xlim([0 xbound]);
title('bone(skull) conducted sound by piezo, singing "I"');
ylabel('Magnitude');
subplot(512)
plot(fni,Yn2_N); xlim([0 xbound]);
title('Ah');
ylabel('Magnitude');
subplot(513)
plot(fni,Yn3_N); xlim([0 xbound]);
title('Mi');
ylabel('Magnitude');
subplot(514)
plot(fni,Yn4_N); xlim([0 xbound]);
title('Oh');
ylabel('Magnitude');
subplot(515)
plot(fni,Yn5_N); xlim([0 xbound]);
title('Me');
xlabel('Frequency (Hz)');
ylabel('Magnitude');

figure(3)
subplot(511)
plot(fni,Hn1); xlim([0 xbound]); %ylim([0 ybound]);
title('transfer funtion from microphone pickup to piezo pickup, singing "I"');
ylabel('Magnitude');
subplot(512)
plot(fni,Hn2); xlim([0 xbound]); %ylim([0 ybound]);
title('Ah');
ylabel('Magnitude');
subplot(513)
plot(fni,Hn3); xlim([0 xbound]); %ylim([0 ybound]); 
title('Mi');
ylabel('Magnitude');
subplot(514)
plot(fni,Hn4); xlim([0 xbound]); %ylim([0 ybound]);
title('Oh');
ylabel('Magnitude');
subplot(515)
plot(fni,Hn5); xlim([0 xbound]); %ylim([0 ybound]);
title('Me');
xlabel('frequency (Hz)');
ylabel('Magnitude');

figure(4)
subplot(511)
plot(fni,Y1_synthN); xlim([0 xbound]);
title('synthesizing bone(skull) conducted sound, I');
ylabel('Magnitude');
subplot(512)
plot(fni,Y2_synthN); xlim([0 xbound]);
title('Ah');
ylabel('Magnitude');
subplot(513)
plot(fni,Y3_synthN); xlim([0 xbound]);
title('Mi');
ylabel('Magnitude');
subplot(514)
plot(fni,Y4_synthN); xlim([0 xbound]);
title('Oh');
ylabel('Magnitude');
subplot(515)
plot(fni,Y5_synthN); xlim([0 xbound]);
title('Me');
xlabel('Frequency (Hz)');
ylabel('Magnitude');

xlow = 0;
xbound = 4000;

figure(5)
subplot(511)
semilogx(fni,Xn1_dBN); xlim([xlow xbound]);
title('air conducted sound by microphone(N=8192), singing "I"');
ylabel('Magnitude (dB)');
subplot(512)
semilogx(fni,Xn2_dBN); xlim([xlow xbound]);
title('Ah');
ylabel('Magnitude (dB)');
subplot(513)
semilogx(fni,Xn3_dBN); xlim([xlow xbound]);
title('Mi');
ylabel('Magnitude (dB)');
subplot(514)
semilogx(fni,Xn4_dBN); xlim([xlow xbound]);
title('Oh');
ylabel('Magnitude (dB)');
subplot(515)
semilogx(fni,Xn5_dBN); xlim([xlow xbound]);
title('Me');
xlabel('Log-scale Frequency (Hz)');
ylabel('Magnitude (dB)');

figure(6)
subplot(511)
semilogx(fni,Yn1_dBN); xlim([xlow xbound]);
title('bone(skull) conducted sound by piezo, singing "I"');
ylabel('Magnitude (dB)');
subplot(512)
semilogx(fni,Yn2_dBN); xlim([xlow xbound]);
title('Ah');
ylabel('Magnitude (dB)');
subplot(513)
semilogx(fni,Yn3_dBN); xlim([xlow xbound]);
title('Mi');
ylabel('Magnitude (dB)');
subplot(514)
semilogx(fni,Yn4_dBN); xlim([xlow xbound]);
title('Oh');
ylabel('Magnitude (dB)');
subplot(515)
semilogx(fni,Yn5_dBN); xlim([xlow xbound]);
title('Me');
xlabel('Log-scale Frequency (Hz)');
ylabel('Magnitude (dB)');