我正试图复制纸张的结果。为什么我在2D fft中得到零行?
“沿着恒定纬度(东西)和经度(南北)的空间和时间的二维傅里叶变换(2D-FT)被用来表征模拟的通量变化的频谱在40degS以南“。 - 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)
所述英尺的每一个第二排完全由的零。为什么是这样? 将一个随机小数字添加到信号中以避免这种情况是可以接受的。
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]等)的值非常低。
谢谢。所以问题在于我故意构建了一个重复十次的时间序列。因此,对应于时间序列上十个周期的正弦曲线的log_ps [10,:]将非常亮/高,其他所有像素将接近零或零。对于{0,1,...,N-1}中的k,Y [10k] = 10 * K [k] – nicholaschris
在你的例子中没有高频变化。但在我的数据中应该有高频率变化,这并没有真正显示。或者,与低频变化(十个周期平铺时间序列)相比,高频变化性非常小。 – nicholaschris
是的,高频组件的能量可能很小,而许多真实世界的信号通常都是这样。尝试绘制一个对数刻度值(不是log(1 + x)刻度)。或者改变你的图的vmin/vmax参数。 –