How to find a power spectral density (PSD) of any matrix?

I am here attached a code which I am using, but don't know it is right or wrong.

please someone check it, It is working right for an matrix?

help will be appreciated.

data_01= load('plot_bouguer_20220713013040.txt');
% https://drive.matlab.com/files/plot_bouguer_20220713013040.txt 
%above link for data

A02=data_01(:,3);
A03=reshape(A02,541,541);
A=A03';
bouin=A;
numrows=541; 
numcolumns=541; 

truncation=0.1 ;
% Calculate mean value of the gravity map 
fftbou=fft2(bouin);   %fft2 computes the 2-D FFT of a 2-D matrix (bou, in this case) 
meangravity=(fftbou(1,1)/(numrows*numcolumns)); %the first term of the 2-D fft array divided by the product of the number of rows and columns provides the mean of given matrix

% Demean data 
disp('Data sets demeaned')  %displays the text 'Data sets demeaned' on the screen 
bou=bouin-meangravity;      %input data 

%A cosine Tukey window with a truncation of 10% default is applied 
wrows = tukeywin(numrows,truncation);  %this computes a 1-D cosine Tukey window of the same length as the original matrix input rows and with the truncation defined by the variable 'truncation' 
wcolumns = tukeywin(numcolumns,truncation);  %this computes a 1-D cosine Tukey window of the same length as the original matrix input columns and with the truncation defined by the variable 'truncation' 

w2 =wrows(:)*wcolumns(:).';  %this generates a 2-D cosine Tukey window multipliying the 1-D windows 
bou=bou.*w2';                %the original gravity input matrix, previously demeaned, is multiplied by the cosine window 

mapabou=bou';   %the original input matrix after demeaning is transposed 
fftbou=fft2(mapabou);  %the 2-D FFT of the gravity input matrix is computed after demeaning 
spectrum=abs(fftbou(1:numrows/2, 1:numcolumns/2));  %this computes the amplitude spectrum 
st11=spectrum;                                       %amplitude spectrum
p_s=log10(st11'.^2);                                % power spectrum 
%st=max(p_s);                                         % mean of each rows  of amplitude spectrum
st=mean(p_s) ;
number_psd = length(st);
st1=reshape(st,number_psd,1);  

wn_k=linspace(0,0.03 ,number_psd);                     % wave number k
dwn_k=wn_k';
plot(wn_k,st1,'.','MarkerSize',15)
grid on
grid minor
xlabel('wavenumber(1/km)')
ylabel('PSD (Log(Amplitude.^2))')