function [f_magstft,f_phase,sideinfo] = audio_to_STFT(f_audio,parameter,sideinfo)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name: audio_to_STFT
% Version: 2
% Date: 2009
% Programmer: Sebastian Ewert, Sebastian Kreuzer
%
% - Computing and saving STFT (short-time fourier transform) features
% - can also do spectral whitening
%
% Input: f_audio
%
%        parameter.win_len = 4096
%        parameter.win_ov = 1891
%        parameter.window = @hamming
%        parameter.useMatlabSpectrogram = 1
%        parameter.applySpectralWhitening = 0
%        parameter.SWnumFilter = 30
%        parameter.SWamountOfSpecWhite = 0.33  ( see Klapuri Paper and
%                                                below: nu)
%
%        parameter.vis = 0;
%        parameter.visSpectralWhitening = 0;
%        parameter.save = 0;
%        parameter.save_dir;
%        parameter.save_filename;
%
%        parameter.fs
%
% Output: f_pitch
%
%        sideinfo.STFT.win_len
%        sideinfo.STFT.win_ov
%        sideinfo.STFT.featureRate
%        sideinfo.STFT.freqInfo
%
%        Information f_magstft and sideinfo is stored in mat-file:
%        save(strcat(parameter.dir_write,filename),'f_magstft','sideinfo');
%
% Notes:
%   regarding useMatlabSpectrogram: Matlabs spectrogram is very flexible,
%   but not available everywhere. This parameter can be used to switch to our old code
%   based on standard fft, which should be on every matlab system. This
%   should produce the same values as spectrogram for common parameters,
%   except that spectrogram usually returns one feature vector less
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
    sideinfo=[];
end

if nargin<2
    parameter=[];
end

if nargin<1
    error('Please specify input data f_audio');
end

if isfield(parameter,'save')==0
    parameter.save = 0;
end

if isfield(parameter,'vis')==0
    parameter.vis = 0;
end

if isfield(parameter,'visSpectralWhitening')==0
    parameter.visSpectralWhitening = 0;
end

if isfield(parameter,'win_ov')==0
    parameter.win_ov = 1891;
end

if isfield(parameter,'win_len')==0
    parameter.win_len = 4096;
end

if isfield(parameter,'window')==0
    parameter.window = @hamming;
end

if isfield(parameter,'useMatlabSpectrogram')==0
    parameter.useMatlabSpectrogram = 1;
end

if isfield(parameter,'applySpectralWhitening')==0
    parameter.applySpectralWhitening = 0;
end

if isfield(parameter,'SWnumFilter')==0
    parameter.SWnumFilter = 30;
end

if isfield(parameter,'SWamountOfSpecWhite')==0
    parameter.SWamountOfSpecWhite = 0.33;
end

if isfield(parameter,'fs')==0
    parameter.fs = 22050;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some pre calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


win_len = parameter.win_len;
win_ov = parameter.win_ov;
featureRate = parameter.fs/(win_len-win_ov);
wav_size = length(f_audio);
win = window(parameter.window,win_len);
num_coeffs = floor(win_len/2)+1;
first_win = floor(win_len/2);
step_size = win_len-win_ov;
num_frames = ceil((wav_size-win_len+first_win)/step_size);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spectrogram calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if parameter.useMatlabSpectrogram
    % first window should center at 0 seconds
    first_win = floor(win_len/2);
    f_audio = [zeros(first_win,1); f_audio];

    % The frequency information in f seems has been checked. It gives the exact
    % frequency in Hz of the analysis function for the corresponding bin 
    [Xs,f,t] = spectrogram(f_audio,parameter.window(win_len),win_ov,[],parameter.fs);
    f_magstft = abs(Xs); % magnitude spectrum
    f_phase = angle(Xs);
    clear Xs

    % we added half a window. We need to adjust the time info
    t = t - first_win/parameter.fs;
else

    % allocate memory
    f_magstft = zeros(num_coeffs,num_frames);

    % first window should center at 0 seconds
    frame = (1:win_len)-first_win;
    x = [zeros(first_win,1); f_audio(1:(win_len-first_win))].*win;

    i = 1;
    while frame(end) < wav_size
        %[Xs,f] = computeDFT(x,win_len,parameter.fs);
        Xs = fft(x,win_len);
        %f = [0:size(Xs,1)-1]'/(size(S,1)-1) * (22050 / 2);
        f_magstft(:,i) = abs(Xs(1:num_coeffs)); % magnitude spectrum
        f_phase(:,i) = angle(Xs(1:num_coeffs));
        clear Xs;
        frame = frame+step_size;
        if frame(end) > wav_size
            frame = frame(1):wav_size;
            win = win(1:length(frame));
        end
        x = f_audio(frame).*win;
        i = i+1;
    end

    % calculate last window
    ft = abs(fft(x,win_len));
    f_magstft(:,i) = ft(1:num_coeffs);  % check here if num_coeffs makes sense

    t = [0:size(f_magstft,2)-1] * step_size/parameter.fs;
    f = [0:size(f_magstft,1)-1]'/(size(f_magstft,1)-1) * (22050 / 2);
