[y fS nbits] = wavread('preparetodie.wav');

%------------------------------------------------------------
% Prepare spectogram
%------------------------------------------------------------
t = 1/Fs;      % Time per sample
Y = length(y);
T = Y*t;       % Total signal length

N = floor(Y/4096);
s = split(y,N);
L = length(s(1,:));
NFFT = 2^nextpow2(L);

f = Fs/2*linspace(0,1,NFFT/2);
S = zeros(N,NFFT);

for i = 1:N
    S(i,:) = fft(s(i,:),NFFT)/L; 
        
    %plot(f,2*S(i,1:NFFT/4));
    %drawnow;
end

% Visualize spectogram
t = linspace(0,T,N);
[X Y] = meshgrid(t,f);

surf(X,Y,2*abs(S(:,1:NFFT/2))');
xlabel('t (in seconds)');
ylabel('f (in Hz)');
title('Prior to filter');

%------------------------------------------------------------
% Prepare filter
%------------------------------------------------------------
F = zeros(1,NFFT);

% Pass region
F = zeros(1,length(S(i,:)));
a1 = floor(1/(Fs/NFFT));
a2 = floor(Fs/(Fs/NFFT));
F(max(a1,1):a2) = 1;
for i = NFFT/2+1:NFFT-1
    F(i) = F(NFFT-i);
end

%figure;
%plot(F);

%------------------------------------------------------------
% Apply filter
%------------------------------------------------------------
for i = 1:N
    S(i,:) = F.*S(i,:);
end

% Visualize processed spectogram
figure
surf(X,Y,2*abs(S(:,1:NFFT/2))');
xlabel('t (in seconds)');
ylabel('f (in Hz)');
title('After applying filter');

%------------------------------------------------------------
% Resynthesize
%------------------------------------------------------------
SS = zeros(N,NFFT);
for i = 1:N
    SS(i,:) = L*ifft(S(i,:),NFFT); 
    %plot(SS(i,:));
    %drawnow;
end
y = append(SS);
figure
subplot(211);
plot(real(y));
hold on;
plot(imag(y),'r');
subplot(212);
plot(handel,'g');
sound(real(y),Fs);