%% energyhist.m % 2D energy histograms % compute and display distribution of energy by radial frequency % S. Allen Broughton - 1 Oct 10 %% get and display image and shifted DFT load mandrill % load earth % load earth; X = X-mean(X(:)); % eliminate the DC coefficient % X = zeros(128,128); X(61:80,41:80)=1; % rectangle fX = fft2(X); sfX = fftshift(fX); figure(1) subplot(1,2,1) imagesc(X), colormap(gray), axis equal, axis tight title('original image'); subplot(1,2,2) imagesc(log(1+abs(sfX))), colormap(gray), axis equal, axis tight title('shifted log(1+|dft|)'); %% get image data and set up shifted frequency matrices [m,n] = size(X); N = m*n; hm = round(m/2); hn = round(n/2); maxfreq = sqrt(hm^2+hn^2); V = ((0:(m-1))-hm)'*ones(1,n); % vertical frequencies H = ones(m,1)*((0:(n-1))-hn); % horizontal frequencies R = sqrt(V.^2+H.^2); % radial frequencies %% display frequency maps figure(2) subplot(1,3,1) imagesc([-hm,hm],[-hm,hm],V), colormap(gray), axis equal, axis tight title('vertical frequencies'); subplot(1,3,2) imagesc([-hm,hm],[-hm,hm],H), colormap(gray), axis equal, axis tight title('horizontal frequencies') subplot(1,3,3) imagesc([-hm,hm],[-hm,hm],R), colormap(gray), axis equal, axis tight title('radial frequencies') %% set up freqency data and calculate energy histogram nbins = 50; cEhist = [ ]; df = maxfreq/nbins; E = energy(X) % display this to indicate calulation is proceeding for j = 1: nbins, J = find (R <= j*df); Y = zeros(size(sfX)); Y(J) = sfX(J); cEhist = [cEhist, energy(Y)/N]; end Ehist = cEhist -[0,cEhist(1:nbins-1)]; % show energies %% display results figure(3) freqs = (1: nbins)*df; subplot(2,1,1) stem(freqs,Ehist), xlabel('frequency') ylabel('energy') title('frequency-energy distribution') subplot(2,1,2) plot(freqs,cEhist,'o',freqs,E*ones(size(freqs)),'k-',freqs,0,'k-') % put in zero to include axis xlabel('frequency') ylabel('energy') title('frequency-energy cumulative distribution')