Python Spectrum模块介绍与使用

Spectrum是一个用于信号处理的Python模块，专门用于频谱分析和谱估计。它提供了多种经典和现代的谱估计方法。

安装

pip install spectrum

主要功能

1. 经典谱估计方法

import spectrum
import numpy as np
import matplotlib.pyplot as plt
from spectrum import data_cosine, arma2psd

# 生成测试信号
N = 1024
fs = 1000  # 采样频率
t = np.arange(N) / fs
f1, f2 = 50, 120
signal = np.sin(2*np.pi*f1*t) + 0.5*np.sin(2*np.pi*f2*t) + 0.1*np.random.randn(N)

# 周期图法 (Periodogram)
from spectrum import Periodogram
p = Periodogram(signal, sampling=fs)
p.plot()
plt.title('Periodogram')
plt.show()

# Welch方法
from spectrum import pwelch
pw = pwelch(signal, sampling=fs)
pw.plot()
plt.title('Welch Method')
plt.show()

2. 现代谱估计方法

# AR模型谱估计
from spectrum import aryule, arma2psd, marple_data

# 使用Yule-Walker方法估计AR参数
order = 30  # AR模型阶数
ar_params, variance, coeffs = aryule(signal, order)
psd = arma2psd(ar_params, NFFT=1024)
frequencies = np.linspace(0, fs/2, len(psd)//2)

plt.figure()
plt.plot(frequencies, 10*np.log10(psd[:len(psd)//2]))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (dB)')
plt.title('AR Model Spectrum (Yule-Walker)')
plt.grid(True)
plt.show()

3. 多 taper谱估计

from spectrum import pmtm

# 多锥形谱估计
[psd, weights, eigenvalues] = pmtm(signal, NW=4, NFFT=1024, sampling=fs)
frequencies = np.linspace(0, fs/2, len(psd)//2)

plt.figure()
plt.plot(frequencies, 10*np.log10(psd[:len(psd)//2]))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (dB)')
plt.title('Multitaper Spectrum')
plt.grid(True)
plt.show()

4. 参数化方法

# MUSIC算法
from spectrum import music

# 生成包含多个频率成分的信号
freqs = [50, 120, 200]
signal_music = sum([np.sin(2*np.pi*f*t) for f in freqs]) + 0.1*np.random.randn(N)

psd_music, weights = music(signal_music, 15, 4, NFFT=1024, sampling=fs)
frequencies = np.linspace(0, fs/2, len(psd_music)//2)

plt.figure()
plt.plot(frequencies, 10*np.log10(psd_music[:len(psd_music)//2]))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (dB)')
plt.title('MUSIC Algorithm')
plt.grid(True)
plt.show()

5. 时频分析

from spectrum import STFT

# 短时傅里叶变换
stft = STFT(signal, 256, 128, sampling=fs)
stft.run()
stft.plot()
plt.title('Short-Time Fourier Transform')
plt.show()

完整示例

import numpy as np
import matplotlib.pyplot as plt
from spectrum import *

# 设置参数
fs = 1000  # 采样频率
N = 2000   # 数据点数
t = np.arange(N) / fs

# 创建测试信号：两个正弦波加噪声
f1, f2 = 50, 150
signal = (np.sin(2*np.pi*f1*t) + 
          0.5*np.sin(2*np.pi*f2*t) + 
          0.2*np.random.randn(N))

# 比较不同的谱估计方法
methods = [
    ('Periodogram', Periodogram(signal, sampling=fs)),
    ('Welch', pwelch(signal, sampling=fs)),
]

# 创建子图
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 绘制原始信号
axes[0,0].plot(t[:500], signal[:500])
axes[0,0].set_xlabel('Time (s)')
axes[0,0].set_ylabel('Amplitude')
axes[0,0].set_title('Original Signal (first 500 points)')
axes[0,0].grid(True)

# 绘制经典谱估计方法
for i, (name, method) in enumerate(methods):
    ax = axes[0,1] if i == 0 else axes[1,0]
    method.plot(ax=ax)
    ax.set_title(f'{name} Spectrum')
    ax.grid(True)

# AR模型谱估计
order = 30
ar_params, variance, coeffs = aryule(signal, order)
psd_ar = arma2psd(ar_params, NFFT=1024)
freqs = np.linspace(0, fs/2, len(psd_ar)//2)

axes[1,1].plot(freqs, 10*np.log10(psd_ar[:len(psd_ar)//2]))
axes[1,1].set_xlabel('Frequency (Hz)')
axes[1,1].set_ylabel('Power (dB)')
axes[1,1].set_title('AR Model Spectrum')
axes[1,1].grid(True)

plt.tight_layout()
plt.show()

# 性能比较
print("谱估计方法比较完成")

主要类和函数

• Periodogram: 经典周期图法
• pwelch: Welch方法（改进的周期图法）
• aryule: Yule-Walker AR参数估计
• pmtm: 多锥形谱估计
• music: MUSIC算法
• STFT: 短时傅里叶变换
• arma2psd: 将ARMA参数转换为功率谱密度

优点

1. 算法丰富: 提供多种经典和现代谱估计方法
2. 易于使用: 简单的API设计
3. 可视化支持: 内置绘图功能
4. 文档完善: 提供详细的文档和示例

注意事项

1. 选择合适的模型阶数对参数化方法很重要
2. 不同方法适用于不同的信号特性
3. 对于非平稳信号，考虑使用时频分析方法

Spectrum模块是信号处理领域非常有用的工具，特别适合频谱分析和谱估计任务。