% This m-file was written by Justin Zito and Steven Wright to comb filter
% the trombone 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';

Fs = 44100; Fo = 110;  Q = 1.4; BW = (Fo/(Fs/2))/Q;  
[b1,a1] = iircomb(round(Fs/Fo),BW,6,'notch');

Y1=filter(b1,a1,aa8);

Fo = 123;  Q = 1.4; BW = (Fo/(Fs/2))/Q;  
[b2,a2] = iircomb(round(Fs/Fo),BW,6,'notch');
Y2=filter(b2,a2,Y1);

Fo = 131;  Q = 1.4; BW = (Fo/(Fs/2))/Q;  
[b3,a3] = iircomb(round(Fs/Fo),BW,6,'notch');
Y3=filter(b3,a3,Y2);

Fo = 147;  Q = 1.4; BW = (Fo/(Fs/2))/Q;  
[b4,a4] = iircomb(round(Fs/Fo),BW,6,'notch');
Y4=filter(b4,a4,Y3);

Fo = 165;  Q = 1.4; BW = (Fo/(Fs/2))/Q;  
[b6,a6] = iircomb(round(Fs/Fo),BW,6,'notch');
Y6=filter(b6,a6,Y4);

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

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

Fo=160; BW = (Fo/(Fs/2))/Q;
[b9,a9] = iircomb(round(Fs/Fo),BW,6,'notch');
Y9=filter(b9,a9,aa8);

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

Fo=250; BW = (Fo/(Fs/2))/Q;
[b11,a11] = iircomb(round(Fs/Fo),BW,6,'notch');
%fvtool(b8,a8);
Y11=filter(b11,a11,Y10);

Y=Y8;

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,'cymbals.wav');
taa8=(1:length(Y))/fs;

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