% IMF_TOIMAGE.M
%
% 
% IMF_TOIMAGE(A,f,Fs,scale,varargin) displays frequencies and amplitudes
%             of imf.
%
% inputs:	
%   - A(mandatory): amplitude matrix
%   - f(mandatory): frequency matrix
%   - Fs(mandatory): sampling frequency
%   - scale(mandatory): linear or logaritmic scale
%                       'lin' for linear
%                       'log' for logaritmic
%   - threshold (optional): amplitude threshold.
%                           amplitedes below threshold won't be displayed,
%                           measured in percentage of maximum amplitude.
%                           ignored if set to 0.
%   - [fmin fmax](optional): minimum and maximum frequency that will be
%                            displayed, ignored if set to [0 0].
%   - [tmin tmax](optional): minimum and maximum sample that will be
%                            displayed
%
%examples:
%
% displays all frequencies in linear scale
%>>imf_toimage(A,f,44100,'lin');
%
% displays frequencies from 100 Hz to 1000 Hz in linear scale
%>>imf_toimage(A,f,44100,'lin', 0, [100 1000]);
%
% displays frequencies from 100 Hz to 1000 Hz in logaritmic scale
%>>imf_toimage(A,f,44100,'log', 0, [100 1000]);
%
% displays frequencies(100Hz, 1000Hz) in logaritmic scale from
% sample 1000 to sample 2000.
%>>imf_toimage(A,f,44100,'log', 0, [100 1000], [1000 2000]);
%
% displays all frequencies in linear scale from sample 1000 to
% sample 2000.
%>>imf_toimage(A,f,44100, 'lin', 0, [0 0], [100 2000]);
%
% displays all frequencies in linear scale from sample 1000 to
% sample 2000 for amplitudes above 50% of maximum amplitude.
%>>imf_toimage(A,f,44100, 'lin', 0.5, [0 0], [100 2000]);
%


function imf_toimage(A,f,Fs,scale,varargin)
res_x = 1024;
res_y = 768;
threshold = 0; % amplitudes below threshold(%) won't be displayed
im = zeros(res_y,res_x);

f = f*Fs;

if(nargin > 4)
  threshold = varargin{1};
  
  if(length(threshold)>1)
     error('Invalid threshold');
  end
end

Amax = max(max(A));
indA = find(A >= threshold*Amax);
fdisp = f(indA);

fmax = max(max(fdisp));
fmin = min(min(fdisp));
tmin = 1;
tmax = length(f);
t = 1:length(f);
f_range = [];


if(nargin < 4)
  error('Incorrect input arguments.');
end

if (nargin > 3)
  if (scale == 'log')
    log = 1;
  else
    log = 0;
  end
end

if (nargin > 5)
  f_range = varargin{2};
  if any(f_range)
    fmin = f_range(1);
    fmax = f_range(2);
    
    if(fmin > max(max(f)) || fmax < min(min(f)))
      error('Frequency out of displaying range.');
    end
  end
end
if (nargin > 6)
  t_range = varargin{3};
  tmin = t_range(1);
  tmax = t_range(2);
  if (~any(f_range) && threshold == 0)
    fmax = max(max(f(:,tmin:tmax)));
    fmin = min(min(f(:,tmin:tmax)));
  end
end


step_x = (res_x-1)/(tmax - tmin);

if(log)
 step_y = (res_y-1)/(log10(fmax) - log10(fmin));
else
 step_y = (res_y-1)/(fmax - fmin);
end

for i=1:size(f,1)
    for j=tmin:tmax
      if(A(i,j)>threshold*Amax)
        if(f(i,j)>= fmin && f(i,j) <= fmax)   
          x = floor(step_x*(j - tmin)) + 1;
          if(log)
            y = floor(step_y*(log10(f(i,j))-log10(fmin))) + 1;
          else
            y = floor(step_y*(f(i,j)-fmin)) + 1;
          end
          im(y,x)=im(y,x) + A(i,j);
        end
      end
    end
end

if(log)
 correction = (fmax-fmin)/(2*(res_y));
else
 correction = 0;
end


h = imagesc([ceil(tmin):ceil(tmax)],[fmin + correction , fmax - correction], im);grid;
a = get(h,'Parent');
if(log)
 set(a, 'YScale', 'log');
 ticks = get(a, 'YTick');
 min_tick = ticks(1)/10;
 max_tick = ticks(end)*10;
 curr_tick = 0;
 i = 0;
 while (curr_tick < max_tick)
  while (curr_tick + min_tick/2 < min_tick*10)
    i = i+ 1;
    curr_tick = curr_tick + min_tick;
    freq(i) = curr_tick;
  end
  min_tick = min_tick*10;
 end
end
set(a, 'YDir', 'normal');
if(log)
  set(a, 'YTick', freq);
end




