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

조회 수: 2(최근 30일)
Priyank Pathak
Priyank Pathak 2022년 7월 13일
편집: dpb 2022년 7월 13일
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');
%above link for data
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);
wn_k=linspace(0,0.03 ,number_psd); % wave number k
grid on
grid minor
ylabel('PSD (Log(Amplitude.^2))')


Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by