计算多维数组中每个数组的每列的傅里叶变换

Question

在下面的4D数组中，每一列代表机器学习模型开发的一个属性。

import numpy as np
from scipy.fft import fft, fftfreq, fftshift
A = np.array([
[[[0, 1, 2, 3],
[3, 0, 1, 2],
[2, 3, 0, 1],
[1, 3, 2, 1],
[1, 2, 3, 0]]],

[[[9, 8, 7, 6],
[5, 4, 3, 2],
[0, 9, 8, 3],
[1, 9, 2, 3],
[1, 0, -1, 2]]],

[[[0, 7, 1, 2],
[1, 2, 1, 0],
[0, 2, 0, 7],
[-1, 3, 0, 1],
[1, 0, 1, 0]]]
])
A.shape
(3, 1, 5, 4)

在本例中，A有3个实例，每个实例由4个特征表示(1个特征表示内部数组(形状为(1,5,4))的每列)。特征用time-domain表示。

我需要一种方便的方法来计算每个实例的每个特征的frequency-domain(傅里叶变换)。

考虑到给定的例子，我这样做。

1. 获取所有特性的数组(每个特性对应1)，如下所示:

   feat1 = A[..., [0]] # the first feature
   feat2 = A[..., [1]] # the second feature
   feat3 = A[..., [2]] # 3rd feature etc...
   feat4 = A[..., [3]] # 4th feature etc...
   

   因此示例数据中的第一个特征是:

   feat1
   array([[[[ 0],
   [ 3],
   [ 2],
   [ 1],
   [ 1]]],
   
   [[[ 9],
   [ 5],
   [ 0],
   [ 1],
   [ 1]]],
   
   [[[ 0],
   [ 1],
   [ 0],
   [-1],
   [ 1]]]])
   

2. 检索每个实例的特性集，如下所示:

   # 1st instance
   inst1_1 = feat1[0].ravel()  # instance1 1st feature
   inst1_2 = feat2[0].ravel()  # instance1 2nd feature
   inst1_3 = feat3[0].ravel()  # instance1 3rd feature
   inst1_4 = feat4[0].ravel()  # instance1 4th feature
   # 2nd instance 
   inst2_1 = feat1[1].ravel()  # instance2 1st feature
   inst2_2 = feat2[1].ravel()  # instance2 2nd feature
   inst2_3 = feat3[1].ravel()  # instance2 3rd feature
   inst2_4 = feat4[1].ravel()  # instance2 4th feature
   # 3rd instance 
   inst3_1 = feat1[2].ravel()  # instance3 1st feature
   inst3_2 = feat2[2].ravel()  # instance3 2nd feature
   inst3_3 = feat3[2].ravel()  # instance3 3rd feature
   inst3_4 = feat4[2].ravel()  # instance3 4th feature
   

3. 然后计算傅里叶变换(对于每个实例的每个特征)。

   inst1_1_fft = fft(inst1_1)
   inst1_2_fft = fft(inst1_2)
   inst1_3_fft = fft(inst1_3)
   inst1_4_fft = fft(inst1_4)
   inst2_1_fft = fft(inst2_1)
   inst2_2_fft = fft(inst2_2)
   inst2_3_fft = fft(inst2_3)
   inst2_4_fft = fft(inst2_4)
   inst3_1_fft = fft(inst3_1)
   inst3_2_fft = fft(inst3_2)
   inst3_3_fft = fft(inst3_3)
   inst3_4_fft = fft(inst3_4)
   

这绝对是dummies approach。此外，我不能在我的真实数据集上这样做，因为我有超过60K个实例。

问题:

1. 如何一次性为整个数组A执行此操作?

2. 对于每个实例的特征，频率似乎是相同的。我这样做对吗?我就是这样做的:

sampling_rate = A.shape[2] # each instance is fixed-sized A.shape[2] (5 in this e.g)
N = inst1_1_fft.size
inst1_1_fft_freq = fftfreq(N, d = 1 / sampling_rate)
inst1_1_fft_freq
array([ 0.,  1.,  2., -2., -1.])

Answer 1

为了一次沿着所有实例的特征执行快速傅里叶变换，您可以沿着axis=2执行。作为一种替代方案，您可以在计算fft之前，先通过使用np.reshape和np.swapaxes来重塑A，一种可能的解决方案是:

A.shape
>>> (3, 1, 5, 4)
N = 5 # number of time samples for the feature axis
inst_reshape = A.swapaxes(1,3).reshape((-1,N))
inst_reshape.shape
>>> (12, 5) # (number of instances x number of features, time sampling)

，它遵循与你的问题中相同的实例/特性顺序。您也可以通过将结果inst_reshape数组与手动组装的对应数组进行比较来检查这一点:

inst_vstack = np.vstack((inst1_1,inst1_2,inst1_3,inst1_4,
inst2_1,inst2_2,inst2_3,inst2_4,
inst3_1,inst3_2,inst3_3,inst3_4))
np.allclose(inst_vstack,inst_reshape)
>>> True

这样做之后，您现在可以沿着新重构数组的axis=1列轴调用np.fft.fft一次:

sampling_rate = A.shape[2]
sample_spacing = 1/sampling_rate
freq = np.linspace(0.0, 1.0/(2.0*sample_spacing), N//2)
X = np.fft.fft(inst_reshape, axis=1)
amp = 2/N*np.abs(X[:N//2])
from matplotlib import pyplot as plt
plt.plot(freq, amp)

我强烈建议你看一下下面的问题，其中很好地解释了这些数组转换背后的一般思想。