function [f_CENS,sideinfo] = pitchSTMSP_to_CENS(f_pitch,parameter,sideinfo);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name: pitchSTMSP_to_CENS 
% Version: 1
% Date: 12.07.2007
% Programmer: Meinard Mueller, ported by Sebastian Ewert
%
% Description:
% Normalized statistical chroma-based energy distribution feature (CENS).
% The following is computed:
% - energy for each chroma band
% - normalisation of the chroma vectors
% - a CENS feature series for each pair of statistic window and
%   downsampling factor. This includes:
%     - component-wise quantisation of the normalized chroma vectors 
%     - upfirdn filters and downsamples each column of f_stat_help
%     - normalize each vector with its l^2 norm
%
% Input:  
%         f_pitch
%         parameter.midi_min = 21
%         parameter.midi_max = 108
%         parameter.win_res = 10      
%         parameter.vec_energy = [40 20 10 5] / 100;
%         parameter.vec_weight = [ 1 1 1 1]/4;
%         parameter.stat_thresh = 0.001;
%         parameter.stat_window_length_array  = 41;
%         parameter.stat_downsample_array     = 10;
%
%         parameter.vis = 0;
%            parameter.vis_ChromaNorm = 1;
%            parameter.vis_CENS = 0;
%            parameter.vis_ChromaCENS = 0;
%            parameter.vis_segment  (not optional)
%         parameter.save = 0;
%            parameter.save_dir = '';
%            parameter.save_filename = '';
%
% Output: 
%         f_CENS
%         sideinfo.CENS.version
%         sideinfo.CENS.midi_min
%         sideinfo.CENS.midi_max
%         sideinfo.CENS.stat_res
%         sideinfo.CENS.win_res
%         sideinfo.CENS.vec_energy
%         sideinfo.CENS.vec_weight
%         sideinfo.CENS.stat_thresh
%         sideinfo.CENS.stat_window_length_array
%         sideinfo.CENS.stat_downsample_array
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

if nargin<3
   sideinfo=[]; 
end

if nargin<2
   parameter=[]; 
end

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

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

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


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

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

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

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

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

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

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

if isfield(parameter,'vec_energy')==0
   parameter.vec_energy = [40 20 10 5] / 100;
end

if isfield(parameter,'vec_weight')==0
   parameter.vec_weight = [ 1 1 1 1]/4;
end

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update sideinfo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sideinfo.CENS.version = 1;
sideinfo.CENS.midi_min = parameter.midi_min;
sideinfo.CENS.midi_max = parameter.midi_max; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

seg_num = size(f_pitch,1);

% calculate energy for each chroma band
f_chroma_energy = zeros(seg_num,12);
for p=parameter.midi_min:parameter.midi_max
    chroma = mod(p,12)+1;
    f_chroma_energy(:,chroma) = f_chroma_energy(:,chroma)+f_pitch(:,p);
end

% normalize the chroma vectors
f_chroma_energy_distr = zeros(seg_num,12);
for k=1:seg_num
   if sum(f_chroma_energy(k,:)>parameter.stat_thresh)>0
        seg_energy_square = sum(f_chroma_energy(k,:));
        %f_chroma_energy_distr(k,:) = 100*((f_chroma_energy(k,:))/seg_energy_square);        
        f_chroma_energy_distr(k,:) = ((f_chroma_energy(k,:))/seg_energy_square);        
    end
end

% calculate a CENS feature series for each pair of statistic window and downsampling factor  
num_stat_window = length(parameter.stat_window_length_array);
f_CENS_cell = cell(num_stat_window,1);
sideinfo.CENS.stat_res = zeros(1,num_stat_window);
for w=1:num_stat_window
    stat_window_length = parameter.stat_window_length_array(w);
    stat_downsample = parameter.stat_downsample_array(w);
    stat_window = hanning(stat_window_length);
    stat_window = stat_window/sum(stat_window);
    sideinfo.CENS.stat_res(w) = parameter.win_res/stat_downsample;

    % component-wise quantisation of the normalized chroma vectors 
    f_stat_help = zeros(seg_num,12);
    for n=1:length(parameter.vec_energy)
        f_stat_help = f_stat_help + (f_chroma_energy_distr>parameter.vec_energy(n))*parameter.vec_weight(n);
    end
    
    % upfirdn filters and downsamples each column of f_stat_help
    f_chroma_energy_stat = zeros(seg_num,12);
    f_chroma_energy_stat = upfirdn(f_stat_help,stat_window,1,stat_downsample);
    stat_num = ceil(seg_num/stat_downsample);
    cut = floor((stat_window_length-1)/(2*stat_downsample));
    f_chroma_energy_stat = f_chroma_energy_stat((1+cut:stat_num+cut),:);            %adjust group delay
    
    % last step: normalize each vector with its l^2 norm
    unit_vec = ones(1,12);
    unit_vec = unit_vec/norm(unit_vec);
    f_CENS = zeros(stat_num,12);
    for k=1:stat_num
        n = norm(f_chroma_energy_stat(k,:));
        if n==0
            f_CENS(k,:) = unit_vec;
        else   
            f_CENS(k,:) = f_chroma_energy_stat(k,:)/n;
        end
    end    
    f_CENS_cell{w} = f_CENS;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update sideinfo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
sideinfo.CENS.win_res = parameter.win_res;
sideinfo.CENS.vec_energy = parameter.vec_energy;
sideinfo.CENS.vec_weight = parameter.vec_weight;
sideinfo.CENS.stat_thresh = parameter.stat_thresh;
sideinfo.CENS.stat_window_length_array = parameter.stat_window_length_array;
sideinfo.CENS.stat_downsample_array = parameter.stat_downsample_array;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Saving to file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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,'_CENS_',num2str(parameter.stat_window_length_array),'_',num2str(parameter.stat_downsample_array));
    save(strcat(parameter.save_dir,filename),'f_CENS','sideinfo');
