信号处理基础：信号的时域和频域分析_（28）.信号处理实验与实践
发布时间
阅读量：
阅读量
信号处理实验与实践
在上一节内容中，我们阐述了信号时域与频域分析的基础理论及其应用方法。本节内容将结合实验与实践环节,深入阐述上述概念的本质内涵；并通过具体的案例演示展示理论知识的实际运用。我们将利用Python语言及其科学计算库（如NumPy、SciPy和Matplotlib）来进行相关信号处理操作的演示教学。
1. 信号的时域分析实验
1.1 信号的生成与可视化
在此章节中, 我们将深入探讨如何生成多种多样的典型信号（如正弦波, 方波和噪声信号）, 并通过Matplotlib库实现这些信号的可视化展示.
1.1.1 生成正弦波信号
import numpy as np
import matplotlib.pyplot as plt
# 生成正弦波信号
def generate_sine_wave(freq, amplitude, duration, sampling_rate):
"""
生成正弦波信号
参数:
freq (float): 信号的频率 (Hz)
amplitude (float): 信号的振幅
duration (float): 信号的持续时间 (秒)
sampling_rate (int): 采样率 (Hz)
返回:
t (np.array): 时间向量
y (np.array): 信号向量
"""
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
y = amplitude * np.sin(2 * np.pi * freq * t)
return t, y
# 参数设置
freq = 10 # 频率 (Hz)
amplitude = 1 # 振幅
duration = 1 # 持续时间 (秒)
sampling_rate = 1000 # 采样率 (Hz)
# 生成信号
t, y = generate_sine_wave(freq, amplitude, duration, sampling_rate)
# 可视化信号
plt.figure(figsize=(10, 4))
plt.plot(t, y)
plt.title('正弦波信号')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')
plt.grid(True)
plt.show()1.1.2 生成方波信号
from scipy import signal
# 生成方波信号
def generate_square_wave(freq, amplitude, duration, sampling_rate):
"""
生成方波信号
参数:
freq (float): 信号的频率 (Hz)
amplitude (float): 信号的振幅
duration (float): 信号的持续时间 (秒)
sampling_rate (int): 采样率 (Hz)
返回:
t (np.array): 时间向量
y (np.array): 信号向量
"""
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
y = amplitude * signal.square(2 * np.pi * freq * t)
return t, y
# 参数设置
freq = 5 # 频率 (Hz)
amplitude = 1 # 振幅
duration = 1 # 持续时间 (秒)
sampling_rate = 1000 # 采样率 (Hz)
# 生成信号
t, y = generate_square_wave(freq, amplitude, duration, sampling_rate)
# 可视化信号
plt.figure(figsize=(10, 4))
plt.plot(t, y)
plt.title('方波信号')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')
plt.grid(True)
plt.show()1.1.3 生成噪声信号
# 生成噪声信号
def generate_noise(duration, sampling_rate):
"""
生成噪声信号
参数:
duration (float): 信号的持续时间 (秒)
sampling_rate (int): 采样率 (Hz)
返回:
t (np.array): 时间向量
y (np.array): 噪声信号向量
"""
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
y = np.random.normal(0, 1, t.size)
return t, y
# 参数设置
duration = 1 # 持续时间 (秒)
sampling_rate = 1000 # 采样率 (Hz)
# 生成信号
t, y = generate_noise(duration, sampling_rate)
# 可视化信号
plt.figure(figsize=(10, 4))
plt.plot(t, y)
plt.title('噪声信号')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')
plt.grid(True)
plt.show()1.2 信号的时域特性分析
本节内容中, 我们将深入研究信号的时间域特征, 其中涉及的主要统计指标包括均值、方差以及自相关与互相关函数。
1.2.1 计算信号的均值和方差
# 计算信号的均值和方差
def calculate_mean_variance(y):
"""
计算信号的均值和方差
参数:
y (np.array): 信号向量
返回:
mean (float): 信号的均值
variance (float): 信号的方差
"""
mean = np.mean(y)
variance = np.var(y)
return mean, variance
# 生成正弦波信号
t, y = generate_sine_wave(10, 1, 1, 1000)
# 计算均值和方差
mean, variance = calculate_mean_variance(y)
print(f'均值: {mean}')
print(f'方差: {variance}')1.2.2 计算信号的自相关函数
# 计算信号的自相关函数
def calculate_autocorrelation(y, lags):
"""
计算信号的自相关函数
参数:
y (np.array): 信号向量
lags (int): 自相关函数的滞后数
返回:
acf (np.array): 自相关函数值
"""
acf = np.correlate(y, y, mode='full')
acf = acf[acorn.size - lags: acorn.size + lags + 1] / np.var(y)
return acf
# 生成正弦波信号
t, y = generate_sine_wave(10, 1, 1, 1000)
# 计算自相关函数
lags = 100
acf = calculate_autocorrelation(y, lags)
# 可视化自相关函数
plt.figure(figsize=(10, 4))
plt.plot(acf)
plt.title('自相关函数')
plt.xlabel('滞后数')
plt.ylabel('自相关值')
plt.grid(True)
plt.show()1.2.3 计算信号的互相关函数
# 计算信号的互相关函数
def calculate_crosscorrelation(y1, y2):
"""
计算两个信号的互相关函数
参数:
y1 (np.array): 信号1向量
y2 (np.array): 信号2向量
返回:
ccf (np.array): 互相关函数值
"""
ccf = np.correlate(y1, y2, mode='full') / (np.sqrt(np.var(y1) * np.var(y2)))
return ccf
# 生成两个正弦波信号
t1, y1 = generate_sine_wave(10, 1, 1, 1000)
t2, y2 = generate_sine_wave(10, 1, 1, 1000)
# 计算互相关函数
ccf = calculate_crosscorrelation(y1, y2)
# 可视化互相关函数
plt.figure(figsize=(10, 4))
plt.plot(ccf)
plt.title('互相关函数')
plt.xlabel('滞后数')
plt.ylabel('互相关值')
plt.grid(True)
plt.show()2. 信号的频域分析实验
2.1 信号的傅里叶变换
在本节中, 我们将深入探讨利用傅里叶变换将时域信号转换为频域信号的同时进行频谱分析。
2.1.1 计算信号的傅里叶变换
import scipy.fft as fft
# 计算信号的傅里叶变换
def calculate_fft(y, sampling_rate):
"""
计算信号的傅里叶变换
参数:
y (np.array): 信号向量
sampling_rate (int): 采样率 (Hz)
返回:
freqs (np.array): 频率向量
fft_vals (np.array): 傅里叶变换值
"""
N = y.size
y_fft = fft.fft(y)
freqs = fft.fftfreq(N, 1/sampling_rate)
fft_vals = np.abs(y_fft)
return freqs, fft_vals
# 生成正弦波信号
t, y = generate_sine_wave(10, 1, 1, 1000)
# 计算傅里叶变换
freqs, fft_vals = calculate_fft(y, 1000)
# 可视化频谱
plt.figure(figsize=(10, 4))
plt.plot(freqs, fft_vals)
plt.title('频谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.grid(True)
plt.show()2.1.2 计算信号的功率谱密度
from scipy.signal import welch
# 计算信号的功率谱密度
def calculate_psd(y, sampling_rate, nperseg):
"""
计算信号的功率谱密度
参数:
y (np.array): 信号向量
sampling_rate (int): 采样率 (Hz)
nperseg (int): 每段的长度
返回:
freqs (np.array): 频率向量
psd (np.array): 功率谱密度值
"""
freqs, psd = welch(y, fs=sampling_rate, nperseg=nperseg)
return freqs, psd
# 生成正弦波信号
t, y = generate_sine_wave(10, 1, 1, 1000)
# 计算功率谱密度
freqs, psd = calculate_psd(y, 1000, 100)
# 可视化功率谱密度
plt.figure(figsize=(10, 4))
plt.semilogy(freqs, psd)
plt.title('功率谱密度')
plt.xlabel('频率 (Hz)')
plt.ylabel('功率密度')
plt.grid(True)
plt.show()2.2 信号的滤波器设计与应用
在此部分中, 我们将深入探讨如何掌握各类滤波器的设计与应用方法, 并以处理信号为例进行详细讲解。
2.2.1 设计和应用低通滤波器
from scipy.signal import butter, filtfilt
# 设计低通滤波器
def design_lowpass_filter(cutoff, fs, order=5):
"""
设计低通滤波器
参数:
cutoff (float): 截止频率 (Hz)
fs (int): 采样率 (Hz)
order (int): 滤波器阶数
返回:
b (np.array): 滤波器分子系数
a (np.array): 滤波器分母系数
"""
nyquist = 0.5 * fs
normal_cutoff = cutoff / nyquist
b, a = butter(order, normal_cutoff, btype='low', analog=False)
return b, a
# 应用低通滤波器
def apply_lowpass_filter(y, b, a):
"""
应用低通滤波器
参数:
y (np.array): 信号向量
b (np.array): 滤波器分子系数
a (np.array): 滤波器分母系数
返回:
y_filtered (np.array): 滤波后的信号向量
"""
y_filtered = filtfilt(b, a, y)
return y_filtered
# 生成正弦波信号
t, y = generate_sine_wave(10, 1, 1, 1000)
# 添加高频噪声
y_noise = y + 0.5 * np.random.normal(0, 1, y.size)
# 设计低通滤波器
cutoff = 20 # 截止频率 (Hz)
fs = 1000 # 采样率 (Hz)
b, a = design_lowpass_filter(cutoff, fs)
# 应用低通滤波器
y_filtered = apply_lowpass_filter(y_noise, b, a)
# 可视化原始信号和滤波后的信号
plt.figure(figsize=(10, 4))
plt.plot(t, y_noise, label='带噪声的信号')
plt.plot(t, y_filtered, label='低通滤波后的信号', color='red')
plt.title('低通滤波器效果')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')
plt.legend()
plt.grid(True)
plt.show()2.2.2 设计和应用高通滤波器
# 设计高通滤波器
def design_highpass_filter(cutoff, fs, order=5):
"""
设计高通滤波器
参数:
cutoff (float): 截止频率 (Hz)
fs (int): 采样率 (Hz)
order (int): 滤波器阶数
返回:
b (np.array): 滤波器分子系数
a (np.array): 滤波器分母系数
"""
nyquist = 0.5 * fs
normal_cutoff = cutoff / nyquist
b, a = butter(order, normal_cutoff, btype='high', analog=False)
return b, a
# 应用高通滤波器
def apply_highpass_filter(y, b, a):
"""
应用高通滤波器
参数:
y (np.array): 信号向量
b (np.array): 滤波器分子系数
a (np.array): 滤波器分母系数
返回:
y_filtered (np.array): 滤波后的信号向量
"""
y_filtered = filtfilt(b, a, y)
return y_filtered
# 生成低频噪声信号
t, y = generate_sine_wave(1, 1, 1, 1000)
# 添加高频信号
y_noise = y + 0.5 * generate_sine_wave(10, 1, 1, 1000)[1]
# 设计高通滤波器
cutoff = 5 # 截止频率 (Hz)
fs = 1000 # 采样率 (Hz)
b, a = design_highpass_filter(cutoff, fs)
# 应用高通滤波器
y_filtered = apply_highpass_filter(y_noise, b, a)
# 可视化原始信号和滤波后的信号
plt.figure(figsize=(10, 4))
plt.plot(t, y_noise, label='带低频噪声的信号')
plt.plot(t, y_filtered, label='高通滤波后的信号', color='red')
plt.title('高通滤波器效果')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')
plt.legend()
plt.grid(True)
plt.show()2.2.3 设计和应用带通滤波器
# 设计带通滤波器
def design_bandpass_filter(lowcut, highcut, fs, order=5):
"""
设计带通滤波器
参数:
lowcut (float): 低截止频率 (Hz)
highcut (float): 高截止频率 (Hz)
fs (int): 采样率 (Hz)
order (int): 滤波器阶数
返回:
b (np.array): 滤波器分子系数
a (np.array): 滤波器分母系数
"""
nyquist = 0.5 * fs
low = lowcut / nyquist
high = highcut / nyquist
b, a = butter(order, [low, high], btype='band', analog=False)
return b, a
# 应用带通滤波器
def apply_bandpass_filter(y, b, a):
"""
应用带通滤波器
参数:
y (np.array): 信号向量
b (np.array): 滤波器分子系数
a (np.array): 滤波器分母系数
返回:
y_filtered (np.array): 滤波后的信号向量
"""
y_filtered = filtfilt(b, a, y)
return y_filtered
# 生成多个频率的信号
t = np.linspace(0, 1, 1000, endpoint=False)
y = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.2 * np.sin(2 * np.pi * 100 * t)
# 设计带通滤波器
lowcut = 20 # 低截止频率 (Hz)
highcut = 80 # 高截止频率 (Hz)
fs = 1000 # 采样率 (Hz)
b, a = design_bandpass_filter(lowcut, highcut, fs)
# 应用带通滤波器
y_filtered = apply_bandpass_filter(y, b, a)
# 可视化原始信号和滤波后的信号
plt.figure(figsize=(10, 4))
plt.plot(t, y, label='原始信号')
plt.plot(t, y_filtered, label='带通滤波后的信号', color='red')
plt.title('带通滤波器效果')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')
plt.legend()
plt.grid(True)
plt.show()3. 信号的时频分析实验
3.1 短时傅里叶变换 (STFT)
本节将介绍短时傅里叶变换（Short-Time Fourier Transform, STFT）这一工具如何用于分析信号的时频特性。该方法通过将信号划分为有限时间窗并对其进行Fourier Transform处理，在不同时间段上揭示信号的频率组成信息。
3.1.1 计算信号的短时傅里叶变换
from scipy.signal import stft
# 计算信号的短时傅里叶变换
def calculate_stft(y, fs, window='hann', nperseg=256, noverlap=None):
"""
计算信号的短时傅里叶变换
参数:
y (np.array): 信号向量
fs (int): 采样率 (Hz)
window (str): 窗函数类型
nperseg (int): 每段的长度
noverlap (int): 段之间的重叠长度
返回:
f (np.array): 频率向量
t (np.array): 时间向量
Zxx (np.array): 短时傅里叶变换值 (复数)
"""
f, t, Zxx = stft(y, fs=fs, window=window, nperseg=nperseg, noverlap=noverlap)
return f, t, Zxx
# 生成复合信号
t = np.linspace(0, 1, 1000, endpoint=False)
y = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.2 * np.sin(2 * np.pi * 100 * t)
# 计算短时傅里叶变换
fs = 1000 # 采样率 (Hz)
f, t, Zxx = calculate_stft(y, fs)
# 可视化短时傅里叶变换
plt.figure(figsize=(10, 6))
plt.pcolormesh(t, f, np.abs(Zxx), shading='gouraud')
plt.title('短时傅里叶变换 (STFT)')
plt.xlabel('时间 (秒)')
plt.ylabel('频率 (Hz)')
plt.colorbar(label='幅度')
plt.grid(True)
plt.show()3.1.2 短时傅里叶变换的逆变换 (ISTFT)
逆短时傅里叶变换（ISTFT）可实现地将短时傅里叶变换的结果准确地还原为原始时域信号，在信号处理领域具有重要意义。它不仅能够有效检验STFT方法的正确性，还可以用于信号重构过程中的关键验证步骤。
from scipy.signal import istft
# 计算信号的逆短时傅里叶变换
def calculate_istft(Zxx, fs, window='hann', nperseg=256, noverlap=None):
"""
计算信号的逆短时傅里叶变换
参数:
Zxx (np.array): 短时傅里叶变换值 (复数)
fs (int): 采样率 (Hz)
window (str): 窗函数类型
nperseg (int): 每段的长度
noverlap (int): 段之间的重叠长度
返回:
t (np.array): 时间向量
y (np.array): 逆变换后的信号向量
"""
t, y = istft(Zxx, fs=fs, window=window, nperseg=nperseg, noverlap=noverlap)
return t, y
# 生成复合信号
t = np.linspace(0, 1, 1000, endpoint=False)
y = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.2 * np.sin(2 * np.pi * 100 * t)
# 计算短时傅里叶变换
fs = 1000 # 采样率 (Hz)
f, t, Zxx = calculate_stft(y, fs)
# 计算逆短时傅里叶变换
t_istft, y_istft = calculate_istft(Zxx, fs)
# 可视化原始信号和逆变换后的信号
plt.figure(figsize=(10, 4))
plt.plot(t, y, label='原始信号')
plt.plot(t_istft, y_istft, label='逆变换后的信号', color='red')
plt.title('逆短时傅里叶变换 (ISTFT) 效果')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')
plt.legend()
plt.grid(True)
plt.show()3.2 小波变换 (WT)
小波变换（Wavelet Transform, WT）主要用于对信号进行多尺度分析，并能有效揭示其时间与局部频率特征。相较于短时傅里叶变换（STFT）而言，在处理非平稳信号方面展现出显著的优势。
3.2.1 计算信号的小波变换
import pywt
# 计算信号的小波变换
def calculate_wavelet_transform(y, wavelet='db4', level=5):
"""
计算信号的小波变换
参数:
y (np.array): 信号向量
wavelet (str): 小波基函数
level (int): 分解层数
返回:
coeffs (list): 小波变换系数列表
"""
coeffs = pywt.wavedec(y, wavelet, level=level)
return coeffs
# 生成复合信号
t = np.linspace(0, 1, 1000, endpoint=False)
y = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.2 * np.sin(2 * np.pi * 100 * t)
# 计算小波变换
wavelet = 'db4' # 小波基函数
level = 5 # 分解层数
coeffs = calculate_wavelet_transform(y, wavelet, level)
# 可视化小波变换系数
plt.figure(figsize=(10, 6))
for i, coeff in enumerate(coeffs):
plt.plot(coeff, label=f'Level {i}')
plt.title('小波变换系数')
plt.xlabel('样本点')
plt.ylabel('系数值')
plt.legend()
plt.grid(True)
plt.show()3.2.2 小波变换的逆变换 (IWT)
逆小波变换（Inverse Wavelet Transform, IWT）能够将小波变换的结果还原为原始时间域信号，并在数学上表示为y(t) = IWT\{Y(w)\}或者y = IWT(Y)。这种技术在信号处理领域同样具有重要意义，并且能够用于验证小波正交性以及实现信号的重构过程。
# 计算信号的逆小波变换
def calculate_inverse_wavelet_transform(coeffs, wavelet='db4'):
"""
计算信号的逆小波变换
参数:
coeffs (list): 小波变换系数列表
wavelet (str): 小波基函数
返回:
y (np.array): 逆变换后的信号向量
"""
y = pywt.waverec(coeffs, wavelet)
return y
# 生成复合信号
t = np.linspace(0, 1, 1000, endpoint=False)
y = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.2 * np.sin(2 * np.pi * 100 * t)
# 计算小波变换
wavelet = 'db4' # 小波基函数
level = 5 # 分解层数
coeffs = calculate_wavelet_transform(y, wavelet, level)
# 计算逆小波变换
y_ist = calculate_inverse_wavelet_transform(coeffs, wavelet)
# 可视化原始信号和逆变换后的信号
plt.figure(figsize=(10, 4))
plt.plot(t, y, label='原始信号')
plt.plot(t, y_ist, label='逆变换后的信号', color='red')
plt.title('逆小波变换 (IWT) 效果')
plt.xlabel('时间 (秒)')
plt.ylabel('振幅')
plt.legend()
plt.grid(True)
plt.show()3.3 信号的时频联合分析
本次将深入研究如何采用时频联合分析方法（例如梅尔频率倒谱系数（MFCC）和时频图）来研究信号的时频特性。
3.3.1 计算梅尔频率倒谱系数 (MFCC)
梅尔频率倒谱系数（Mel-Frequency Cepstral Coefficients, MFCC）是音频信号处理中广泛应用的一种重要特征提取工具，在语音识别和音乐分类等领域发挥着关键作用。
import librosa
import librosa.display
# 生成音频信号
y, sr = librosa.load(librosa.ex('trumpet'), duration=1)
# 计算MFCC
mfcc = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)
# 可视化MFCC
plt.figure(figsize=(10, 4))
librosa.display.specshow(mfcc, x_axis='time', sr=sr)
plt.title('梅尔频率倒谱系数 (MFCC)')
plt.xlabel('时间 (秒)')
plt.ylabel('MFCC 系数')
plt.colorbar(format='%+2.0f dB')
plt.grid(True)
plt.show()3.3.2 生成时频图 (Spectrogram)
时频图（Spectrogram）是一种显示信号在时间-频率域中分布状态的方式，并能直观地呈现其频率成分在其发展过程中的变化情况。
# 生成时频图
def generate_spectrogram(y, fs, nperseg=256, noverlap=None):
"""
生成时频图
参数:
y (np.array): 信号向量
fs (int): 采样率 (Hz)
nperseg (int): 每段的长度
noverlap (int): 段之间的重叠长度
返回:
f (np.array): 频率向量
t (np.array): 时间向量
Sxx (np.array): 时频图值 (复数)
"""
f, t, Sxx = stft(y, fs=fs, window='hann', nperseg=nperseg, noverlap=noverlap)
return f, t, np.abs(Sxx)
# 生成复合信号
t = np.linspace(0, 1, 1000, endpoint=False)
y = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.2 * np.sin(2 * np.pi * 100 * t)
# 生成时频图
fs = 1000 # 采样率 (Hz)
f, t, Sxx = generate_spectrogram(y, fs)
# 可视化时频图
plt.figure(figsize=(10, 6))
plt.pcolormesh(t, f, Sxx, shading='gouraud')
plt.title('时频图 (Spectrogram)')
plt.xlabel('时间 (秒)')
plt.ylabel('频率 (Hz)')
plt.colorbar(label='幅度')
plt.grid(True)
plt.show()4. 总结
基于一系列实验与实践的基础上
全部评论 (0)
还没有任何评论哟~