来自 FFT 的频率与幅度的二维热图

Question

如何创建 2D 热图（或类似的），其中 x 轴为频率，y 轴为幅度？

import numpy as np
import matplotlib.pyplot as plt

fs = 1000
t = np.arange(0, 0.1, 1/fs)
N = len(t)
f_bin = fs / N
f = np.arange(0, fs, f_bin)
X = [np.fft.fft(np.sin(2 * np.pi * 100 * t)), np.fft.fft(np.sin(2 * np.pi * 200 * t)), np.fft.fft(np.sin(2 * np.pi * 300 * t))]
M = np.absolute(X)

fig1, (ax1) = plt.subplots(nrows=1, ncols=1)
ax1.plot(f[:len(f)//2], M[0,:len(M[0][:])//2], "r+", label="100 Hz")
ax1.plot(f[:len(f)//2], M[1,:len(M[1][:])//2], "g*", label="200 Hz")
ax1.plot(f[:len(f)//2], M[2,:len(M[2][:])//2], "b.", label="300 Hz")
ax1.set_xlabel("Frequency (Hz)")
ax1.set_ylabel("Magnitude")
ax1.set_title(f"DFT")
ax1.grid(True, which="both")
ax1.legend()
plt.show()

总的来说，我寻找一种好方法来在单张图片中可视化多个（例如，1000 个，而不仅仅是本例中的 3 个）FFT 的输出（例如，使用

matplotlib.imshow()

？、

numpy.histogram2d()

？）。

Answer 1

这听起来像是与聚类、异常检测或异常值检测相关的任务。

sklearn

有用于此类任务的工具。

一种方法是使用

BayesianMixtureModel

获得幅度与频率的密度图。如果新线落在较暗（低密度）区域，您可以将其标记为异常。这些数字的代码如下。

原始数据：

用于异常检测的密度估计数据：

import numpy as np
import matplotlib.pyplot as plt

fs = 10_000
t = np.arange(0, 0.5, 1/fs)
N = len(t)
f_bin = fs / N
f = np.arange(0, fs, f_bin)
v = [np.sinc(2 * np.pi * f * t)**2*(np.random.rand(len(t))>0.5) for f in range(100, 110)]
X = [np.fft.fft(time) for time in v]
X_arr = np.stack(X, axis=0)
M = np.absolute(X_arr)

fig1, ax = plt.subplots()
ax.plot(f[:N//2], M[:, :N//2].T)
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Magnitude")
ax.set_title("DFT")
plt.show()

#Fit clusterer to get density plot
from sklearn.mixture import BayesianGaussianMixture
f_rpt = np.tile(f[:N//2].reshape(-1, 1), [len(v), 1])
M_flat = M[:, :N//2].reshape(-1, 1)

f_and_M = np.concatenate([f_rpt, M_flat], axis=1)
bgm = BayesianGaussianMixture(n_components=20, random_state=0)
bgm.fit(f_and_M)

#Plotting code adapted from:
#   github.com/ageron/handson-ml3/blob/main/09_unsupervised_learning.ipynb
from matplotlib.colors import LogNorm
X = f_and_M
clusterer = bgm
resolution = 1000
mins = X.min(axis=0) - 0.1
maxs = X.max(axis=0) + 0.1
xx, yy = np.meshgrid(np.linspace(mins[0], maxs[0], resolution),
                     np.linspace(mins[1], maxs[1], resolution))
Z = -clusterer.score_samples(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)

plt.contourf(xx, yy, Z,
            norm=LogNorm(vmin=1.0, vmax=30.0),
            levels=np.logspace(0, 2, 12), cmap='plasma_r')
plt.contour(xx, yy, Z,
            norm=LogNorm(vmin=1.0, vmax=30.0),
            levels=np.logspace(0, 2, 12),
            linewidths=1, colors='k')

plt.plot(X[:, 0], X[:, 1], 'w.', markersize=2)
plt.xlabel("frequency")
plt.ylabel("magnitude")
plt.title('Density plot (white lines are the original data)')