end

if parameter.applySpectralWhitening
    % see 'Multiple Fundamental Frequency Estimation by Summing Harmonic
    % Amplitudes' (Klapuri)

    nFilters = parameter.SWnumFilter;

    % c_b and c_b_indices have nFilters+2 entries. Filter center
    % frequencies are at c_b(2:end-1). c_b(1) and c_b(end) are for later
    % usuage
    c_b = 229 * (10.^(([0:nFilters+1]+1)/21.4) - 1);
    c_b_indices = zeros(nFilters+2,1);
    for k=0:nFilters+1
        [ignoreMe c_b_indices(k+1)] = min(abs(f-c_b(k+1)));
    end

    numCoeffs = size(f_magstft,1);
    numFrames = size(f_magstft,2);

    % H_b Filterbank and calculate the sigma value for each filter
    % and frame
    sigmasPerFilter = zeros(nFilters,numFrames);
    H_b = zeros(numCoeffs,nFilters);
    for k=1:nFilters
        startIndex = c_b_indices(k);
        centerIndex = c_b_indices(k+1);
        stopIndex = c_b_indices(k+2);

        H_b(startIndex:centerIndex,k) = (0:centerIndex-startIndex)/ (centerIndex-startIndex);
        H_b(centerIndex:stopIndex,k) = (stopIndex-centerIndex:-1:0)/ (stopIndex-centerIndex);

        sigmasPerFilter(k,:) = sqrt(mean(repmat(H_b(:,k),1,numFrames).* (f_magstft.^2),1)); 
    end

    % Now construct the spectral whitening filter in spectral domain
    nu = parameter.SWamountOfSpecWhite;
    specWhiteFilter = zeros(size(f_magstft));
    for k=0:nFilters
        curIndex = c_b_indices(k+1);
        nextIndex = c_b_indices(k+2);
        if k==0
            curValues = zeros(1,numFrames);
            nextValues = sigmasPerFilter(1,:).^(nu-1);
        elseif k==nFilters
            curValues = sigmasPerFilter(end,:).^(nu-1);
            nextValues = zeros(1,numFrames);
        else
            curValues = sigmasPerFilter(k,:).^(nu-1);
            nextValues = sigmasPerFilter(k+1,:).^(nu-1);
        end
        diffValues = nextValues - curValues;
        specWhiteFilter(curIndex:nextIndex,:) = repmat(curValues,nextIndex-curIndex+1,1) + ...
            ([0:nextIndex-curIndex]/(nextIndex-curIndex))' * diffValues;
    end

    if parameter.visSpectralWhitening && parameter.vis
        f_magstft_org = f_magstft;
    end
    f_magstft = specWhiteFilter .* f_magstft;

end

sideinfo.STFT.win_len = win_len;
sideinfo.STFT.win_ov = win_ov;
sideinfo.STFT.featureRate = featureRate;
sideinfo.STFT.freqInfo = f;
sideinfo.STFT.timeInfo = t;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save stft for each window size separately as f_magstft
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if parameter.save == 1
    if isfield(parameter,'save_dir')==0
        parameter.save_dir = '';
    end
    if isfield(parameter,'save_filename')==0
        parameter.save_filename = '';
    end

    filename = strcat(parameter.save_filename,'_STFT_',num2str(win_len));
    save(strcat(parameter.save_dir,filename),'f_magstft','f_phase','sideinfo');
else

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if parameter.vis
    figure;
    hot2 = hot;
    hot2 = [hot2; hot2(32:64,:); hot2(32:64,:)];
    hot2 = sort(hot2);
    imagesc(t,f,f_magstft);
    set(gca,'YDir','normal');
    colormap(hot2);
    colorbar;

    if parameter.visSpectralWhitening==1 && parameter.applySpectralWhitening==1
        figure;
        imagesc(t,f,f_magstft_org);
        set(gca,'YDir','normal');
        colormap(hot2);
        colorbar;

        figure;
        imagesc(t,f,specWhiteFilter);
        axis xy;
        colorbar

        figure;
        bar(H_b(:,1))
        hold on;
        for k=1:nFilters
            hPlot = plot(H_b());
        end
    end
end

end
