Problem with 2D fftshift

Hello, everyone,

I hope you can help me

i have a problem with 2d fftshift.

I have to get the noise power spectrum from 18 Rois. I use fftshift but in the center of the picture are higher frequencies that don't belong. as far as i know there are zero frequencies in the middle. can someone help me with this problem.

My Code:---------------------------------------------------------------------------------------------------

clear all

close all

S = load('matlab.mat');

b= cell2mat(struct2cell(S));

[l m n]=size(b);

mm =mean(b,3);

figure(1)

subplot(1,3,1)

imshow(mm,'DisplayRange',[]);title('Original Image');

for i=1:n

spe1(:,:,i)= fftshift(fftn(b(:,:,i)-mm));

spe(:,:,i)=(abs(spe1(:,:,i)).^2);

end

w = (mean(spe,3));

subplot(1,3,2)

imshow(w,'DisplayRange',[]);

title('DFT');

%Radial average

[XX, YY] = meshgrid( (1:m)-m/2, (1:m)-m/2);

R = sqrt(XX.^2 + YY.^2);

profile = [];

for i = 1:1:(m+1)

mask = (i-1<R & R<i+1);

values = w(mask);

profile(i)= mean( values(:) );

end

%Frequenz

frequency=(1/((m)*0.68));

sumb=0:frequency:(((m))*frequency);

subplot(1,3,3)

plot(sumb, profile, '-')

xlabel('Spatial Frquency mm^-1')

ylabel('Noise power spectrum HU^2 mm^2')

Answer 1

[Dev]

Hi @Latifa Bouguessaa,

I am assuming that you want to apply FFT shift to a 2D FFT to get the zero-frequency (DC) component at the centre of the spectrum. Also, the highest frequencies should be at the corners. Please find some relevant changes below to the code provided in the question which can help us achieve the same-

• Meshgrid for Frequency Mapping- The code provided in the question is not centered correctly for even-sized images. For even m, zero-frequency is between pixels (m/2, m/2) and (m/2+1, m/2+1), as mentioned in the code below-

[XX, YY] = meshgrid( (-m/2):(m/2-1), (-m/2):(m/2-1) );

• Radial Profile Binning- The binning in the code may not be robust. Instead, we can use the “accumarray” function in MATLAB for radial averaging, as follows-

R = sqrt(XX.^2 + YY.^2);
R = round(R); % integer radius
maxR = max(R(:));
profile = accumarray(R(:)+1, w(:), [maxR+1 1], @mean, NaN);

For more information on the usage of “accumarray” function, please refer to the documentation link below-https://www.mathworks.com/help/releases/R2024b/matlab/ref/accumarray.html

• Frequency Axis- The frequency axis calculation in the code provided can be modified. For a pixel size of 0.68 mm, the frequency step is:

dx = 0.68; % mm
freq_step = 1/(m*dx);
frequencies = (0:maxR) * freq_step;

After incorporating the above changes, we can visualise the frequencies as follows-

I hope the above explanation helps to solve the question.