为什么我的2D FFT会出现全零行？

Question

我正在尝试复现一篇论文中的结果。

论文中提到：“在南纬40度以南的区域，使用空间和时间上的二维傅里叶变换（2D-FT）来描述模拟的流量变化的频谱。” - Lenton等人（2006年）

发布的图表显示了“二维傅里叶变换方差的对数”。

我尝试创建一个数组，这个数组包含了类似数据的季节性循环以及噪声。我把噪声定义为原始数组减去信号数组。

这是我用来绘制信号数组在纬度上平均的二维傅里叶变换的代码：

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]等）的值都很低。

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])，那么对于 k 在 {0, 1, ..., N-1} 中，Y[2k] = 2*X[k]，其他情况下为零。

所以，我建议你看看你的 signal_time_lon 是怎么处理的。问题可能就出在这里。