使用Matplotlib在Python中绘制相位谱|极客教程

使用Matplotlib在Python中绘制相位谱

参考：Plot the phase spectrum in Python using Matplotlib

相位谱是信号处理和频域分析中的重要工具，它可以帮助我们了解信号在不同频率下的相位特性。在Python中，我们可以使用强大的Matplotlib库来绘制相位谱，从而直观地展示信号的相位信息。本文将详细介绍如何使用Matplotlib在Python中绘制相位谱，包括基本概念、数据准备、绘图技巧以及高级应用等方面。

1. 相位谱的基本概念

相位谱是描述信号在不同频率下相位变化的图形表示。它与幅度谱一起构成了信号的完整频域表示。相位谱的横轴表示频率，纵轴表示相位角度（通常以弧度或度为单位）。

在进行相位谱分析时，我们通常会使用傅里叶变换将时域信号转换为频域表示。对于离散信号，我们使用离散傅里叶变换（DFT）或快速傅里叶变换（FFT）来实现这一转换。

以下是一个简单的示例，展示如何使用NumPy和Matplotlib计算并绘制一个简单信号的相位谱：

import numpy as np
import matplotlib.pyplot as plt

# 生成信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)

# 计算FFT
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(t), t[1] - t[0])

# 计算相位谱
phase_spectrum = np.angle(fft_result)

