SciPy - FFTpack

SciPy - 傅里叶变换 (scipy.fft)

傅里叶变换（Fourier Transform）是一种强大的数学工具，可以将信号（通常是时域信号 time-domain signal）分解为其构成的频率分量。在频域（frequency domain）分析信号可以揭示在时域中隐藏的特性。它广泛应用于信号处理、图像分析、音频工程以及许多其他科学和工程领域。SciPy 通过 scipy.fft 模块（它是旧版 scipy.fftpack 的现代替代）提供了全面的傅里叶变换功能。

快速傅里叶变换 (FFT)

快速傅里叶变换（Fast Fourier Transform, FFT）是计算离散傅里叶变换（Discrete Fourier Transform, DFT）及其逆变换的高效算法。

一维离散傅里叶变换

DFT 将一个 N 个复数序列 $x_0, x_1, \ldots, x_{N-1}$ 转换为另一个 N 个复数序列 $X_0, X_1, \ldots, X_{N-1}$，后者通常称为频率 bin。scipy.fft 中的 fft() 函数计算 DFT，ifft() 函数计算逆 DFT（Inverse DFT）。

import numpy as np
from scipy.fft import fft, ifft

# Create a simple input array (time-domain signal)
x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])
print(f"Original signal x: {x}")

# Apply the FFT
y = fft(x)
print(f"FFT of x (frequency-domain signal y):\n{y}")

# Apply the Inverse FFT (IFFT)
y_inverted = ifft(y)
print(f"IFFT of y (should be close to original x):\n{y_inverted}")
# Note: IFFT result might have very small imaginary parts due to floating point precision.
# We can take the real part if the original signal was real.
print(f"Real part of IFFT of y: {np.real(y_inverted)}")

输出将是：

Original signal x: [ 1.   2.   1.  -1.   1.5]
FFT of x (frequency-domain signal y):
[ 4.50000000+0.j          2.08155948-1.65109876j -1.83155948+1.60822041j
 -1.83155948-1.60822041j  2.08155948+1.65109876j]
IFFT of y (should be close to original x):
[ 1.0+0.j   2.0+0.j   1.0+0.j  -1.0+0.j   1.5+0.j ]
Real part of IFFT of y: [ 1.   2.   1.  -1.   1.5]

我们来看一个更实际的例子：一个带噪声的正弦波。

import numpy as np
from scipy.fft import fft, fftfreq
import matplotlib.pyplot as plt # For plotting (optional)

# Signal parameters
sampling_rate = 100  # Hz (samples per second)
duration = 5.0       # seconds
frequency = 2.0      # Hz (frequency of the sine wave)

num_samples = int(sampling_rate * duration)
time_vector = np.linspace(0.0, duration, num_samples, endpoint=False)

# Generate a sine wave with some noise
signal = np.sin(2 * np.pi * frequency * time_vector) + 0.5 * np.random.randn(num_samples)

print(f"Number of samples: {signal.size}")

# Compute the FFT
signal_fft = fft(signal)

# Generate the frequency bins for the x-axis of the FFT plot
# fftfreq(number_of_samples, sample_spacing)
sample_spacing = 1.0 / sampling_rate
frequency_bins = fftfreq(num_samples, d=sample_spacing)

# Plotting (optional, requires matplotlib)
# plt.figure(figsize=(12, 6))
# plt.subplot(1, 2, 1)
# plt.plot(time_vector, signal)
# plt.title("Time-Domain Signal")
# plt.xlabel("Time (s)")
# plt.ylabel("Amplitude")

# plt.subplot(1, 2, 2)
# # We plot the magnitude of the FFT (abs value) and usually only the positive frequencies
# positive_freq_indices = frequency_bins >= 0
# plt.plot(frequency_bins[positive_freq_indices], np.abs(signal_fft[positive_freq_indices]))
# plt.title("Frequency-Domain Signal (FFT Magnitude)")
# plt.xlabel("Frequency (Hz)")
# plt.ylabel("Magnitude")
# plt.xlim(0, sampling_rate / 2) # Show up to Nyquist frequency
# plt.grid(True)
# plt.tight_layout()
# plt.show()

# For text output, let's show some FFT values and corresponding frequencies
print("\nSample FFT magnitudes and their frequencies:")
for i in range(5):
    print(f"Freq: {frequency_bins[i]:.2f} Hz, Magnitude: {np.abs(signal_fft[i]):.2f}")
# Find the dominant frequency component
dominant_freq_idx = np.argmax(np.abs(signal_fft[frequency_bins >= 0])) # consider positive freqs
print(f"Dominant positive frequency: {frequency_bins[dominant_freq_idx]:.2f} Hz")

此代码生成一个带噪声的正弦波。fftfreq 函数计算与 FFT 输出中每个点对应的频率。FFT 输出 signal_fft 包含复数；np.abs(signal_fft) 给出它们的幅值（magnitude），表示每个频率分量的强度。对于实数输入信号，FFT 输出是对称的，因此我们通常只分析正频率。

一个示例文本输出（由于噪声，具体值会有所不同）：

Number of samples: 500

Sample FFT magnitudes and their frequencies:
Freq: 0.00 Hz, Magnitude: 1.89
Freq: 0.20 Hz, Magnitude: 6.09
Freq: 0.40 Hz, Magnitude: 6.37
Freq: 0.60 Hz, Magnitude: 3.90
Freq: 0.80 Hz, Magnitude: 10.47
Dominant positive frequency: 2.00 Hz

找到的主导频率应接近我们为正弦波设置的原始 frequency 2.0 Hz。

离散余弦变换 (DCT)

离散余弦变换（Discrete Cosine Transform, DCT）与 DFT 相关，但只使用余弦函数。它将一组有限的数据点表示为不同频率余弦函数的总和。DCT 广泛应用于信号和图像压缩（例如 JPEG）。

SciPy 的 scipy.fft 模块提供 dct 用于正向 DCT，idct 用于逆向 DCT (IDCT)。

from scipy.fft import dct, idct
import numpy as np

# Input array
x_dct = np.array([4., 3., 5., 10., 5., 3.])
print(f"Original array for DCT: {x_dct}")

# Compute DCT (default is Type-II DCT)
dct_coeffs = dct(x_dct)
print(f"DCT coefficients:\n{dct_coeffs}")

# Reconstruct the original signal using IDCT (default is Type-III IDCT for Type-II DCT)
x_reconstructed = idct(dct_coeffs)
print(f"Reconstructed array using IDCT:\n{x_reconstructed}")

输出将是：

Original array for DCT: [ 4.  3.  5. 10.  5.  3.]
DCT coefficients:
[ 60.          -3.48476592 -13.85640646  11.3137085   6.
  -6.31319305]
Reconstructed array using IDCT:
[ 4.  3.  5. 10.  5.  3.]

DCT 有不同的类型（Type-I, Type-II, Type-III, Type-IV）。scipy.fft.dct 和 scipy.fft.idct 允许使用 type 参数指定类型。默认的 dct 是 Type-II，其对应的逆变换 idct 是 Type-III。这些是压缩等应用中的常见选择。

scipy.fft 模块还提供了用于二维（和 N 维）FFT 和 DCT 的函数（例如，fft2、dctn），这对于图像处理和其他多维数据分析至关重要。