Почему я получаю ряды нулей в моем 2D FFT? - PullRequest
2 голосов
/ 16 ноября 2011

Я пытаюсь воспроизвести результаты из статьи.

"Двумерное преобразование Фурье (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] соответствует синусоиде с одним циклом по всему сигналу в строках начального изображения.

Вот начальное изображение с использованием следующего кода: Starting image

plt.pcolormesh(signal_time_lon); plt.colorbar(); plt.axis('tight')

Вот результат с использованием следующего кода: enter image description here

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] и т. Д.) Имеют очень низкие значения.

1 Ответ

3 голосов
/ 17 ноября 2011

Рассмотрим следующий пример:

In [59]: from scipy import absolute, fft

In [60]: absolute(fft([1,2,3,4]))
Out[60]: array([ 10.        ,   2.82842712,   2.        ,   2.82842712])

In [61]: absolute(fft([1,2,3,4, 1,2,3,4]))
Out[61]: 
array([ 20.        ,   0.        ,   5.65685425,   0.        ,
         4.        ,   0.        ,   5.65685425,   0.        ])

In [62]: absolute(fft([1,2,3,4, 1,2,3,4, 1,2,3,4]))
Out[62]: 
array([ 30.        ,   0.        ,   0.        ,   8.48528137,
         0.        ,   0.        ,   6.        ,   0.        ,
         0.        ,   8.48528137,   0.        ,   0.        ])

Если X[k] = fft(x) и Y[k] = fft([x x]), то Y[2k] = 2*X[k] для k in {0, 1, ..., N-1} и ноль в противном случае.

Поэтому я бы посмотрелв то, как ваш signal_time_lon выложен плиткой.Возможно, именно в этом проблема.

...