% 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 2004 function [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');