function [specData, featureRate, F, T] = audio_to_spectrogram(audio, parameter)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name: audio_to_spectrogram
% Version: 1
% Date: 20090304
% Programmer: Peter Grosche
%
% Computes various spectrograms
%
%
%  specData    K x N  matrix
%              N = columns (time)
%              K = rows (frequency)
%  featureRate     features per second
%
%  parameter
%       vis_spect   :
%       win_len     : in samples, may be an array of length(parameter.F), than each of the
%                     coefficients is computed using the corresponding
%                     win_len, overlap is adjusted to match hop_size,
%                     resulting in the same featureRate for all coefficients
%       hop_size    : in samples
%       spectMethod : 0;1;2  : complex, magnitude, power
%       F           : Frequencies in Hz, only certain coefficients are computed. If
%                     omitted, standard DFT spacing of coefficients is used.
%       nFFT        : number of coefficients
%       blockSize   : signal is devided into blocks of blockSize length in
%                     sec to allow longer files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

if nargin<2
    parameter=[];
end

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

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

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

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

if isfield(parameter,'spectMethod')==0
    parameter.spectMethod = 1;
end
if isfield(parameter,'fs')==0
    parameter.fs = 1;
end
if isfield(parameter,'F')==0
    parameter.F = 0;
end
if isfield(parameter,'blockSize')==0
    parameter.blockSize = 180;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

win_len = parameter.win_len;
featureRate = parameter.fs./parameter.hop_size;
numBlocks = ceil( (length(audio)+win_len(1)) /(parameter.blockSize.*22050) * (win_len(1)/parameter.hop_size)/2 ); % one block is 3 min of audio at 22050 with overlap = 0.5

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute spectrogam
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if numBlocks == 1 || length(parameter.win_len) > 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % single block
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if (isempty(parameter.F) || ~any(parameter.F))
        audio = [zeros(round(win_len./2),1); audio;zeros(round(win_len./2),1)];
        overlap = win_len - parameter.hop_size;
        [specData,F,T] = spectrogram(audio,hann(win_len),overlap, parameter.nFFT, parameter.fs );
%         specData = specData./sqrt(win_len) / sum(hann(win_len)) * win_len;
        T = T - T(1);     % first window is centered around zero
    elseif length(parameter.win_len) == 1
        audio = [zeros(round(win_len./2),1); audio;zeros(round(win_len./2),1)];
        overlap = win_len - parameter.hop_size;
        [specData,F,T] = compute_fourierCoefficients(audio, hann(win_len),overlap, parameter.F(:), parameter.fs);
%         specData = specData./sqrt(win_len) / sum(hann(win_len)) * win_len;
        T = T - T(1);     % first window is centered around zero
    elseif length(win_len) == length(parameter.F)
        win_num = floor(length(audio)./parameter.hop_size)+1;
        specData = zeros(length(parameter.F),win_num);
        F = zeros(length(parameter.F),1);
        for idx = 1:length(parameter.F)
            f_audio = [zeros(round(win_len(idx)./2),1); audio;zeros(round(win_len(idx)./2),1)];
            overlap = win_len(idx) - parameter.hop_size;
            [temp,F(idx),T] = compute_fourierCoefficients(f_audio, hann(win_len(idx)),overlap, parameter.F(idx), parameter.fs);
            specData(idx,1:size(temp,2)) = temp./sqrt(win_len(idx))/sum(hann(win_len(idx))) * win_len(idx);
        end
        T = T - T(1);     % first window is centered around zero
        clear f_audio
    else
        error('Wrong parameter!');
    end


    switch parameter.spectMethod
        case 0, % complex

        case 1, % magnitude
            specData = abs(specData);
        case 2, % power
            specData = abs(specData).^2;
    end




else
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % multiple blocks
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    audio = [zeros(round(win_len./2),1); audio;zeros(round(win_len./2),1)];

    overlap = win_len - parameter.hop_size;
    win_num = fix((length(audio)-overlap)/(win_len-parameter.hop_size));
    winPerBlock = ceil(win_num/numBlocks);
    if (isempty(parameter.F) || ~any(parameter.F))
        specData = zeros(round(parameter.nFFT/2)+1,win_num);
    else
        specData = zeros(length(parameter.F),win_num);
    end

    for block = 1:numBlocks

        start = max(1,((block-1)*winPerBlock.*parameter.hop_size)+1);
        stop  = min(length(audio),block*winPerBlock.*parameter.hop_size + round(win_len./2));

        startWin = max(1,(block-1)*winPerBlock+1);
        stopWin = min(win_num,block*winPerBlock);

        if (isempty(parameter.F) || ~any(parameter.F))
            switch parameter.spectMethod
                case 0, % complex
                    [specData(:,startWin:stopWin)] = spectrogram(audio(start:stop),hann(win_len),overlap, parameter.nFFT, parameter.fs );
                case 1, % magnitude
                    [specData(:,startWin:stopWin)] = abs(spectrogram(audio(start:stop),hann(win_len),overlap, parameter.nFFT, parameter.fs ));
                case 2, % power
                    [specData(:,startWin:stopWin)] = abs(spectrogram(audio(start:stop),hann(win_len),overlap, parameter.nFFT, parameter.fs )).^2;
            end
            F = (0:size(specData,1)-1)'./(size(specData,1)-1) .*parameter.fs/2;
        elseif length(parameter.win_len) == 1
            switch parameter.spectMethod
                case 0, % complex
                    specData(:,startWin:stopWin) = compute_fourierCoefficients(audio(start:stop), hann(win_len),overlap, parameter.F(:), parameter.fs);
                case 1, % magnitude
                    specData(:,startWin:stopWin) = abs(compute_fourierCoefficients(audio(start:stop), hann(win_len),overlap, parameter.F(:), parameter.fs));
                case 2, % power
                    specData(:,startWin:stopWin) = abs(compute_fourierCoefficients(audio(start:stop), hann(win_len),overlap, parameter.F(:), parameter.fs)).^2;
            end
            F = parameter.F(:);
        else
            error('Block processing and different window lengths doesnt work')
        end
        T = (0:size(specData,2)-1)./parameter.fs.*parameter.hop_size;
    end

end






%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Visualization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if parameter.vis_spect
    figure
    imagesc(T,F,abs(specData)); title('Spectrogram');
    colorbar;
    set(gca,'YDir','normal')
end
