import osimport numpy as npimport matplotlib.pyplot as pltfrom matplotlib import cmfrom matplotlib.colors import Normalize# 0. 输出目录OUTDIR = "waterfall_demo_output_clean"os.makedirs(OUTDIR, exist_ok=True)# 1. 构造模拟的时间-频率-能量数据def gaussian(x, mu, sigma): return np.exp(-0.5 * ((x - mu) / sigma) ** 2)def moving_average_1d(x, w=5): if w <= 1: return x.copy() kernel = np.ones(w) / w return np.convolve(x, kernel, mode="same")def generate_tf_energy(nt=90, nf=260, seed=2026): rng = np.random.default_rng(seed) t = np.linspace(0, 40, nt) # time f = np.linspace(0, 120, nf) # frequency E = np.zeros((nt, nf), dtype=float) for i, ti in enumerate(t): # 成分1:低频起伏带 c1 = 18 + 10 * np.sin(2 * np.pi * ti / 18.0) a1 = 1.2 + 0.5 * np.cos(2 * np.pi * ti / 11.0) band1 = a1 * gaussian(f, c1, 4.5) # 成分2:中高频摆动带 c2 = 68 + 13 * np.cos(2 * np.pi * (ti - 5) / 15.0) a2 = 1.0 + 0.4 * np.sin(2 * np.pi * ti / 9.5) band2 = a2 * gaussian(f, c2, 6.5) # 成分3:斜向 chirp 型能量脊 c3 = 10 + 1.7 * ti a3 = 0.8 + 0.25 * np.sin(2 * np.pi * ti / 8.0) band3 = a3 * gaussian(f, c3, 3.0) # 局部强脉冲热点 burst = ( 1.8 * gaussian(ti, 10.0, 1.2) * gaussian(f, 88, 5.0) + 2.2 * gaussian(ti, 21.0, 1.5) * gaussian(f, 52, 4.0) + 1.6 * gaussian(ti, 31.0, 1.0) * gaussian(f, 97, 3.8) ) # 背景起伏 background = ( 0.15 + 0.10 * np.sin(0.15 * f + 0.55 * ti) + 0.06 * np.cos(0.10 * f - 0.35 * ti) ) # 噪声 noise = 0.05 * rng.normal(size=nf) E[i, :] = band1 + band2 + band3 + burst + background + noise E = np.clip(E, 0, None) # 沿时间方向平滑 for j in range(nf): E[:, j] = moving_average_1d(E[:, j], w=5) return t, f, E# 2. 基础版瀑布图def plot_basic_waterfall(t, f, E, out_png): fig = plt.figure(figsize=(12, 8), dpi=180) ax = fig.add_subplot(111, projection="3d") idx = np.arange(0, len(t), 2) norm = Normalize(vmin=E.min(), vmax=E.max()) cmap = cm.viridis for i in idx: color = cmap(norm(E[i].max())) ax.plot( f, np.full_like(f, t[i]), E[i], color=color, lw=1.4, alpha=0.95 ) ax.plot( f, np.full_like(f, t[i]), np.zeros_like(f), color=color, lw=0.6, alpha=0.12 ) ax.set_xlabel("Frequency", fontsize=12, labelpad=10) ax.set_ylabel("Time", fontsize=12, labelpad=10) ax.set_zlabel("Energy", fontsize=12, labelpad=8) ax.view_init(elev=28, azim=-62) ax.set_xlim(f.min(), f.max()) ax.set_ylim(t.min(), t.max()) ax.set_zlim(0, E.max() * 1.10) ax.xaxis.pane.set_alpha(0.06) ax.yaxis.pane.set_alpha(0.06) ax.zaxis.pane.set_alpha(0.02) ax.grid(True, alpha=0.25) mappable = cm.ScalarMappable(norm=norm, cmap=cmap) mappable.set_array([]) cbar = plt.colorbar(mappable, ax=ax, pad=0.08, shrink=0.78) cbar.set_label("Slice peak energy", fontsize=11) plt.tight_layout() plt.savefig(out_png, bbox_inches="tight")# 3. 重新整理后的进阶版瀑布图def plot_advanced_waterfall(t, f, E, out_png): fig = plt.figure(figsize=(12, 8), dpi=180) ax = fig.add_subplot(111, projection="3d") # 只抽取较少切片,避免太密 idx = np.arange(0, len(t), 4) norm = Normalize(vmin=E.min(), vmax=E.max()) cmap = cm.plasma peak_each_time = E.max(axis=1) top_idx = np.argsort(peak_each_time)[-3:] # 最强3个切片高亮 for i in idx: color = cmap(norm(E[i].max())) lw = 1.0 alpha = 0.45 if i in top_idx: lw = 2.2 alpha = 0.95 ax.plot( f, np.full_like(f, t[i]), E[i], color=color, lw=lw, alpha=alpha ) # 主能量脊线:每个时刻最大能量对应的频率 ridge_f = f[np.argmax(E, axis=1)] ridge_z = E.max(axis=1) ridge_f_smooth = moving_average_1d(ridge_f, w=7) ridge_z_smooth = moving_average_1d(ridge_z, w=7) # 主脊线本体 ax.plot( ridge_f_smooth, t, ridge_z_smooth + 0.03 * E.max(), color="cyan", lw=2.4, alpha=0.95, label="Dominant energy ridge" ) # 加少量散点,让这条线更有 3D 路径感 step = 6 ax.scatter( ridge_f_smooth[::step], t[::step], ridge_z_smooth[::step] + 0.03 * E.max(), s=18, color="cyan", edgecolors="black", linewidths=0.4, alpha=0.95 ) # 只标出一个全局峰值,不再写 Peak1 Peak2 g_idx = np.argmax(peak_each_time) g_f = f[np.argmax(E[g_idx])] g_z = E[g_idx].max() ax.scatter( [g_f], [t[g_idx]], [g_z], s=52, color="yellow", edgecolors="black", linewidths=0.8, zorder=10 ) ax.set_xlabel("Frequency", fontsize=12, labelpad=10) ax.set_ylabel("Time", fontsize=12, labelpad=10) ax.set_zlabel("Energy", fontsize=12, labelpad=8) ax.set_xlim(f.min(), f.max()) ax.set_ylim(t.min(), t.max()) ax.set_zlim(0, E.max() * 1.10) # 视角略调整 ax.view_init(elev=27, azim=-58) ax.xaxis.pane.set_alpha(0.05) ax.yaxis.pane.set_alpha(0.05) ax.zaxis.pane.set_alpha(0.02) ax.grid(True, alpha=0.18) mappable = cm.ScalarMappable(norm=norm, cmap=cmap) mappable.set_array([]) cbar = plt.colorbar(mappable, ax=ax, pad=0.08, shrink=0.78) cbar.set_label("Slice peak energy", fontsize=11) ax.legend(loc="upper right", frameon=True, fontsize=10) plt.tight_layout() plt.savefig(out_png, bbox_inches="tight")# 4. 主程序def main(): t, f, E = generate_tf_energy() basic_png = os.path.join(OUTDIR, "waterfall_basic.png") adv_png = os.path.join(OUTDIR, "waterfall_advanced_clean.png") plot_basic_waterfall(t, f, E, basic_png) plot_advanced_waterfall(t, f, E, adv_png)if __name__ == "__main__": main()