rapsd2d.m

The radially averaged power spectrum (RAPS) is the direction-independent mean spectrum, i.e. the ave

% Computes and plots radially averaged power spectral density (power
% spectrum) of image IMG.
%
% (C)E. Ruzanski, RCG, 2009
function raPsd2d(img)

%% Process image size information
[N M] = size(img);
maxDim = max(N,M);
halfDim = floor(maxDim/2) + 1; % Only consider one half of spectrum

%% Compute power spectrum
imgf = fftshift(fft2(img));
imgfp = (abs(imgf)/(N*M)).^2; % Normalize

%% Radially average power spectrum
[X Y] = meshgrid(-N/2:N/2-1, -N/2:N/2-1);

[theta rho] = cart2pol(X, Y); % Convert to polar coordinate axes
rho = round(rho);
i = cell(floor(N/2) + 1, 1);
for r = 0:N/2
    i{r + 1} = find(rho == r);
end
Pf = zeros(1, floor(N/2) + 1);
for r = 0:floor(N/2)
    Pf(1, r + 1) = nanmean( imgfp( i{r+1} ) );
end

%% Setup plot
fontSize = 14;
maxX = 10^(ceil(log10(halfDim)));
f1 = 1:(maxX/halfDim):maxX; % Set abscissa

% Find axes boundaries
xMin = 0; % No negative image dimension
xMax = ceil(log10(halfDim));
xRange = (xMax:-1:xMin); % Frequencies decrease L-to-R

yMin = floor(log10(min(Pf)));
yMax = ceil(log10(max(Pf)));
yRange = (yMin:yMax);

% Create plot axis labels
for i = 1:length(xRange)
    xRangeS = num2str(10^(xRange(i)));
    xCell(i) = cellstr(xRangeS);
end

for i = 1:length(yRange)
    yRangeS = num2str(yRange(i));
    yCell(i) = strcat('10e',cellstr(yRangeS));
end

%% Generate plot
figure
loglog(f1,Pf,'k-','LineWidth',2.0)
set(gcf,'color','white')
set(gca,'FontSize',fontSize,'FontWeight','bold',...
    'YTickLabel',yCell,...
    'XTickLabel',xCell,'XGrid','on','YAxisLocation','right');
xlabel('Wavelength (km)','FontSize',fontSize,'FontWeight','Bold');
ylabel('Power','FontSize',fontSize,'FontWeight','Bold');
title('Radially averaged power spectrum of reflectivity image','FontSize',fontSize,'FontWeight','Bold')