时域信号转PSD谱实战python代码

文章手把手教你把时域信号转成 PSD 谱手把手教你把时域信号转成 PSD 谱(补充)对应的python代码如下，供大家参考，欢迎指正：

import numpy as np

import matplotlib.pyplot as plt

from scipy import signal

from numpy.fft import fft, fftfreq

import matplotlib.ticker as ticker

from matplotlib.ticker import FuncFormatter

# ==========================================

# 0. 全局设置 (隔离策略)

# ==========================================

# 1. 字体顺序：数字/西文优先用 Arial，中文用 SimHei

plt.rcParams['font.sans-serif'] = ['Arial', 'SimHei', 'Microsoft YaHei']

# 2. 负号处理：使用 ASCII 短横线，保证图1-4线性坐标正常显示

plt.rcParams['axes.unicode_minus'] = False 

# 3. 关闭全局 MathText，只在图6局部开启

plt.rcParams['axes.formatter.use_mathtext'] = False

plt.close('all') 

# ==========================================

# 1. 数据生成

# ==========================================

fs = 10240  # 采样率

T = 2.0     # 时长

N = int(T * fs)

t = np.linspace(0, T, N, endpoint=False)

np.random.seed(42)

noise = np.random.normal(scale=0.5, size=N)

sig1 = 2.0 * np.sin(2 * np.pi * 400 * t) 

sig2 = 0.8 * np.sin(2 * np.pi * 1500 * t)

x = noise + sig1 + sig2

nperseg = 2048   

noverlap = nperseg // 2

# ==========================================

# Step 1: 图 1 - 原始时域信号

# ==========================================

print("正在生成图 1...")

plt.figure(figsize=(12, 4))

plt.plot(t, x, color='grey', alpha=0.8, linewidth=0.3)

plt.title(f'图1：原始时域加速度信号 (时长 T={T}s)')

plt.xlabel('时间 [s]')

plt.ylabel('幅度 [g]')

plt.grid(True, linestyle='--')

plt.xlim(0, T)

plt.show()

# ==========================================

# Step 2: 图 2 - 单个数据片段

# ==========================================

segment_data = x[0:nperseg]

segment_time = t[0:nperseg]

print("正在生成图 2...")

plt.figure(figsize=(12, 4))

plt.plot(segment_time, segment_data, color='#1f77b4')

plt.title(f'图2：被切分出的单个数据段 (点数 N={nperseg})')

plt.xlabel('时间 [s]')

plt.ylabel('幅度 [g]')

plt.grid(True, linestyle='--')

plt.xlim(segment_time[0], segment_time[-1])

plt.show()

# ==========================================

# Step 3: 图 3 - 汉宁窗

# ==========================================

window = signal.windows.hann(nperseg)

print("正在生成图 3...")

plt.figure(figsize=(12, 3))

plt.plot(np.arange(nperseg),window,color='#2ca02c', linewidth=2)

plt.fill_between(np.arange(nperseg),window,color='#2ca02c',alpha=0.2)

plt.title('图3：准备施加的“汉宁窗 (Hanning Window)”')

plt.grid(True, linestyle='--')

plt.show()

# ==========================================

# Step 4: 图 4 - 加窗效果

# ==========================================

windowed_segment = segment_data * window

print("正在生成图 4...")

plt.figure(figsize=(12, 4))

plt.plot(segment_time,segment_data,color='grey',alpha=0.3,label='原始')

plt.plot(segment_time,windowed_segment,color='#d62728',linewidth=1.5, label='加窗后')

plt.title('图4：加窗后的数据段 (两端归零，防止泄露)')

plt.legend()

plt.grid(True, linestyle='--')

plt.xlim(segment_time[0], segment_time[-1])

plt.show()

# ==========================================

# Step 5: 图 5 - 单片段 FFT (全频段)

# ==========================================

fft_vals = fft(windowed_segment)

raw_spectrum = np.abs(fft_vals[:nperseg//2])**2

freqs_temp = fftfreq(nperseg, 1/fs)[:nperseg//2]

print("正在生成图 5...")

plt.figure(figsize=(12, 5))

plt.semilogy(freqs_temp,raw_spectrum,color='#9467bd',alpha=0.6, linewidth=0.8)

plt.title('图5：仅对一个片段做 FFT (充满了随机噪声/毛刺)')

plt.xlabel('频率 [Hz]')

plt.grid(True, which='both', linestyle='--')

plt.xlim(0, fs/2) # 这里是 5120Hz

plt.show()

# ==========================================

# Step 6: 图 6 - 最终成果 (改为全频段一致)

# ==========================================

# 计算 Welch PSD

f_psd,Pxx_den=signal.welch(x,fs,window='hann',nperseg=nperseg, noverlap=noverlap, scaling='density')

# 计算 RMS

rms_time = np.sqrt(np.mean(x**2))

rms_freq = np.sqrt(np.trapezoid(Pxx_den, f_psd))

error_pct = abs(rms_time - rms_freq) / rms_time * 100

print("正在生成图 6...")

fig, ax = plt.subplots(figsize=(12, 5))

# 绘图

ax.semilogy(f_psd, Pxx_den, color='black', linewidth=2, label='Welch 平均 PSD')

ax.semilogy(freqs_temp,raw_spectrum/nperseg/fs*4,color='#9467bd', alpha=0.15, label='单片段背景')

# --- 【关键修改】横坐标范围改为 fs/2 ---

# 这样就和图5保持一致，显示完整的 0-5120Hz

ax.set_xlim(0, fs/2) 

ax.set_ylim(1e-7, 1) 

# --- 纵坐标刻度格式化 (局部强制 LaTeX) ---

def format_log_ticks(x, pos):

    exponent = int(np.log10(x))

    return f'$10^{{{exponent}}}$'

ax.yaxis.set_major_formatter(FuncFormatter(format_log_ticks))

ax.yaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20))

# --- 非均匀密集网格 ---

loc_minor = ticker.LogLocator(base=10.0, subs=np.arange(2, 10), numticks=100)

ax.yaxis.set_minor_locator(loc_minor)

ax.minorticks_on()

ax.grid(True,which='major',linestyle='--', linewidth=0.8, alpha=0.8)

ax.grid(True, which='minor', axis='y', linestyle='--', linewidth=0.5, alpha=0.4)

ax.grid(True, which='minor', axis='x', linestyle=':', linewidth=0.5, alpha=0.3)

# 标签

ax.set_title('图6：最终成果 —— Welch 平滑 PSD 与 RMS 校验')

ax.set_xlabel('频率 [Hz]')

ax.set_ylabel('PSD [$g^2$/Hz]') 

# --- 布局：左右排列 ---

text_str = (f"Time RMS: {rms_time:.4f} g\n"

            f"Freq RMS: {rms_freq:.4f} g\n"

            f"Error:    {error_pct:.2f}%")

props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

# 文本框 (最右)

ax.text(0.98, 0.95, text_str, transform=ax.transAxes, fontsize=12,

        verticalalignment='top', horizontalalignment='right', bbox=props)

# 图例 (文本框左侧)

ax.legend(loc='upper right', bbox_to_anchor=(0.72, 0.97))

plt.tight_layout()

plt.show()