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

2 views (last 30 days)
Priyank Pathak on 13 Jul 2022
Edited: dpb on 13 Jul 2022
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.
% https://drive.matlab.com/files/plot_bouguer_20220713013040.txt
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))')