% This m-file was written by Justin Zito and Steven Wright to extract a
% drum track from a single channel recording with two instruments.

close all
clear all
clc

[aa8,fs,nbits]=wavread('funky2.wav');
taa8=(1:length(aa8))/fs;

figure
plot(taa8,aa8);
xlabel('Time (sec)');
ylabel('Amplitude');
title('Original Recording')
% figure
% subplot(2,1,1);
% plot(abs(fft(aa8,fs)));
% title('Magnitude Spectrum of Phoneme /aa8/ by J. Zito');
% xlabel('Frequency (Hz)');
% ylabel('Amplitude');
%
% subplot(2,1,2);
% plot(20*log(abs(fft(aa8,fs))));
% title('Magnitude Spectrum of Phoneme /aa8/ by J. Zito');
% xlabel('Frequency (Hz)');
% ylabel('Amplitude');
%aa8=aa8';
N=3;
scale=44100;
Wn1=[140/scale 190/scale];
Wn2=[290/scale 360/scale];
Wn3=[580/scale 710/scale];  
%Wn4=[1300/scale 1500/scale];
%Wn5=[2700/scale 2900/scale];
Wn6=[5500/scale 5700/scale];
Wn7=[10500/scale 12000/scale];
Wn4=[14000/scale 22049/scale]; 
Wn5=[1/scale 11000/scale]; %
%Wn5=[1700/scale 1800/scale];
Wn9=[50/scale 200/scale];

[B1,A1] = BUTTER(N,Wn1);%,'stop');
[B2,A2] = BUTTER(N,Wn2);%,'stop');
[B3,A3] = BUTTER(N,Wn3);%,'stop');
[B4,A4] = BUTTER(N,Wn4);%,'stop');
[B5,A5] = BUTTER(N,Wn5);%,'stop');

Y1 = filter(B1,A1,aa8);
Y2 = filter(B2,A2,aa8);
Y3 = filter(B3,A3,aa8);
Y4 = filter(B4,A4,aa8);
Y5 = filter(B5,A5,aa8);

Fs = 44100; Fo = 147;  Q = 1.4; BW = (Fo/(Fs/2))/Q;  %176,155,147,160,137.5    Q=1.4 interesting
[b6,a6] = iircomb(round(Fs/Fo),BW,6,'notch');
%fvtool(b6,a6);
Y6=filter(b6,a6,aa8);

Fo=150; BW = (Fo/(Fs/2))/Q;
[b7,a7] = iircomb(round(Fs/Fo),BW,6,'notch'); %'peak' 'notch'
%fvtool(b7,a7);
Y7=filter(b7,a7,Y6);

Fo=155; BW = (Fo/(Fs/2))/Q;
[b8,a8] = iircomb(round(Fs/Fo),BW,6,'notch');
%fvtool(b8,a8);
Y8=filter(b8,a8,Y7);

Fo=160; BW = (Fo/(Fs/2))/Q;
[b9,a9] = iircomb(round(Fs/Fo),BW,6,'notch');
%fvtool(b8,a8);
Y9=filter(b9,a9,Y8);

Fo=176; BW = (Fo/(Fs/2))/Q;
[b10,a10] = iircomb(round(Fs/Fo),BW,6,'notch');
%fvtool(b8,a8);
Y10=filter(b10,a10,Y9);


Y = Y4+Y10;

t=(1:length(Y))/fs;
figure
subplot(2,1,1);
plot(t,Y);
subplot(2,1,2);
plot(abs(fft(Y,fs)));

%Yn=randn(size(Y,1),size(Y,2))*.0001;
%Y=Y+Yn;
% t=(1:length(Y1+Y2+Y3))/fs;
% figure
% subplot(2,1,1);
% plot(t,y);
% subplot(2,1,2);
% plot(abs(fft(y,fs)));
%Y=Y*10;

soundsc(Y,fs);
%wavwrite(Y,fs,nbits,'final_drum_track.wav');
taa8=(1:length(Y))/fs;

figure
plot(taa8,Y);
xlabel('Time (sec)');
ylabel('Amplitude');
title('Final Drum Track')