Я пытаюсь воспроизвести результаты из статьи.
"Двумерное преобразование Фурье (2D-FT) в пространстве и времени вдоль участков постоянной широты (восток-запад) и долготы (север-юг) были использованы для характеристики спектра смоделированной изменчивости потока к югу от 40 градусов Цельсия ".- Лентон и др. (2006) Опубликованные цифры показывают «журнал отклонений 2D-FT».
Я попытался создать массив, состоящий из сезонного цикла схожих данных, а также шума,Я определил шум как исходный массив минус массив сигналов.
Вот код, который я использовал для построения 2D-FT массива сигналов, усредненного по широте:
import numpy as np
from numpy import ma
from matplotlib import pyplot as plt
from Scientific.IO.NetCDF import NetCDFFile
### input directory
indir = '/home/nicholas/data/'
### get the flux data which is in
### [time(5day ave for 10 years),latitude,longitude]
nc = NetCDFFile(indir + 'CFLX_2000_2009.nc','r')
cflux_southern_ocean = nc.variables['Cflx'][:,10:50,:]
cflux_southern_ocean = ma.masked_values(cflux_southern_ocean,1e+20) # mask land
nc.close()
cflux = cflux_southern_ocean*1e08 # change units of data from mmol/m^2/s
### create an array that consists of the seasonal signal fro each pixel
year_stack = np.split(cflux, 10, axis=0)
year_stack = np.array(year_stack)
signal_array = np.tile(np.mean(year_stack, axis=0), (10, 1, 1))
signal_array = ma.masked_where(signal_array > 1e20, signal_array) # need to mask
### average the array over latitude(or longitude)
signal_time_lon = ma.mean(signal_array, axis=1)
### do a 2D Fourier Transform of the time/space image
ft = np.fft.fft2(signal_time_lon)
mgft = np.abs(ft)
ps = mgft**2
log_ps = np.log(mgft)
log_mgft= np.log(mgft)
Каждый второй ряд футов полностью состоит из нулей.Почему это?Было бы приемлемо добавить случайное малое число к сигналу, чтобы избежать этого.
signal_time_lon = signal_time_lon + np.random.randint(0,9,size=(730, 182))*1e-05
РЕДАКТИРОВАТЬ: Добавление изображений и уточнение значения
Вывод rfft2 по-прежнему представляется сложным массивом,Использование fftshift смещает края изображения к центру;У меня все еще есть спектр мощности независимо.Я ожидаю, что причина, по которой я получаю ряды нулей, заключается в том, что я заново создал временные ряды для каждого пикселя.Пиксель ft [0, 0] содержит среднее значение сигнала.Таким образом, ft [1, 0] соответствует синусоиде с одним циклом по всему сигналу в строках начального изображения.
Вот начальное изображение с использованием следующего кода:
plt.pcolormesh(signal_time_lon); plt.colorbar(); plt.axis('tight')
Вот результат с использованием следующего кода:
ft = np.fft.rfft2(signal_time_lon)
mgft = np.abs(ft)
ps = mgft**2
log_ps = np.log1p(mgft)
plt.pcolormesh(log_ps); plt.colorbar(); plt.axis('tight')
На изображении может быть не ясно, но только каждая вторая строка содержит полностью нули.Каждый десятый пиксель (log_ps [10, 0]) является высоким значением.Другие пиксели (log_ps [2, 0], log_ps [4, 0] и т. Д.) Имеют очень низкие значения.