%--------------------------------------------------
% Main code for spectogram generation
% Example usage:
%   [orig synth Fs] = main('phonetician.wav',200,500);
%--------------------------------------------------

function [y z Fs] = main(filename,filterlow,filterhigh)
%--------------------------------------------------
% Load sound file
%--------------------------------------------------
[y Fs nBits] = wavread(filename);
sound(y,Fs);
Y = length(y);

%--------------------------------------------------
% Knobs
%--------------------------------------------------
N_desired = 80; 

L = floor(Y/N_desired);
L = 2^(nextpow2(L)-1);

%--------------------------------------------------
% Cut the sound stream into parts each of length L
%--------------------------------------------------
s = split(y,L);
N = length(s(:,1)); % Actual number of partitions

%--------------------------------------------------
% Prepare spectogram
%--------------------------------------------------
f = Fs/2*linspace(0,1,L/2);

S = zeros(N,L);
for i = 1:N
    S(i,:) = fft(s(i,:),L)/L;
    
    plot(f,log(2*abs(S(i,1:L/2))));
    drawnow;
end

%--------------------------------------------------
% Visualize spectogram
%--------------------------------------------------
t     = linspace(0,Y/Fs,N);
[X Y] = meshgrid(t,f);
% Decide on linear vs. log plot
%surf(X,Y,2*abs(S(:,1:L/2))');
surf(X,Y,log(2*abs(S(:,1:L/2)))');
xlabel('t (in seconds)');
ylabel('f (in Hz)');
title('Spectogram of original sound file');

%--------------------------------------------------
% Apply filter
%--------------------------------------------------
F = zeros(1,L);

% Pass region
F = zeros(1,L);
a1 = floor(filterlow/(Fs/L))
a2 = floor(filterhigh/(Fs/L))
if (a1 ~= a2)
    F(max(a1,1):a2) = 1;
end

for i = L/2+1:L-1
    F(i) = F(L-i);
end

for i = 1:N
    S(i,:) = F.*S(i,:);
end

figure;
surf(X,Y,log(2*abs(S(:,1:L/2)))');
xlabel('t (in seconds)');
ylabel('f (in Hz)');
title('After applying filter');

%--------------------------------------------------
% Resynthesized
%--------------------------------------------------
SS = zeros(N,L);
for i = 1:N
    SS(i,:) = L*ifft(S(i,:),L);
end
z = appendv(SS);

figure
subplot(211);
plot(real(z));
hold on;
plot(imag(z),'r');

subplot(212);
plot(y);
sound(real(z),Fs);
z = real(z);

% Save the result
filename = filename(1:length(filename)-4);
outFile = strcat(filename,'_',int2str(filterlow),'_',int2str(filterhigh));
wavwrite(z,Fs,outFile);