% SCALOGRAM - Calculates phase and amplitude scalogram of 1D signal.
%
% Usage: 
% [amplitude, phase] = scalogram(signal, minwavelength, mult, nscales,...
%                                sigmaOnf, map, threeD)
%
% Function to calculate the phase and amplitude scalograms of a 1D signal
% Analysis is done using quadrature pairs of log Gabor filters
%
% Arguments:
%        signal        - a 1D vector to be analyzed (for maximum speed 
%                        length should be a power of 2)
%        minwavelength - wavelength of smallest scale filter to use
%        mult          - scaling factor between successive filters
%        nscales       - No of filtering scales to use
%        sigmaOnf      - Shape factor of log Gabor filter controlling bandwidth
%                       .35 - bandwidth of approx 3 ocatves
%                       .55 - bandwidth of approx 2 ocatves
%                       .75 - bandwidth of approx 1 ocatve
%        map           - Optional specificalion of (cyclic) colour map. If
%                        not supplied or is empty a white -blue - white red
%                        cyclic map is used (cmap('C4'))
%        threeD        - An optional argument (0 or 1) indicating whether
%                        to display additional 3D visualisations of the
%                        result (this can take a while to display).
%        displayPlots  - Option flag 0/1 indicating whether to display plots
%                        or not.  Defaults to 1
%
% Output:
%         amplitude    - image of amplitude responses
%         phase        - image of phase responses
%
% Suggested values:
% [amplitude, phase] = scalogram(signal, 4, 1.05, 128, .55)
%

% Copyright (c) 1997-2014 Peter Kovesi
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in 
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% March     1997  - Original version
% September 2001  - Code tidied and plots improved.
% January   2013  - More Octave compatibility (Carnë Draug)
% February  2014  - Ensure input signal is double! This avoids problems with
%                   MATLAB's stupid convention of casting data types *down*
%                   to the lowest class when they are combined.
% October  2014   - Option to supply cyclic colour map
% January  2015   - Some modifications to plotting and plotting options.

