CODEWITHSUNDEEP
matlabmatrix
Taitai
Taitai

Reputation: 616

how to calculate the spectral density of a matrix of data use matlab

I am not doing signal processing. But in my area, I will use the spectral density of a matrix of data. I get quite confused at a very detailed level.

%matrix H is given.
corr=xcorr2(H);  %get the correlation 
spec=fft2(corr); % Wiener-Khinchin Theorem

In matlab, xcorr2 will calculate the correlation function of this matrix. The lag will range from -N+1 to N-1. So if size of matrix H is N by N, then size of corr will be 2N-1 by 2N-1. For discretized data, I should use corr or half of corr?

Another problem is I think Wiener-Khinchin Theorem is basically for continuous function. I have always thought that Discretized FT is an approximation to Continuous FT, or you can say it is a tool to calculate Continuous FT. If you use matlab build in function 'fft', you should divide the final result by \delta x.

Any kind soul who knows this area well there to share some matlab code with me?

Upvotes: 1

Views: 1174

Answers (1)

BillBokeey
BillBokeey

Reputation: 3511

Basically, approximating a continuous FT by a Discretized FT is the same as approximating an integral by a finite sum.

We will first discuss the 1D case, then we'll discuss the 2D case.

Let's look at the Wiener-Kinchin theorem (for example here).

It states that :

"For the discrete-time case, the power spectral density of the function with discrete values x[n], is :

where

Is the autocorrelation function of x[n]."

1) You can see already that the sum is taken from -infty to +infty in the calculation of S(f)

2) Now considering the Matlab fft - You can see (command 'edit fft' in Matlab), that it is defined as :

 X(k) =       sum_{n=1}^N  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.

which is exactly what you want to be done in order to calculate the power spectral density for a frequency f.

Note that, for continuous functions, S(f) will be a continuous function. For Discretized function, S(f) will be discrete.

Now that we know all that, it can easily be extended to the 2D case. Indeed, the structure of fft2 matches the structure of the right hand side of the Wiener-Kinchin Theorem for the 2D case.

Though, it will be necessary to divide your result by NxM, where N is the number of sample points in x and M is the number of sample points in y.

Upvotes: 1

Related Questions