How to transform data into wavenumber-frequency space?

Im trying to understand how to transform OLR netcdf data into wavenumber frequency space. I used before bandpass function, but that only works to extract only certain time-frequencies, as 20 days to 60 days window. But, how does this work for wavenumber domain? Which function does that?

Answer 1

[Sudarsanan A K]

Hello Luis,

To transform OLR (Outgoing Longwave Radiation) data from the time-space domain to the wavenumber-frequency domain in MATLAB, you can use the 3D Fast Fourier Transform (FFT). This will allow you to analyze the spatial variations of the OLR data in terms of wavenumbers.

Here's a step-by-step approach on how to perform this transformation:

1. Load the OLR data from the netCDF file into MATLAB using appropriate functions like "ncread" or "ncinfo".
2. Assuming your OLR data is a 3D field with dimensions (latitude, longitude, and time), you can calculate the spatial mean and remove it from the data. This step helps to remove the large-scale variations and focus on the smaller-scale features.
3. Apply a windowing function to reduce the spectral leakage that can occur during the Fourier transform. Common windowing functions include the Hann window ("hann" function) or the cosine taper ("tukeywin" function). Apply the chosen window function to the OLR data.
4. Perform 3D FFT (using n dimensional FFT function "fftn") to compute the Fourier transform of the windowed OLR data. This will yield the OLR data in wavenumber space.
5. Shift the zero-frequency component (DC component) to the center of the spectrum using the "fftshift" function. This step is necessary to correctly visualize the wavenumber spectrum.
6. Calculate the wavenumbers corresponding to the shifted spectrum using the the grid spacing of your OLR data.

Here's an example code snippet that demonstrates these steps:

% Load OLR data from netCDF file
olrData = ncread('olr_data.nc', 'olr');

% Calculate spatial mean and remove it
meanOLR = mean(olrData, [1, 2]);
olrData = olrData - meanOLR;

% Apply Hann window
windowedOLR = hann(size(olrData, 1)) * hann(size(olrData, 2))' .* olrData;

% Perform 3D FFT
fftOLR = fftshift(fftn(windowedOLR));

% Calculate wavenumbers
[nlat, nlon, ntime] = size(olrData);
dx = 360 / nlon; % Assuming a global grid
dy = 180 / nlat; % Assuming a global grid
dkx = 1 / (nlon * dx); % Wavenumber increment in x direction
dky = 1 / (nlat * dy); % Wavenumber increment in y direction
kx = fftshift((-nlon/2:nlon/2-1) * dkx);
ky = fftshift((-nlat/2:nlat/2-1) * dky);

% Plot the wavenumber spectrum at a specific time index (e.g., 1)
imagesc(kx, ky, abs(fftOLR(:,:,1))); % Use magnitude for visualization
colorbar;
xlabel('Zonal Wavenumber (kx)');
ylabel('Meridional Wavenumber (ky)');

This code calculates the 3D FFT of the windowed OLR data and plots the resulting wavenumber spectrum. The "k" variable represents the wavenumbers in both the zonal (x) and meridional (y) directions. Note that this code assumes a regular global grid with a longitude range of 0 to 360 degrees. Adjustments may be needed depending on the specific characteristics of your OLR data.

For better understanding of the "fftn" and "fftshif" functions, you can additionally refer to the MathWorks documentation in the links:

https://in.mathworks.com/help/matlab/ref/fftn.html

https://in.mathworks.com/help/matlab/ref/fftshift.html

I hope this helps you transform your OLR data into the wavenumber domain in MATLAB!

[Luis Jesús Olvera Lazcano]

Wow, thank you so much for the detailed explanation! Im gonna try it tomorrow so I'll let you know how it goes. Nevertheless Im sorry for taking so long to answer

You just helped me a lot with this!

[Ya]

Hello Sudarsanan,

Thank you so much for the elaboration and script for the wavenumber-frequency spectrum.

1) For the frequency scale, woud you use:

f = fftshift((-ntime/2:ntime/2-1) * df);

where df = 1/(ntime*sampling_rate)

This would mean the zeor frequnecy is at the centre and one would read from the middle of the freq scale?

2) To obtain the correct physical quantity, i.e. magnitude, of the fft, would the following to correct?

mag(2:( ntime/2-1)) = abs(fftOLR(2:(ntime/2-1))) ./ (ntime/2);

3) Could you elaborate on how to obtain the single sided Power Spectrum Desnity fucntion form the fftn reuslts?

PSD = (1/((1/dt)*ntime)) * abs(fftOLR).^2

Kind regards,

Ya