% ODOT - Demonstrates odot and oslash operators on 1D signal%% Usage:  [smooth, energy] = odot(f, K)%% Arguments:    f - a 1D signal%               K - optional 'Weiner' type factor to condition the results%                   where division by 0 occurs in the 'oslash' operation.%                   K defaults to 'eps', If oscillations appear in the%                   plots try increasing the value of K%% Returns:     energy - the Local Energy of the signal.%              smooth - the smooth component of the signal obtained by%                       performing the 'oslash' operator between the%                       signal and its Local Energy.% Plots:%           Signal            Hilbert Transform of Signal%           Local Energy      Hilbert Transform of Energy%           Smooth Component  Reconstruction%% Smooth         = signal 'oslash' energy% Reconstruction = energy  'odot'  smooth%% Glitches in the results will be seen at points where the Local Energy% peaks - these points cause numerical grief.  These problems can be% alleviated by smoothing the signal slightly and/or increasing the% parameter K.%% This code only works for 1D signals - I am not sure how you would% implement it for 2D images...% Reference:  Robyn Owens.  "Feature-Free Images", Pattern Recognition% Letters. Vol 15. pp 35-44, 1994.% Copyright (c) 2004 Peter Kovesi% School of Computer Science & Software Engineering% The University of Western Australia% http://www.csse.uwa.edu.au/% % 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 2004function [smooth, energy] = odot(f, K)        if nargin == 1	K = eps;    end        N = length(f);        if rem(N,2) > 0  % odd No of elements - trim the last one off to make the		     % number of elements even for simplicity.	N = N - 1;	f = f(1:N);    end        F = fft(f);        F(1) = 0;                        % Kill DC component    f = real(ifft(F));               % Reconstruct signal minus DC component        % Perform 90 degree phase shift on the signal by multiplying +ve    % frequencies of the fft by i and the -ve frequencies by -i, and then    % inverting.         phaseshift = [ ones(1,N/2)*i ones(1,N/2)*(-i) ];        % Hilbert Transform of signal    h = real(ifft(F.*phaseshift));                energy = sqrt(f.^2 + h.^2);           % Energy    % Hilbert Transform of Energy    energyh = real(ifft(fft(energy).*phaseshift)); 					      % smooth = signal 'oslash' energy    divisor = energy.^2 + energyh.^2;       % Where divisor << K,  weinercorrector -> 0/K    % Where divisor >> K,  weinercorrector -> 1        weinercorrector = divisor.^2 ./ ((divisor.^2)+K);    smooth = (f.*energy + energyh.*h)./divisor .* weinercorrector;        % Hilbert transform of smooth component    smoothh=real(ifft(fft(smooth).*phaseshift));        % Reconstruction = energy odot smooth    recon = (smooth.*energy - smoothh.*energyh);             subplot(3,2,1), plot(f),title('Signal');     subplot(3,2,2), plot(h),title('Hilbert Transform of Signal');     subplot(3,2,3), plot(energy),title('Local Energy');     subplot(3,2,4), plot(energyh),title('Hilbert Transform of  Energy');     subplot(3,2,5), plot(smooth),title('Smooth Component');     subplot(3,2,6), plot(recon),title('Reconstruction');