python - 为什么我的 2D fft 中会出现几行零？

Question

我正在尝试复制论文的结果。

“沿恒定纬度(东西)和经度(南北)部分的空间和时间二维傅立叶变换 (2D-FT) 用于表征 40 度以南模拟通量变化的频谱。 ” - Lenton 等人(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)

ft 的每隔一行完全由零组成。为什么是这样？ 向信号中添加一个随机的小数字以避免这种情况是否可以接受？

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] 等)具有非常低的值。

Answer 1

考虑以下示例:

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] 对于 {0, 1, ..., N-1} 中的 k，否则为零。

因此，我会研究一下您的 signal_time_lon 是如何平铺的。这可能就是问题所在。