# 绘制相位谱
plt.figure(figsize=(10, 6))
plt.plot(frequencies[:len(frequencies)//2], phase_spectrum[:len(frequencies)//2])
plt.title("Phase Spectrum - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase (radians)")
plt.grid(True)
plt.show()

Output:

使用Matplotlib在Python中绘制相位谱

在这个示例中，我们首先生成了一个包含两个不同频率正弦波的信号。然后，我们使用NumPy的FFT函数计算信号的傅里叶变换，并使用np.angle()函数计算相位谱。最后，我们使用Matplotlib绘制了相位谱图。

2. 数据准备和预处理

在绘制相位谱之前，我们通常需要对原始数据进行一些预处理。这可能包括信号去噪、窗函数应用、零填充等操作。以下是一些常见的预处理步骤：

2.1 信号去噪

信号去噪可以帮助我们获得更清晰的相位谱。一种简单的去噪方法是使用移动平均滤波器：

import numpy as np
import matplotlib.pyplot as plt

def moving_average(data, window_size):
    return np.convolve(data, np.ones(window_size)/window_size, mode='valid')

# 生成带噪声的信号
t = np.linspace(0, 1, 1000, endpoint=False)
clean_signal = np.sin(2 * np.pi * 10 * t)
noise = np.random.normal(0, 0.1, len(t))
noisy_signal = clean_signal + noise

# 应用移动平均滤波器
filtered_signal = moving_average(noisy_signal, 5)

# 绘制原始信号和滤波后的信号
plt.figure(figsize=(12, 6))
plt.plot(t, noisy_signal, label='Noisy Signal')
plt.plot(t[:len(filtered_signal)], filtered_signal, label='Filtered Signal')
plt.title("Signal Denoising - how2matplotlib.com")
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.legend()
plt.grid(True)
plt.show()

Output:

使用Matplotlib在Python中绘制相位谱

这个示例展示了如何使用移动平均滤波器对带噪声的信号进行简单的去噪处理。在实际应用中，可能需要根据具体情况选择更复杂的滤波方法。

2.2 窗函数应用

窗函数可以帮助减少频谱泄漏，提高相位谱的准确性。以下是一个使用汉宁窗的示例：

import numpy as np
import matplotlib.pyplot as plt

# 生成信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t)

# 应用汉宁窗
window = np.hanning(len(signal))
windowed_signal = signal * window

# 计算FFT
fft_original = np.fft.fft(signal)
fft_windowed = np.fft.fft(windowed_signal)

# 计算相位谱
phase_original = np.angle(fft_original)
phase_windowed = np.angle(fft_windowed)

# 绘制相位谱对比
frequencies = np.fft.fftfreq(len(t), t[1] - t[0])
plt.figure(figsize=(12, 6))
plt.plot(frequencies[:len(frequencies)//2], phase_original[:len(frequencies)//2], label='Original')
plt.plot(frequencies[:len(frequencies)//2], phase_windowed[:len(frequencies)//2], label='Windowed')
plt.title("Phase Spectrum Comparison - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase (radians)")
plt.legend()
plt.grid(True)
plt.show()

Output:

使用Matplotlib在Python中绘制相位谱

这个示例展示了如何应用汉宁窗函数，并比较了原始信号和加窗信号的相位谱。

2.3 零填充

零填充可以增加频率分辨率，使相位谱更加平滑。以下是一个零填充的示例：

import numpy as np
import matplotlib.pyplot as plt

# 生成信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t)

# 零填充
zero_padded_signal = np.pad(signal, (0, len(signal)), 'constant')

# 计算FFT
fft_original = np.fft.fft(signal)
fft_zero_padded = np.fft.fft(zero_padded_signal)

# 计算相位谱
phase_original = np.angle(fft_original)
phase_zero_padded = np.angle(fft_zero_padded)

# 绘制相位谱对比
frequencies_original = np.fft.fftfreq(len(t), t[1] - t[0])
frequencies_zero_padded = np.fft.fftfreq(len(zero_padded_signal), t[1] - t[0])

plt.figure(figsize=(12, 6))
plt.plot(frequencies_original[:len(frequencies_original)//2], phase_original[:len(frequencies_original)//2], label='Original')
plt.plot(frequencies_zero_padded[:len(frequencies_zero_padded)//2], phase_zero_padded[:len(frequencies_zero_padded)//2], label='Zero Padded')
plt.title("Phase Spectrum Comparison with Zero Padding - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase (radians)")
plt.legend()
plt.grid(True)
plt.show()

Output:

使用Matplotlib在Python中绘制相位谱

这个示例展示了如何使用零填充来增加频率分辨率，并比较了原始信号和零填充信号的相位谱。

3. 相位谱绘制技巧

在使用Matplotlib绘制相位谱时，有一些技巧可以帮助我们创建更加清晰、信息丰富的图表。

3.1 相位展开

相位谱通常被限制在-π到π的范围内，这可能导致相位跳变。我们可以使用相位展开技术来消除这些跳变，使相位谱更加连续：

import numpy as np
import matplotlib.pyplot as plt

# 生成信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)

# 计算FFT
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(t), t[1] - t[0])

# 计算相位谱
phase_spectrum = np.angle(fft_result)

# 相位展开
unwrapped_phase = np.unwrap(phase_spectrum)

# 绘制原始相位谱和展开后的相位谱
plt.figure(figsize=(12, 6))
plt.plot(frequencies[:len(frequencies)//2], phase_spectrum[:len(frequencies)//2], label='Original')
plt.plot(frequencies[:len(frequencies)//2], unwrapped_phase[:len(frequencies)//2], label='Unwrapped')
plt.title("Phase Spectrum Comparison - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase (radians)")
plt.legend()
plt.grid(True)
plt.show()

Output:

使用Matplotlib在Python中绘制相位谱

这个示例展示了如何使用NumPy的unwrap()函数来展开相位谱，使其更加连续。

3.2 对数频率轴

在某些情况下，使用对数频率轴可以更好地展示低频和高频的相位特性：

import numpy as np
import matplotlib.pyplot as plt

# 生成信号
t = np.linspace(0, 1, 10000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 100 * t) + 0.2 * np.sin(2 * np.pi * 1000 * t)

# 计算FFT
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(t), t[1] - t[0])

# 计算相位谱
phase_spectrum = np.angle(fft_result)

# 绘制相位谱（对数频率轴）
plt.figure(figsize=(12, 6))
plt.semilogx(frequencies[:len(frequencies)//2], phase_spectrum[:len(frequencies)//2])
plt.title("Phase Spectrum with Log Frequency Axis - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase (radians)")
plt.grid(True)
plt.show()

Output:

使用Matplotlib在Python中绘制相位谱

这个示例使用了Matplotlib的semilogx()函数来创建一个横轴为对数刻度的相位谱图。

3.3 相位谱平滑

有时，相位谱可能会出现噪声或不规则波动。我们可以使用平滑技术来改善相位谱的可读性：

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter

# 生成信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + np.random.normal(0, 0.1, len(t))

# 计算FFT
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(t), t[1] - t[0])

# 计算相位谱
phase_spectrum = np.angle(fft_result)

# 应用Savitzky-Golay滤波器进行平滑
smoothed_phase = savgol_filter(phase_spectrum, window_length=51, polyorder=3)

# 绘制原始相位谱和平滑后的相位谱
plt.figure(figsize=(12, 6))
plt.plot(frequencies[:len(frequencies)//2], phase_spectrum[:len(frequencies)//2], label='Original', alpha=0.5)
plt.plot(frequencies[:len(frequencies)//2], smoothed_phase[:len(frequencies)//2], label='Smoothed')
plt.title("Phase Spectrum Smoothing - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase (radians)")
plt.legend()
plt.grid(True)
plt.show()

Output:

使用Matplotlib在Python中绘制相位谱

这个示例使用了SciPy的Savitzky-Golay滤波器来平滑相位谱。

4. 高级相位谱分析

除了基本的相位谱绘制，我们还可以进行一些高级的相位谱分析。

4.1 相位延迟和群延迟

相位延迟和群延迟是描述信号在不同频率下延迟特性的重要指标。我们可以从相位谱计算这些指标：

import numpy as np
import matplotlib.pyplot as plt

# 生成信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)

# 计算FFT
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(t), t[1] - t[0])

# 计算相位谱
phase_spectrum = np.unwrap(np.angle(fft_result))

# 计算相位延迟和群延迟
phase_delay = -phase_spectrum / (2 * np.pi * frequencies)
group_delay = -np.diff(phase_spectrum) / np.diff(2 * np.pi * frequencies)

# 绘制相位延迟和群延迟
plt.figure(figsize=(12, 10))

plt.subplot(3, 1, 1)
plt.plot(frequencies[:len(frequencies)//2], phase_spectrum[:len(frequencies)//2])
plt.title("PhaseSpectrum - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase (radians)")
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(frequencies[1:len(frequencies)//2], phase_delay[1:len(frequencies)//2])
plt.title("Phase Delay - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Delay (seconds)")
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(frequencies[1:len(frequencies)//2-1], group_delay[:len(frequencies)//2-1])
plt.title("Group Delay - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Delay (seconds)")
plt.grid(True)

plt.tight_layout()
plt.show()

这个示例展示了如何计算和绘制相位谱、相位延迟和群延迟。相位延迟表示信号在每个频率上的延迟时间，而群延迟表示信号包络的延迟时间。

4.2 相位一致性分析

相位一致性分析可以帮助我们识别信号中的周期性成分。以下是一个简单的相位一致性分析示例：

import numpy as np
import matplotlib.pyplot as plt

def phase_consistency(signal, window_size, overlap):
    num_windows = (len(signal) - window_size) // (window_size - overlap) + 1
    fft_size = window_size
    consistency = np.zeros(fft_size // 2 + 1)

    for i in range(num_windows):
        start = i * (window_size - overlap)
        end = start + window_size
        window = signal[start:end] * np.hanning(window_size)
        fft = np.fft.rfft(window)
        phase = np.angle(fft)
        consistency += np.exp(1j * phase)

    consistency = np.abs(consistency) / num_windows
    return consistency

# 生成信号
t = np.linspace(0, 10, 10000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + np.random.normal(0, 0.1, len(t))

# 计算相位一致性
window_size = 1000
overlap = 500
consistency = phase_consistency(signal, window_size, overlap)

# 绘制相位一致性
frequencies = np.fft.rfftfreq(window_size, t[1] - t[0])
plt.figure(figsize=(12, 6))
plt.plot(frequencies, consistency)
plt.title("Phase Consistency Analysis - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase Consistency")
plt.ylim(0, 1)
plt.grid(True)
plt.show()

这个示例展示了如何计算和绘制相位一致性。相位一致性值越接近1，表示该频率成分在信号中的相位越稳定。

4.3 时频相位分析

时频相位分析可以帮助我们了解信号相位随时间和频率的变化。短时傅里叶变换（STFT）是一种常用的时频分析方法：

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import stft

# 生成信号
t = np.linspace(0, 10, 10000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * (20 + t) * t)

# 计算STFT
f, t, Zxx = stft(signal, fs=1000, nperseg=1000, noverlap=900)

# 提取相位信息
phase = np.angle(Zxx)

# 绘制时频相位图
plt.figure(figsize=(12, 6))
plt.pcolormesh(t, f, phase, cmap='hsv', shading='gouraud')
plt.title("Time-Frequency Phase Analysis - how2matplotlib.com")
plt.xlabel("Time (s)")
plt.ylabel("Frequency (Hz)")
plt.colorbar(label="Phase (radians)")
plt.ylim(0, 50)  # 限制频率范围以便更好地观察
plt.show()

Output:

使用Matplotlib在Python中绘制相位谱

这个示例使用了SciPy的stft函数来计算短时傅里叶变换，并绘制了时频相位图。颜色表示不同时间和频率下的相位值。

5. 实际应用案例

相位谱分析在许多领域都有重要应用。以下是一些实际应用案例：

5.1 音频信号处理

在音频信号处理中，相位谱可以用于分析音频的谐波结构、识别音高和音色特征等。以下是一个简单的音频信号相位谱分析示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile

# 读取音频文件
sample_rate, audio_data = wavfile.read("example_audio.wav")  # 请替换为实际的音频文件路径

# 选择一个短时间段进行分析
start_time = 1  # 起始时间（秒）
duration = 0.1  # 持续时间（秒）
start_sample = int(start_time * sample_rate)
end_sample = start_sample + int(duration * sample_rate)
audio_segment = audio_data[start_sample:end_sample]

# 计算FFT
fft_result = np.fft.fft(audio_segment)
frequencies = np.fft.fftfreq(len(audio_segment), 1/sample_rate)

# 计算相位谱
phase_spectrum = np.angle(fft_result)

# 绘制相位谱
plt.figure(figsize=(12, 6))
plt.plot(frequencies[:len(frequencies)//2], phase_spectrum[:len(frequencies)//2])
plt.title("Audio Signal Phase Spectrum - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase (radians)")
plt.grid(True)
plt.show()

这个示例展示了如何读取音频文件，选择一个短时间段，并计算其相位谱。

5.2 图像处理

在图像处理中，相位谱可以用于边缘检测、图像配准等任务。以下是一个简单的图像相位谱分析示例：

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

# 读取图像
image = np.array(Image.open("example_image.jpg").convert("L"))  # 请替换为实际的图像文件路径

# 计算2D FFT
fft_result = np.fft.fft2(image)
fft_shifted = np.fft.fftshift(fft_result)

# 计算相位谱
phase_spectrum = np.angle(fft_shifted)

# 绘制原始图像和相位谱
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.imshow(image, cmap='gray')
plt.title("Original Image - how2matplotlib.com")
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(phase_spectrum, cmap='hsv')
plt.title("Phase Spectrum - how2matplotlib.com")
plt.axis('off')

plt.tight_layout()
plt.show()

这个示例展示了如何读取图像，计算其2D FFT，并绘制相位谱。

5.3 通信系统分析

在通信系统中，相位谱分析可以用于评估信道特性、检测相位噪声等。以下是一个简单的通信信号相位谱分析示例：

import numpy as np
import matplotlib.pyplot as plt

# 生成QPSK信号
num_symbols = 1000
symbols = np.random.choice([1+1j, 1-1j, -1+1j, -1-1j], num_symbols)
t = np.arange(num_symbols)
carrier_freq = 0.1
signal = np.real(symbols * np.exp(1j * 2 * np.pi * carrier_freq * t))

# 添加相位噪声
phase_noise = np.cumsum(np.random.normal(0, 0.01, len(signal)))
noisy_signal = signal * np.exp(1j * phase_noise)

# 计算FFT
fft_result = np.fft.fft(noisy_signal)
frequencies = np.fft.fftfreq(len(t), t[1] - t[0])

# 计算相位谱
phase_spectrum = np.angle(fft_result)

# 绘制原始信号和相位谱
plt.figure(figsize=(12, 10))

plt.subplot(3, 1, 1)
plt.plot(t, signal)
plt.title("Original QPSK Signal - how2matplotlib.com")
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(t, noisy_signal)
plt.title("QPSK Signal with Phase Noise - how2matplotlib.com")
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(frequencies[:len(frequencies)//2], phase_spectrum[:len(frequencies)//2])
plt.title("Phase Spectrum - how2matplotlib.com")
plt.xlabel("Frequency")
plt.ylabel("Phase (radians)")
plt.grid(True)

plt.tight_layout()
plt.show()