function [amplitude, phase, tickLocations, tickValues] = ...
        scalogram(signal, minwavelength, mult, nscales, sigmaOnf, map, threeD, ...
                  displayPlots)

    Octave = exist('OCTAVE_VERSION', 'builtin') == 5; % Are we running under Octave
    
    assert(isvector(signal), 'Input vector must be 1D');
    signal = signal(:)';     % Ensure it is a row vector
    if ~isa(signal,'double'), signal = double(signal); end
    if ~exist('map', 'var') || isempty(map)
        map = cmap('C4','s',.75); 
    end 
    if ~exist('threeD', 'var'), threeD = 0; end 
    if ~exist('displayPlots', 'var'), displayPlots = 1; end 


    ndata = length(signal);
    
    % Trim off the last signal element if the array length is odd.  The basic
    % filter construction code assumes an even number of data elements.
    if mod(ndata,2), signal = signal(1:end-1); ndata = ndata-1; end
    
    signalfft = fft(signal);           % Take FFT of signal
    
    amplitude = zeros(nscales,ndata);  % Pre-allocate memory for speed
    phase     = zeros(nscales,ndata);
    logGabor  = zeros(1,ndata);
    EO        = zeros(1,ndata);
    
    radius =  [0:fix(ndata/2)]/fix(ndata/2)/2;  % Frequency values 0 - 0.5
    radius(1) = 1;                     % Fudge to stop log function complaining at origin.

    wavelength = minwavelength;

    for row = 1:nscales
        fo = 1.0/wavelength;            % Centre frequency of filter.
        
        % Construct filter
        logGabor(1:ndata/2+1) = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2));  
        logGabor(1) = 0;                % Set value at zero frequency to 0 (undo fudge).

        % Multiply filter and FFT of signal, then take inverse FFT.
        EO = ifft(signalfft .* logGabor);
        
        amplitude(row,:) = abs(EO);     % Record the amplitude of the result
        phase(row,:) = angle(EO);       % .. and the phase.
        
        wavelength = wavelength * mult; % Increment the filter wavelength.
    end 


    % Set up axis range for plotting the results.
    minsig = min(signal); maxsig = max(signal); r = maxsig-minsig;  
    rmin = minsig-r/10; rmax = maxsig+r/10;

    % Set up data for labelling the scale axis
    Nticks = 8;
    dtick = fix(nscales/Nticks);
    tickLocations = 1:dtick:nscales;
    tickValues = round(minwavelength*mult.^(tickLocations-1));

    if displayPlots
        
    % Amplitude only scalogram.
    figure(1), clf, colormap(gray)
    subplot(2,1,1), plot(signal), axis([0,ndata,rmin,rmax]), title('signal');
    h = subplot(2,1,2); imagesc(-amplitude), colormap(gray(256))
    title('amplitude scalogram');
    set(h,'YTick',tickLocations);
    if ~Octave
        set(h,'YTickLabel',tickValues);
    end
    
    ylabel('scale');

    % Phase only scalogram: phase encoded by colour map, saturation uniform.
    figure(2), clf

    subplot(2,1,1), plot(signal), axis([0,ndata,rmin,rmax]), title('signal');
    h = subplot(2,1,2); 
    image(showangularim(phase, map, 'cycle', 2*pi, 'bw', 1, 'fig', 0))
    title('phase scalogram');
    set(h,'YTick',tickLocations);
    if ~Octave
        set(h,'YTickLabel',tickValues);
    end
    
    ylabel('scale');
    
    % Phase and amplitude scalogram: phase encoded by colour map, amplitude
    % encoded by saturation.   Note we use the log of the amplitude to
    % to give a 'flatter' modulation of saturation 
    figure(3), clf
    amp = log1p(amplitude);
    subplot(2,1,1), plot(signal),  axis([0,ndata,rmin,rmax]), title('signal');
    h = subplot(2,1,2); 
    image(showangularim(phase, map, 'amp', amp, 'cycle', 2*pi, 'bw', 1, ...
                        'fig',0))
    
    title('phase-amplitude scalogram');
    set(h,'YTick',tickLocations);
    if ~Octave
        set(h,'YTickLabel',tickValues);
    end

    ylabel('scale');

    
    % Figure showing phase and phase-amplitude scalograms
    figure(10), clf
    
    subplot(2,2,1), plot(signal), axis([0,ndata,rmin,rmax]), title('signal');
    
    h = subplot(2,2,2); 
    image(showangularim(phase, map, 'cycle', 2*pi, 'bw', 1, 'fig', 0))
    title('phase scalogram');
    set(h,'YTick',tickLocations);
    ylabel('scale');

    h = subplot(2,2,3); imagesc(-amplitude), colormap(gray(256))
    title('amplitude scalogram');
    set(h,'YTick',tickLocations);
    ylabel('scale');    

    h = subplot(2,2,4); 
    image(showangularim(phase, map, 'amp', amplitude, 'cycle', 2*pi, ...
                        'bw', 1, 'fig', 0))
    title('phase-amplitude scalogram');
    set(h,'YTick',tickLocations);
    ylabel('scale');

%{
    figure(20), clf
    plot(signal),  axis([0,ndata,rmin,rmax]), title('signal');
    ylabel(' ')

    figure(21), clf
    amp = log1p(amplitude);    
    image(showangularim(phase, map, 'amp', amplitude, 'cycle', 2*pi, 'bw', 1, ...
                        'fig',0))
    
    title('phase-amplitude scalogram');
    set(h,'YTick',tickLocations);
    ylabel('scale');
%}

    
    
    if threeD            % Generate 3D plots.
        h = figure(4); clf;
        surfl(amplitude, [30,60]), axis('ij'), view(30,70), box on
        axis([0,ndata,0,nscales,0,max(max(amplitude))]), shading interp, colormap(gray);
        title('amplitude surface'); 
        h = get(h,'CurrentAxes');
        set(h,'YTick',tickLocations);
        if ~Octave
            set(h,'YTickLabel',tickValues);
        end
        ylabel('scale');
        
        h = figure(5); clf;
        hsv(:,:,2) = ones(size(amplitude));           % saturation fixed at 1
        warp(amplitude, hsv2rgb(hsv)), axis('ij'), view(30,70), box on, grid on
        title('amplitude surface, phase encoded by hue');
        h = get(h,'CurrentAxes');
        set(h,'YTick',tickLocations);
        if ~Octave
            set(h,'YTickLabel',tickValues);
        end
        ylabel('scale');
    end

    end % if displayPlots