How to plot the power spectrum

Question

I have an array of data whose columns are solution vectors to a system of ODEs at a specific time. I want to plot the power spectrum of a solution at a specific time, but when I attempt this I get unexpected results. I want to make sure I am plotting correctly.

First of all, I follow Matlab's format for storing the spatial frequency values $k$:

k=[0:L/2-1 -L/2:-1]*2*pi/L where L is the total number of grid points. I do NOT understand why the ks are stored this way, but I digress.

In the following snippit of code, I pick a column, the sixth, take the FFT, then take the absolute-value squared and plot this against k in the following way:

plot(k, abs( fft( udata(:,6) ) )'.^2 )

where udata is the array of solution vectors.

For certain parameters of my code, there is a value of k for which I expect the power-spectrum to be the most extreme. However, I do not see this reflected in my plotting.

For example, in the below code, the value for which I expect the power-spectrum to be maximized at is k=0.336311. However, in the plot below, you can see that I do not observe a big peak there. In fact, there is a peak at 0 which I do not understand.

%% Strang splitting for solving the fDNLS
% i ∂_t u_n = (L_alpha u)_n - |u_n|^2 u_n

% Grid and time information

% Integer lattice
L = 100; 
dt=0.01; tmax=25; % Time step and max time
nmax=round(tmax/dt); % Maximum cycles related to max time
dx=1; % Unit grid spacing
x=[-L/2:dx:L/2-dx]'; % Grid
k=2*pi*[0:L/2-1 -L/2:-1]'/L;
[row col] = size(x);
%% Parameters
A = 1;
alpha =1;
eps=1;
kmax=0.336311;

% Initial condition
%u = sin(x).^2; 
initial_condition = A + exp(1i*kmax*x);
%initial_condition(1:50) = 0;
u = initial_condition; 
%u = A*ones(row,col) + ones(row,col).*cos(pi*x);    
%u = A*ones(row,col);
udata=u; tdata=0;

% Long Range Interaction matrix
LRI_MATRIX = LRI_Matrix_periodic_direct(L/2,alpha); 

for nn=1:nmax
    % Cycle through the dense part of the splitting operators
    % Let v denote the intermediate solution 
    if nn == 1
        % First cycle

        % Half-step of Linear operator 
        v = expm(-1i * eps^(-(1+alpha)) *LRI_MATRIX * dt/2)*u; 
        
        % Full-step of Nonlinear operator
        u = v.*exp(1i * v .* conj(v) * dt); 

    elseif nn ~= 1 && nn ~= nmax
        % Middle cycles
        
        % Full-step of Linear operator
        v =  expm(-1i * eps^(-(1+alpha)) *LRI_MATRIX * dt)*u;  
        % Full-step of Nonlinear operator
        u = v.*exp(1i * v .* conj(v) * dt);
        
    elseif nn == nmax
        % Last cycle
        
        % Half-step of Linear operator
        u = expm(-1i * eps^(-(1+alpha)) * LRI_MATRIX * dt/2)*u; 
    end

   % Store the cycle data
   udata=[udata u]; tdata=[tdata nn*dt];
end


figure
plot(k, abs(    fft(   udata(:,6)   )  )'.^2 )
xlabel("Spatial frequency"), ylabel('$|\hat{u}|^2$','Interpreter','latex'),
title("Power spectrum")

%Long range interaction on the lattice [-N,N] with periodic BC on n in {-N, ..., N-1}
%p = dim = 2N
function M = LRI_Matrix_periodic_direct(N,a)

p = 2*N;
%% Initialize the coefficient matrix for the LRI operator
M = zeros(p,p);

%% Populate the diagonal 
M_DIAG = 2*zeta(1+a)*(1 - 1/(2*N)^(1+a))*eye(p);

%% Populate the dense part of the matrix
M_DENSE = zeros(p);

for i = 1:p
% Go through each row
    for j = 1:p 
        % Go through each column
        if j~=i
            M_DENSE(i,j) = -cnr(i-(N+1),j-(N+1),a,N);
        end
    end
end

M = M_DIAG + M_DENSE;
end
function off_diag = cnr(n,r,a,N)
    off_diag = 1/abs(n-r)^(1+a) + 1/(2*N)^(1+a) * (hurwitzZeta(1+a, (r-n)/(2*N)) + hurwitzZeta(1+a, -(r-n)/(2*N)) - abs((r-n)/(2*N))^(-(1+a)) * (1 + exp(-1i * (1+a) * pi)));
end

Answer 1

I pick a column, the sixth, take the FFT, then take the absolute-value squared

I guess this is not one PROPER method to estimate power spectrum.

The code plot(k, abs( fft( udata(:,6) ) )'.^2 ) you provided uses the Fast Fourier Transform (FFT) to estimate the power spectrum of the signal, but you may need a more refined method.

---

FYI

Here are some commonly used methods for Spectral Density Estimation:

1. Periodogram: This is one of the simplest and most direct methods of estimating the power spectrum. The periodogram is the absolute square of the Fourier transform of a signal. It provides a way to identify the frequency components of a signal and the power at those frequencies.

2. Welch's Method: This method is an improvement over the periodogram method that reduces noise in the estimated power spectra in exchange for reducing the frequency resolution. In Welch's method, the longer time series is split into overlapping segments, each of which is detrended, windowed, Fourier transformed, and then the power spectra of these segments are averaged.

3. Bartlett's Method: Similar to Welch's method, Bartlett's method involves dividing the time series into segments and averaging the periodograms of these segments. The difference is that Bartlett's method does not use overlapping segments or windowing, and therefore, it can yield higher noise in the estimated power spectrum.

4. Blackman-Tukey Method: This method is based on estimating the power spectrum through the Fourier transform of the autocorrelation sequence. It involves windowing the autocorrelation sequence to reduce the effect of sidelobes.

5. Maximum Entropy Method (MEM): This method estimates the power spectrum based on an autoregressive model. The MEM aims to find the model that is most consistent with the data, under the constraint of having the maximum entropy.