No BSD License  

Highlights from
Adaptive Optimal Kernel

from Adaptive Optimal Kernel by Tony Reina
An Adaptive Optimal-Kernel Time-Frequency Representation

polafint(nrad,nphi,ntheta,maxrad,nlag,plag,ptheta,rectafm2)
function polafm2 = polafint(nrad,nphi,ntheta,maxrad,nlag,plag,ptheta,rectafm2)
%  polafint: interpolate AF on polar grid;
% G.A. Reina 16 Jan 2007
% Modified from the C code provided by D. L. Jones and R. G. Baraniuk
% "An Adaptive Optimal-Kernel Time-Frequency Representation"
%   by D. L. Jones and R. G. Baraniuk, IEEE Transactions on Signal 
%   Processing, Vol. 43, No. 10, pp. 2361--2371, October 1995.

half_nphi = floor(nphi / 2);

polafm2 = zeros(max(maxrad), nphi);

for ii = 0:(nphi-1), % for all phi ...
    for jj = 0:(maxrad(ii+1)-1),	% and all radii with |theta| < pi

        ilag = floor(plag(jj+1, ii+1));
        rlag = plag(jj+1, ii+1) - ilag;

        if (ilag <= (nlag - 2))
            
            itheta = floor(ptheta(jj+1, ii+1));
            rtheta = ptheta(jj+1, ii+1) - itheta;

            % Some sort of wrap around going on here. Not quite sure why
            if (ii == (half_nphi))
                
                rtemp =  (rectafm2(itheta+2, ilag+1) ...
                    - rectafm2((nphi - (itheta + 1)), (ilag+1)))*rtheta ...
                    + rectafm2((nphi - (itheta + 1)), (ilag+1));

                rtemp1 =  (rectafm2(itheta+2, ilag+2) ...
                    - rectafm2((nphi - (itheta + 1)), (ilag+2)))*rtheta ...
                    + rectafm2((nphi - (itheta + 1)), (ilag+2));

                polafm2(jj+1, ii+1) = (rtemp1-rtemp)*rlag + rtemp;

            else

                itheta1 = itheta + 1;
                if ( itheta1 >= ntheta )
                    itheta1 = 0;
                end

                rtemp =  (rectafm2(itheta1+1, ilag+1) ...
                    - rectafm2(itheta+1, ilag+1))*rtheta ...
                    + rectafm2(itheta+1, ilag+1);

                rtemp1 =  (rectafm2(itheta1+1, ilag+2)...
                    - rectafm2(itheta+1, ilag+2))*rtheta ...
                    + rectafm2(itheta+1, ilag+2);

                polafm2(jj+1, ii+1) = (rtemp1-rtemp)*rlag + rtemp;
            end

        end
    end
end

Contact us at files@mathworks.com