end


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

if parameter.vis

    chroma_names                 = ['C ';'C#';'D ';'D#';'E ';'F ';'F#';'G ';'G#';'A ';'A#';'B '];
    figure_position              = [655    70   560   461];
    hot2 = hot(64); hot2 = [hot2; hot2(32:64,:)]; hot2 = sort(hot2);

    % gray2 = log(gray(64)+1.01);
    % gray2 = gray2-gray2(1,1);
    % gray2 = gray2/(gray2(end,1)-gray2(1,1));
    %
    % gray3 = log(log(gray(64)+1.01));
    % gray3 = gray3-gray3(1,1);
    % gray3 = gray3/(gray3(end,1)-gray3(1,1));


    if parameter.vis_ChromaNorm
        figure;
        set(gcf,'renderer','painters');
        set(gcf,'Position',figure_position);
        t1 = max(0,floor(parameter.vis_segment(1)*parameter.win_res));
        t2 = min(size(f_chroma_energy_distr,1)-1,ceil(parameter.vis_segment(2)*parameter.win_res));
        t = [t1:t2]/parameter.win_res;
        h = imagesc(t,[1:12],f_chroma_energy_distr(t1+1:t2+1,:)',[0,1]);
        set(gca,'YTick',[1:12]);
        set(gca,'YTickLabel',chroma_names);
        set(gca,'YDir','normal');
        if t(end)>t(1)
            set(gca,'XLim',[t(1),t(end)]);
        end
        %xlabel('Time in seconds');
        %colormap(gray2);
        colormap(hot2);
        h = colorbar;
        %pos = get(h,'Position'); pos(3)=pos(3)/3; set(h,'Position',pos)
        drawnow;
        if parameter.print == 1
            set(gcf,'PaperPosition',parameter.printPaperPositionImage);
            print('-depsc2',strcat(parameter.printDir,parameter.printFileChromaNorm));
        end
    end

    if parameter.vis_CENS
        for w=1:num_stat_window
            figure;
            set(gcf,'renderer','painters');
            set(gcf,'Position',figure_position);
            t1 = max(0,floor(parameter.vis_segment(1)*sideinfo.CENS.stat_res(w)));
            t2 = min(size(f_CENS_cell{w},1)-1,ceil(parameter.vis_segment(2)*sideinfo.CENS.stat_res(w)));
            t = [t1:t2]/sideinfo.CENS.stat_res(w);
            h = imagesc(t,[1:12],f_CENS_cell{w}(t1+1:t2+1,:)',[0,1]);
            set(gca,'YTick',[1:12]);
            set(gca,'YTickLabel',chroma_names);
            set(gca,'YDir','normal');
            if t(end)>t(1)
                set(gca,'XLim',[t(1),t(end)]);
            end
            %xlabel('Time in seconds');
            %colormap(gray2);
            colormap(hot2);
            h = colorbar;
            %pos = get(h,'Position'); pos(3)=pos(3)/3; set(h,'Position',pos)
            drawnow;

            if parameter.print == 1
                set(gcf,'PaperPosition',parameter.printPaperPositionImage);
                print('-depsc2',strcat(parameter.printDir,parameter.printFileCENS));
            end
        end
    end

    if parameter.vis_ChromaCENS
        for w=1:num_stat_window
            figure;
            set(gcf,'renderer','painters');
            set(gcf,'Position',figure_position);
            t1 = max(0,floor(parameter.vis_segment(1)*sideinfo.CENS.stat_res(w)));
            t2 = min(size(f_CENS_cell{w},1)-1,ceil(parameter.vis_segment(2)*sideinfo.CENS.stat_res(w)));
            t = [t1:t2]/sideinfo.CENS.stat_res(w);
            max_value = max(1.1,1.1*max(max(f_CENS_cell{w})));
            %max_value = 1.05;
            t1_chroma = max(0,floor(parameter.vis_segment(1)*parameter.win_res));
            t2_chroma = min(size(f_chroma_energy_distr,1)-1,ceil(parameter.vis_segment(2)*parameter.win_res));
            t_chroma = [t1_chroma:t2_chroma]/parameter.win_res;
            for p=1:12
                subplot(12,1,12-p+1);
                h = bar(t,f_CENS_cell{w}(t1+1:t2+1,p));
                set(h,'EdgeColor','k');
                set(h,'FaceColor','k');
                if t(end)>t(1)
                    axis([t(1),t(end),0,max_value]);
                end
                hylabel = ylabel(chroma_names(p,:));
                set(hylabel,'Rotation',0);
                set(hylabel,'HorizontalAlignment','right');
                set(hylabel,'VerticalAlignment','middle');
                set(gca,'YTick',[0 1]);
                if p~=1
                    set(gca,'XTickLabel',[]);
                end
                if p==1
                    %xlabel('Time in seconds');
                end
                hold on;
                h = plot(t_chroma,f_chroma_energy_distr(t1_chroma+1:t2_chroma+1,p));
                set(h,'LineWidth',2);
                set(h,'Color',[170 170 170]/255);
            end
            drawnow;
            if parameter.print == 1
                set(gcf,'PaperPosition',parameter.printPaperPositionPlot);
                print('-depsc2',strcat(parameter.printDir,parameter.printFileChromaCENS));
            end
        end
    end
end
