如何让 numpy 数组的 FFT 工作？

Question

我正在读取 csv 文件的特定列作为 numpy 数组。当我尝试做这个数组的 fft 时，我得到了一个 NaN 数组。如何让 fft 工作？到目前为止，这是我所拥有的：

#!/usr/bin/env python
from __future__ import division

import numpy as np
from numpy import fft

import matplotlib.pyplot as plt


fileName = '/Users/Name/Documents/file.csv'

#read csv file
df = np.genfromtxt(fileName, dtype = float, delimiter = ',', names = True)

X = df['X'] #get X from file

rate = 1000. #rate of data collection in points per second
Hx = abs(fft.fft(X))
freqX = fft.fftfreq(len(Hx), 1/rate)

plt.plot(freqX,Hx) #plot freqX vs Hx

Answer 1

您的 csv 文件中可能存在一些缺失值。默认情况下，np.genfromtxt 会将缺失值替换为NaN。

如果数组中有任何NaNs 或Infs，则fft 将全部为NaNs 或Infs。

例如：

import numpy as np

x = [0.1, 0.2, np.nan, 0.4, 0.5]
print np.fft.fft(x)

我们会得到：

array([ nan +0.j,  nan+nanj,  nan+nanj,  nan+nanj,  nan+nanj])

但是，由于 FFT 对一系列规则间隔的值进行运算，因此从数组中删除非有限值比仅删除它们要复杂一些。

pandas 有几个专门的操作可以做到这一点，如果你愿意使用它（例如fillna）。不过，用“纯”的 numpy 做起来并不难。

首先，我假设您正在处理一系列连续的数据，因为您正在对这些值进行 FFT。在这种情况下，我们希望根据它们周围的值对 NaN 值进行插值。线性插值 (np.interp) 可能并非在所有情况下都是理想的，但它不是一个糟糕的默认选择：

例如：

import numpy as np

x = np.array([0.1, 0.2, np.nan, 0.4, 0.5])
xi = np.arange(len(x))

mask = np.isfinite(x)
xfiltered = np.interp(xi, xi[mask], x[mask])

我们会得到：

In [18]: xfiltered
Out[18]: array([ 0.1,  0.2,  0.3,  0.4,  0.5])

然后我们可以正常计算FFT：

In [19]: np.fft.fft(xfiltered)
Out[19]: 
array([ 1.50+0.j        , -0.25+0.34409548j, -0.25+0.08122992j,
       -0.25-0.08122992j, -0.25-0.34409548j])

...并得到一个有效的结果。

Answer 2

如果您的数据包含 NaN 值，则需要对它们进行插值。或者，您可以使用傅立叶方程计算频谱，其中np.sum 替换为np.nansum。使用这种方法，您不需要插入 NaN 值，尽管缺失数据的数量会影响频谱。更多的缺失数据会导致光谱噪声，从而导致光谱值不准确。

下面是一个 MWE 来说明这个概念，并用图表显示了结果。 MWE 说明了如何计算包含多个缺失值的简单参考信号的单边幅度谱。

#!/usr/bin/python
# Python code to plot amplitude spectrum of signal containing NaN values
# Python version 2.7.13

from __future__ import division
import numpy as np
import pylab as pl
import random

LW = 2 #line width
AC = 0.5 #alpha channel
pi = np.pi

def periodogramSS(inputsignal,fsamp):
 N = len(inputsignal)
 N_notnan = np.count_nonzero(~np.isnan(inputsignal))
 hr = fsamp/N #frequency resolution
 t = np.arange(0,N*Ts,Ts)
 #flow,fhih = -fsamp/2,(fsamp/2)+hr #Double-sided spectrum
 flow,fhih = 0,fsamp/2+hr #Single-sided spectrum
 #flow,fhih = hr,fsamp/2
 frange = np.arange(flow,fhih,hr)
 fN = len(frange)
 Aspec = np.zeros(fN)
 n = 0
 for f in frange:
  Aspec[n] = np.abs(np.nansum(inputsignal*np.exp(-2j*pi*f*t)))/N_notnan
  n+=1
 Aspec *= 2 #single-sided spectrum
 Aspec[0] /= 2 #DC component restored (i.e. halved)
 return (frange,Aspec)

#construct reference signal:
f1 = 10 #Hz
T = 1/f1
fs = 10*f1
Ts = 1/fs
t = np.arange(0,20*T,Ts)
DC = 3.0
x = DC + 1.5*np.cos(2*pi*f1*t)

#randomly delete values from signal x:
ndel = 10 #number of samples to replace with NaN
random.seed(0)
L = len(x)
randidx = random.sample(range(0,L),ndel)
for idx in randidx:
 x[idx] = np.nan

(fax,Aspectrum) = periodogramSS(x,fs)

fig1 = pl.figure(1,figsize=(6*3.13,4*3.13)) #full screen
pl.ion()

pl.subplot(211)
pl.plot(t, x, 'b.-', lw=LW, ms=2, label='ref', alpha=AC)

#mark NaN values:
for (t_,x_) in zip(t,x):
 if np.isnan(x_):
  pl.axvline(x=t_,color='g',alpha=AC,ls='-',lw=2)

pl.grid()
pl.xlabel('Time [s]')
pl.ylabel('Reference signal')

pl.subplot(212)
pl.stem(fax, Aspectrum, basefmt=' ', markerfmt='r.', linefmt='r-')
pl.grid()
pl.xlabel('Frequency [Hz]')
pl.ylabel('Amplitude spectrum')

fig1name = './signal.png'
print 'Saving Fig. 1 to:', fig1name
fig1.savefig(fig1name)

参考信号（真实）显示为蓝色，缺失值用绿色标记。单边幅度谱以红色显示。 10 Hz 处的直流分量和幅度值清晰可见。其他值是由于参考信号被缺失数据分解造成的。