clear all;
close all;

[y, fs] = audioread('mariuana.wav');

% lower the damping factor the smaller the pass band
damp = 0.05;

% min and max centre cutoff frequency of variable bandpass filter
minf=300;
maxf=2000;
% wah frequency, how many Hz per second are cycled through
Fw = 2000;

% change in centre frequency per sample (Hz)
delta = Fw/fs;

no_points = length(y(:,1));
fwd_back = 0;
ifs = 0;
z = ((maxf-minf)/delta);
yh = 1:2;
yb = 1:2;
y1 = 1:2;
for i=1:no_points
    if(fwd_back==0)
       ifs = ifs +1;
       delay(i) = ((maxf-minf)/z)*ifs+minf;
    end
    
    if(ifs == uint32(z))
      fwd_back = 1;
      ifs = ifs -1;
    end
    
     if(fwd_back==1)
       ifs = ifs -1;
       delay(i) = ((maxf-minf)/z)*ifs+minf;
       if(ifs == 0)
           fwd_back = 0;
       end
     end
end

F1 = 2*sin((pi*delay(1))/fs);
Q1 = 2*damp;
yh(1) = y(1);
yb(1) = F1*yh(1);
yl(1) = F1*yb(1);

%for i=2:no_points
 %   yh(i) = y(i) - yl(i-1) - Q1*yb(i-1);
 %   yb(i) = F1*yh(i) + yb(i-1);
 %   yl(i) = F1*yb(i) + yl(i-1);
 %  F1 = 2*sin((pi*delay(i))/fs);
 %   out_wav(i) = yb(i);
%end

for i=2:no_points
   yh(2) = y(i) - y1(1) - Q1*yb(1);
   yb(2) = F1*yh(2) + yb(1);
   y1(2) = F1*yb(2) + y1(1);
   F1 = 2*sin((pi*delay(i))/fs);
%    if (yb(2) > 32678)
%        yb(2) = 32678
%    elsif (yb(2) < -32678)
%        yb(2) = -32678
%    end
       
   out_wav(i) = yb(2);
   yb(1) = yb(2);
   y1(1) = y1(2);
end

p = audioplayer(out_wav, fs);
play(p);
plot(out_wav);