相关性解码
- 相关匹配/模板匹配 (Correlation/Template Matching) 解码方法特别适用于信号质量较差（有噪声、失真）或需要高可靠性的场景。其核心思想是：将接收到的波形片段与预先定义的“理想”波形模板进行对比，计算它们的相似度（相关值），并选择相似度最高的模板作为解码结果。它的优势在低信噪比、高可靠性要求的场景下尤为突出。随着计算能力的提升，这种经典的信号处理方法在现代智能系统（如自动驾驶、智能医疗设备）中依然扮演着不可或缺的角色：
- 数字调制解调：
- BPSK/QPSK 解调：接收端使用本地生成的载波和已知的符号波形作为模板，与接收到的信号做相关，以判决发送的是
0 还是 1 或哪个相位。 - 扩频通信 (如 CDMA)：利用伪随机码（PN码）作为模板进行相关解扩。只有使用相同码的接收机才能从宽带噪声中“拉出”有用信号，实现多址接入和抗干扰。
- 同步字检测：
- 在数据帧开头插入特定的“同步头”（如
0101... 或巴克码）。接收端用此同步头作为模板，在接收流中滑动匹配，一旦找到高相关值，即确认帧起始位置，完成帧同步。
- 雷达通信一体化：在通信信号中嵌入雷达探测功能，利用相关性检测目标回波。
| 优点 | 缺点 |
|---|
| 抗噪声能力强：即使信号有失真，只要整体形状匹配，仍能正确解码。 | 计算量大：每个位都需要 2 次点积运算（TEMPLATE_LEN 次乘加），对 CPU 要求高。 |
| 可靠性高：相比简单阈值法，误码率更低。 | 需要预知采样率和位率：模板依赖于固定的 SAMPLE_RATE 和 BIT_RATE。 |
| 易于实现：逻辑清晰，代码结构简单。 | 内存占用：需要存储模板数组和中间结果存储。 |
| 可扩展：可轻松添加同步头、状态字等更复杂的模板。 | 对时钟漂移敏感：如果实际位周期与模板不匹配，性能下降。 |
实现逻辑
- 模板生成（
0 和 1 的曼彻斯特编码波形） - 滑动窗口相关匹配
- 同步头检测（3 位同步头）
- 帧解析
- 奇校验验证
- 完整的解码流程演示
python实现
import numpy as np
import matplotlib.pyplot as plt
from typing import List, Tuple, Optional
from scipy import signal
# ========================
# 配置参数
# ========================
# 采样信息
SAMPLE_RATE = 10_000_000 # 10 MHz 采样率
BIT_RATE = 1_000_000 # 1 Mbps (MIL-STD-1553B)
HALF_BIT_SAMPLES = SAMPLE_RATE // BIT_RATE // 2 # 半位采样点数 5
BIT_SAMPLES = 2 * HALF_BIT_SAMPLES # 一位总采样点数 10
# 信号电平（模拟ADC值）
HIGH_LEVEL = 10000
LOW_LEVEL = -10000
# 帧结构
ALL_DATA_BITS = 17 # 16个数据位和一个校验位
SYNC_LENGTH = 15 # 同步头
# ========================
# 模板生成
# ========================
def generate_manchester_templates() -> Tuple[np.ndarray, np.ndarray]:
"""
生成 MIL-STD-1553B 曼彻斯特编码模板
- '1': 正脉冲 + 负脉冲 -> [HIGH, LOW]
- '0': 负脉冲 + 正脉冲 -> [LOW, HIGH]
"""
template_1 = np.concatenate([
np.full(HALF_BIT_SAMPLES, HIGH_LEVEL),
np.full(HALF_BIT_SAMPLES, LOW_LEVEL)
]) # []
template_0 = np.concatenate([
np.full(HALF_BIT_SAMPLES, LOW_LEVEL),
np.full(HALF_BIT_SAMPLES, HIGH_LEVEL)
])
return template_0, template_1
# 生成同步头模板
def generate_sync_template(template_1: np.ndarray, num_sync_bits: int = 3) -> np.ndarray:
"""生成同步头模板"""
return np.concatenate([np.zeros(SYNC_LENGTH, dtype=int), np.ones(SYNC_LENGTH, dtype=int)]) # 共30位，15个低位+15个高位。
# ========================
# 相关匹配查找同步头
# ========================
def find_sync_start_correlation(waveform: np.ndarray,
length: int,
sync_template: np.ndarray,
threshold: float = 75000.0) -> int:
"""
使用互相关查找同步头位置，仅当相关值超过阈值且为局部最大值时才认为是有效同步头。
参数:
waveform (np.ndarray): 输入的波形数据
length (int): 搜索范围长度（从 waveform 开头搜索）
sync_template (np.ndarray): 同步头模板
threshold (float): 相关值阈值
返回:
int: 第一个有效同步头的起始位置（索引），未找到返回 -1
"""
# 截取搜索范围
search_wave = waveform[:length]
# 计算互相关，mode='valid' 表示只有完全重叠的部分
corr = np.correlate(search_wave, sync_template, mode='valid')
# 如果相关结果为空，直接返回-1
if len(corr) == 0:
return -1
# 存储局部最大值的位置
peak_positions = []
# 遍历相关值，寻找局部最大值
for i in range(1, len(corr) - 1): # 忽略边界
if corr[i] > threshold: # 先判断是否超过阈值
# 判断是否为局部最大值（峰值）
if corr[i] > corr[i-1] and corr[i] > corr[i+1]:
peak_positions.append(i)
# 如果没有找到任何峰值，返回 -1
if not peak_positions:
return -1
# 返回第一个检测到的峰值位置（对应 waveform 的起始索引）
# 注意：corr[i] 对应 waveform[i] 到 waveform[i + len(template)-1]
return peak_positions[0]
# ========================
# 相关匹配解码单个位
# ========================
def match_bit(signal_segment: np.ndarray, template_0: np.ndarray, template_1: np.ndarray) -> int:
"""
使用点积相关匹配判断是 0 还是 1
返回: 0, 1, 或 -1（无效）
"""
if len(signal_segment) != len(template_0):
raise ValueError("Signal segment length does not match template length")
corr_0 = np.dot(signal_segment, template_0)
corr_1 = np.dot(signal_segment, template_1)
if corr_0 > corr_1:
return 0
elif corr_1 > corr_0:
return 1
else:
return -1 # 平局，视为错误
def match_bit_r_corr(signal_segment: np.ndarray, template_0: np.ndarray, template_1: np.ndarray) -> int:
"""
使用点积相关匹配判断是 0 还是 1
返回: 0, 1, 或 -1（无效）
"""
if len(signal_segment) != len(template_0):
raise ValueError("Signal segment length does not match template length")
corr_0 = np.dot(signal_segment, template_0) # -1*(signal_segment[0]+signal_segment[1]+signal_segment[2]+signal_segment[3]+signal_segment[4])+(signal_segment[5]+signal_segment[6]+signal_segment[7]+signal_segment[8]+signal_segment[9])
corr_1 = np.dot(signal_segment, template_1) # (signal_segment[0]+signal_segment[1]+signal_segment[2]+signal_segment[3]+signal_segment[4])-1*(signal_segment[5]+signal_segment[6]+signal_segment[7]+signal_segment[8]+signal_segment[9])
if corr_0 > corr_1:
return 0, corr_0
elif corr_1 > corr_0:
return 1, corr_1
else:
return -1,0 # 平局，视为错误
# ========================
# 滑动窗口法寻找极值位
# ========================
def find_peaks_sliding_window(data, window_size=11, threshold=0.3):
"""
使用滑动窗口检测局部极值
参数:
data: 输入时间序列
window_size: 滑动窗口大小（必须为奇数）
threshold: 相对幅度阈值
返回:
positive_peaks: 正峰值位置
negative_peaks: 负峰值位置
"""
if window_size % 2 == 0:
window_size += 1 # 确保窗口大小为奇数
half_window = window_size // 2
positive_peaks = []
negative_peaks = []
# 计算全局统计量用于阈值
data_range = np.max(data) - np.min(data)
abs_threshold = threshold * data_range
for i in range(half_window, len(data) - half_window):
window = data[i - half_window:i + half_window + 1]
center_value = data[i]
# 检查是否是正峰值
if center_value == np.max(window) and center_value > np.min(window) + abs_threshold:
# 确保是严格的局部最大值
is_peak = True
for j in range(1, half_window + 1):
if data[i - j] >= center_value or data[i + j] >= center_value:
is_peak = False
break
if is_peak:
positive_peaks.append(i)
# 检查是否是负峰值
elif center_value == np.min(window) and center_value < np.max(window) - abs_threshold:
# 确保是严格的局部最小值
is_valley = True
for j in range(1, half_window + 1):
if data[i - j] <= center_value or data[i + j] <= center_value:
is_valley = False
break
if is_valley:
negative_peaks.append(i)
return np.array(positive_peaks), np.array(negative_peaks)
# ========================
# 滑动窗口解码位流
# ========================
def decode_bits_from_waveform(waveform: np.ndarray,
template_0: np.ndarray,
template_1: np.ndarray) -> List[int]:
"""从波形中解码出位流"""
bits = []
i = 0
j = 0
while (i <= len(waveform) - BIT_SAMPLES) & (j<17):
segment = waveform[i:i + BIT_SAMPLES]
bit,corr = match_bit_r_corr(segment, template_0, template_1) # match_bit时平局，视为错误，返回-1
if bit != -1:
bits.append(bit)
else:
bits.append(-1) # 标记错误位
i += BIT_SAMPLES # 移动一个位
j += 1 # 解码的数据个数加一
return bits
# ========================
# 数据校验 如果符合偶数校验则返回True，否则返回False
# ========================
def even_parity_check(decoded_bits):
"""
对给定的二进制列表进行偶数奇偶校验。
参数:
decoded_bits (list): 二进制位列表（0或1），最后一位被认为是奇偶校验位
返回:
bool: 如果符合偶数奇偶校验则返回True，否则返回False
"""
# 计算除了校验位外的所有1的数目
count_ones = sum(decoded_bits[:-1])
# 判断包括校验位在内的所有1的数目是否为偶数
return (count_ones + decoded_bits[-1]) % 2 == 0
# ========================
# 绘制数据
# ========================
def plot_data(data):
plt.figure(figsize=(10, 5)) # 设置图表大小
plt.plot(data, marker='o') # 绘制数据点，添加标记使数据点更明显
plt.title('Data Plot') # 添加标题
plt.xlabel('Index') # X轴标签
plt.ylabel('Value') # Y轴标签
plt.grid(True) # 显示网格线
# 显示图表
plt.show()
print(data)
# 绘制同步头的位置
def plot_data_mark(data,sync_positions):
"""
绘制数据，并在 sync_positions 指定的位置添加垂直线标记（如同步头位置）
参数:
data (list or np.ndarray): 要绘制的数据
sync_positions (list of int): 需要标记的位置索引列表，例如 [150, 450]
"""
plt.figure(figsize=(12, 6)) # 设置图表大小
plt.plot(data, marker='o', markersize=3, linewidth=1, label='Signal') # 绘制数据
plt.axvline(x=sync_positions, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.axvline(x=sync_positions+30, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.title('Data Plot with Sync Markers') # 添加标题
plt.xlabel('Sample Index') # X轴标签
plt.ylabel('Amplitude') # Y轴标签
plt.grid(True, alpha=0.6) # 显示网格线
# 图例：只显示一次 "Sync" 标签
handles, labels = plt.gca().get_legend_handles_labels()
if handles:
plt.legend()
# 显示图表
plt.tight_layout()
plt.show()
# 输出数据（可选）
print("Data:", data)
print("Sync positions marked:", sync_positions)
# 绘制多个同步头的位置
def plot_data_mark_mut(data,sync_positions):
"""
绘制数据，并在 sync_positions 指定的位置添加垂直线标记（如同步头位置）
参数:
data (list or np.ndarray): 要绘制的数据
sync_positions (list of int): 需要标记的位置索引列表，例如 [150, 450]
"""
plt.figure(figsize=(12, 6)) # 设置图表大小
plt.plot(data, marker='o', markersize=3, linewidth=1, label='Signal') # 绘制数据
for pos in sync_positions:
plt.axvline(x=pos, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.axvline(x=pos+30, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.title('Data Plot with Sync Markers') # 添加标题
plt.xlabel('Sample Index') # X轴标签
plt.ylabel('Amplitude') # Y轴标签
plt.grid(True, alpha=0.6) # 显示网格线
# 图例：只显示一次 "Sync" 标签
handles, labels = plt.gca().get_legend_handles_labels()
if handles:
plt.legend()
# 显示图表
plt.tight_layout()
plt.show()
# 输出数据（可选）
print("Data:", data)
print("Sync positions marked:", sync_positions)
def plot_data_mark_mut_dot(data, sync_positions):
"""
绘制数据，并在 sync_positions 指定的位置添加垂直线标记（如同步头位置）
在第二个红线后以10为距离添加标记，直到下一个第一条红线为止，标记显示对应index
"""
plt.figure(figsize=(12, 6)) # 设置图表大小
plt.plot(data, marker='o', markersize=3, linewidth=1, label='Signal') # 绘制数据
# 先绘制所有红线
red_lines = []
for pos in sync_positions:
# 第一条红线位置
first_line = pos
# 第二条红线位置
second_line = pos + 30
plt.axvline(x=first_line, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.axvline(x=second_line, color='red', linestyle='--', linewidth=2, alpha=0.8)
red_lines.append((first_line, second_line))
# 为每个区间添加浮动文本标记
for i, (first_line, second_line) in enumerate(red_lines):
# 确定当前区间的结束位置（下一个第一条红线）
if i < len(red_lines) - 1:
end_pos = red_lines[i+1][0] # 下一个第一条红线位置
else:
end_pos = len(data) - 1 # 最后一个区间以数据长度为终点
# 从第二条红线后10个位置开始
start_pos = second_line + 10
# 确保起始位置不超过结束位置
if start_pos > end_pos:
continue
# 以10为间隔添加浮动标记（显示index）
current_pos = start_pos
while current_pos <= end_pos:
# 获取当前位置的数据值（用于定位文本）
data_value = data[current_pos]
# 添加浮动文本标记（显示当前index）
# 位置：在点的右下方（可通过xytext调整偏移）
# 仅第一个标记添加图例标签
label = 'Index Mark' if (i == 0 and current_pos == start_pos) else ""
if(current_pos in score_dict):
plt.text(
current_pos, data_value, # 文本锚点位置（对应数据点）
f'{score_dict[current_pos]}', # 显示的文本内容（index值）
fontsize=8, # 文本大小
color='blue', # 文本颜色
ha='left', va='bottom', # 水平/垂直对齐方式
bbox=dict(facecolor='white', alpha=0.3, boxstyle='round,pad=0.3'), # 白色背景框
label=label
)
# print("----",current_pos)
current_pos += 10
plt.title('Data Plot with Sync Markers') # 添加标题
plt.xlabel('Sample Index') # X轴标签
plt.ylabel('Amplitude') # Y轴标签
plt.grid(True, alpha=0.6) # 显示网格线
# 处理图例，确保只显示唯一标签
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys())
# 显示图表
plt.tight_layout()
plt.show()
# 输出数据（可选）
print("Data:", data)
print("Sync positions marked:", sync_positions)
# ========================
# 数据滤波
# ========================
def lowpass_filter_scipy(data, cutoff_freq, sampling_rate, filter_type='butter', order=4):
"""
使用SciPy对NumPy数组进行低通滤波
参数:
data: numpy数组，输入信号
cutoff_freq: float, 截止频率 (Hz)
sampling_rate: float, 采样频率 (Hz)
filter_type: str, 滤波器类型 ('butter', 'cheby1', 'bessel')
order: int, 滤波器阶数
返回:
filtered_data: 滤波后的numpy数组
"""
# 计算归一化截止频率 (Nyquist频率的倍数)
nyquist_freq = 0.5 * sampling_rate
normal_cutoff = cutoff_freq / nyquist_freq
if filter_type == 'butter':
b, a = signal.butter(order, normal_cutoff, btype='low', analog=False)
elif filter_type == 'cheby1':
b, a = signal.cheby1(order, 0.5, normal_cutoff, btype='low', analog=False)
elif filter_type == 'bessel':
b, a = signal.bessel(order, normal_cutoff, btype='low', analog=False)
else:
raise ValueError("不支持的滤波器类型")
# 应用滤波器
filtered_data = signal.filtfilt(b, a, data)
return filtered_data
def moving_average_filter(data, window_size):
"""
使用移动平均进行低通滤波
参数:
data: numpy数组
window_size: int, 窗口大小
返回:
filtered_data: 滤波后的数组
"""
if window_size <= 1:
return data.copy()
# 使用卷积实现移动平均
window = np.ones(window_size) / window_size
filtered_data = np.convolve(data, window, mode='same')
return filtered_data
res_dict = {}
score_dict = {}
if __name__ == "__main__":
sync_head_num = 0 # 查找到同步头的个数
paritiy_num = 0
paritiy_list = []
not_par_num = 0
not_par_list = []
data = np.loadtxt('***.dat', skiprows=1) # 加载数据文件，跳过第一行
# 预处理数据
data = -1*data
print(data)
startindex = 0
data = data[startindex:startindex+49000] # 100000
# 0.预处理
# data = moving_average_filter(data,2)
data = lowpass_filter_scipy(data,40000,100000)
# 1. 生成模板
template_0, template_1 = generate_manchester_templates()
sync_template = generate_sync_template(template_1, 3)
# 2. 波形数据
test_wave = data
# plot_data(test_wave)
data_len = data.size
index = 0
sync_list = [] # 标记同步头位置的列表
while(index<(data_len-1000)): # 留1000个缓冲空间
# 查找同步头
sync_positions = -1
window_size = 300
window_overlap = 100
while(sync_positions ==-1):
index = index + window_size - window_overlap # 窗口重叠
segment = test_wave[index:-1]
if(segment.size != 0):
sync_positions = find_sync_start_correlation(segment, window_size, sync_template, threshold= 55000.0)
else:
print("已遍历所有数据")
break
if(segment.size != 0):
sync_head_num = sync_head_num +1
print("第一个同步头的位置",index+sync_positions)
# plot_data_mark(test_wave[0:-1],index+sync_positions) # python 不支持函数重载
sync_list.append(index+sync_positions)
# 3. 解码位流
data = test_wave[index+sync_positions+30:index+sync_positions+30+170]
print("+++",index+sync_positions+30)
decoded_bits = decode_bits_from_waveform(data, template_0, template_1) # 解完17个就返回
res_dict[index+sync_positions] = decoded_bits
value = index+sync_positions # +30+10+0*10
# score_dict[value] = r_list
# 4. 对list数据decoded_bits进行奇偶校验
is_valid = even_parity_check(decoded_bits)
print(f"奇偶校验{'通过' if is_valid else '未通过'}")
if is_valid:
paritiy_num = paritiy_num + 1
paritiy_list.append(index+sync_positions)
else:
not_par_num = not_par_num + 1
not_par_list.append(index+sync_positions)
print("-------------统计信息------------")
print("找到同步头",sync_head_num)
print("奇偶校验未通过个数",paritiy_num,paritiy_list)
print("奇偶校验通过个数",not_par_num,not_par_list)
print("---------------------------------")
print("-------------绘制同步头位置------------")
plot_data_mark_mut_dot(test_wave,sync_list)
print("---------------------------------")
互相关解码
import numpy as np
import matplotlib.pyplot as plt
from typing import List, Tuple, Optional
from scipy import signal
# ========================
# 配置参数
# ========================
# 采样信息
SAMPLE_RATE = 10_000_000 # 10 MHz 采样率
BIT_RATE = 1_000_000 # 1 Mbps (MIL-STD-1553B)
HALF_BIT_SAMPLES = SAMPLE_RATE // BIT_RATE // 2 # 半位采样点数 5
BIT_SAMPLES = 2 * HALF_BIT_SAMPLES # 一位总采样点数 10
# 信号电平（模拟ADC值）
HIGH_LEVEL = 10000
LOW_LEVEL = -10000
# 帧结构
ALL_DATA_BITS = 17 # 16个数据位和一个校验位
SYNC_LENGTH = 15 # 同步头
# ========================
# 模板生成
# ========================
def generate_manchester_templates() -> Tuple[np.ndarray, np.ndarray]:
"""
生成 MIL-STD-1553B 曼彻斯特编码模板
- '1': 正脉冲 + 负脉冲 -> [HIGH, LOW]
- '0': 负脉冲 + 正脉冲 -> [LOW, HIGH]
"""
template_1 = np.concatenate([
np.full(HALF_BIT_SAMPLES, HIGH_LEVEL),
np.full(HALF_BIT_SAMPLES, LOW_LEVEL)
]) # []
template_0 = np.concatenate([
np.full(HALF_BIT_SAMPLES, LOW_LEVEL),
np.full(HALF_BIT_SAMPLES, HIGH_LEVEL)
])
return template_0, template_1
# 生成同步头模板
def generate_sync_template(template_1: np.ndarray, num_sync_bits: int = 3) -> np.ndarray:
"""生成同步头模板"""
return np.concatenate([np.zeros(SYNC_LENGTH, dtype=int), np.ones(SYNC_LENGTH, dtype=int)]) # 共30位，15个低位+15个高位。
# ========================
# 相关匹配查找同步头
# ========================
def find_sync_start_correlation(waveform: np.ndarray,
length: int,
sync_template: np.ndarray,
threshold: float = 75000.0) -> int:
"""
使用互相关查找同步头位置，仅当相关值超过阈值且为局部最大值时才认为是有效同步头。
参数:
waveform (np.ndarray): 输入的波形数据
length (int): 搜索范围长度（从 waveform 开头搜索）
sync_template (np.ndarray): 同步头模板
threshold (float): 相关值阈值
返回:
int: 第一个有效同步头的起始位置（索引），未找到返回 -1
"""
# 截取搜索范围
search_wave = waveform[:length]
# 计算互相关，mode='valid' 表示只有完全重叠的部分
corr = np.correlate(search_wave, sync_template, mode='valid')
# 如果相关结果为空，直接返回-1
if len(corr) == 0:
return -1
# 存储局部最大值的位置
peak_positions = []
# 遍历相关值，寻找局部最大值
for i in range(1, len(corr) - 1): # 忽略边界
if corr[i] > threshold: # 先判断是否超过阈值
# 判断是否为局部最大值（峰值）
if corr[i] > corr[i-1] and corr[i] > corr[i+1]:
peak_positions.append(i)
# 如果没有找到任何峰值，返回 -1
if not peak_positions:
return -1
# 返回第一个检测到的峰值位置（对应 waveform 的起始索引）
# 注意：corr[i] 对应 waveform[i] 到 waveform[i + len(template)-1]
return peak_positions[0]
# ========================
# 相关匹配解码单个位
# ========================
def match_bit(signal_segment: np.ndarray, template_0: np.ndarray, template_1: np.ndarray) -> int:
"""
使用点积相关匹配判断是 0 还是 1
返回: 0, 1, 或 -1（无效）
"""
if len(signal_segment) != len(template_0):
raise ValueError("Signal segment length does not match template length")
corr_0 = np.dot(signal_segment, template_0)
corr_1 = np.dot(signal_segment, template_1)
if corr_0 > corr_1:
return 0
elif corr_1 > corr_0:
return 1
else:
return -1 # 平局，视为错误
def match_bit_r_corr(signal_segment: np.ndarray, template_0: np.ndarray, template_1: np.ndarray) -> int:
"""
使用点积相关匹配判断是 0 还是 1
返回: 0, 1, 或 -1（无效）
"""
if len(signal_segment) != len(template_0):
raise ValueError("Signal segment length does not match template length")
corr_0 = np.dot(signal_segment, template_0) # -1*(signal_segment[0]+signal_segment[1]+signal_segment[2]+signal_segment[3]+signal_segment[4])+(signal_segment[5]+signal_segment[6]+signal_segment[7]+signal_segment[8]+signal_segment[9])
corr_1 = np.dot(signal_segment, template_1) # (signal_segment[0]+signal_segment[1]+signal_segment[2]+signal_segment[3]+signal_segment[4])-1*(signal_segment[5]+signal_segment[6]+signal_segment[7]+signal_segment[8]+signal_segment[9])
if corr_0 > corr_1:
return 0, corr_0
elif corr_1 > corr_0:
return 1, corr_1
else:
return -1,0 # 平局，视为错误
# ========================
# 滑动窗口法寻找极值位
# ========================
def find_peaks_sliding_window(data, window_size=11, threshold=0.3):
"""
使用滑动窗口检测局部极值
参数:
data: 输入时间序列
window_size: 滑动窗口大小（必须为奇数）
threshold: 相对幅度阈值
返回:
positive_peaks: 正峰值位置
negative_peaks: 负峰值位置
"""
if window_size % 2 == 0:
window_size += 1 # 确保窗口大小为奇数
half_window = window_size // 2
positive_peaks = []
negative_peaks = []
# 计算全局统计量用于阈值
data_range = np.max(data) - np.min(data)
abs_threshold = threshold * data_range
for i in range(half_window, len(data) - half_window):
window = data[i - half_window:i + half_window + 1]
center_value = data[i]
# 检查是否是正峰值
if center_value == np.max(window) and center_value > np.min(window) + abs_threshold:
# 确保是严格的局部最大值
is_peak = True
for j in range(1, half_window + 1):
if data[i - j] >= center_value or data[i + j] >= center_value:
is_peak = False
break
if is_peak:
positive_peaks.append(i)
# 检查是否是负峰值
elif center_value == np.min(window) and center_value < np.max(window) - abs_threshold:
# 确保是严格的局部最小值
is_valley = True
for j in range(1, half_window + 1):
if data[i - j] <= center_value or data[i + j] <= center_value:
is_valley = False
break
if is_valley:
negative_peaks.append(i)
return np.array(positive_peaks), np.array(negative_peaks)
def find_alternating_peaks(time_series):
"""
查找时间序列中交替出现的正负峰值，要求正负峰值间数据单调
参数:
time_series: 一维数组或列表，代表输入的时间序列
返回:
peaks: 列表，包含检测到的峰值信息，每个元素为元组(位置索引, 峰值类型, 峰值大小)
其中峰值类型为'positive'(正峰)或'negative'(负峰)
"""
if len(time_series) < 3:
return [] # 太短的序列无法形成峰值
# 转换为numpy数组便于处理
ts = np.asarray(time_series)
peaks = []
zero_peak_type = 'positive'
if ts[0]<0:
zero_peak_type = 'negative'
peaks.append((0, zero_peak_type, ts[0]))
n = len(ts)
i = 1 # 从第二个元素开始检查
next_expected_type = 'negative' if zero_peak_type == 'positive' else 'positive'
while i < n - 1:
# 检查是否符合预期类型的峰值
is_positive_peak = (ts[i] > ts[i-1] and ts[i] > ts[i+1])
is_negative_peak = (ts[i] < ts[i-1] and ts[i] < ts[i+1])
if (i>=3) & (i<=158):
is_positive_peak = (ts[i] > ts[i-1] and ts[i] > ts[i+1]) & (ts[i] > ts[i-2] and ts[i] > ts[i+2])
is_negative_peak = (ts[i] < ts[i-1] and ts[i] < ts[i+1]) & (ts[i] < ts[i-2] and ts[i] < ts[i+2])
if (next_expected_type == 'positive' and is_positive_peak) or \
(next_expected_type == 'negative' and is_negative_peak):
peaks.append((i, next_expected_type, ts[i]))
next_expected_type = 'negative' if next_expected_type == 'positive' else 'positive'
i += 1
return peaks # 感觉还是要用滑动窗口方法
# ========================
# 滑动窗口解码位流
# ========================
def decode_bits_from_waveform(waveform: np.ndarray,
template_0: np.ndarray,
template_1: np.ndarray) -> List[int]:
"""从波形中解码出位流"""
bits = []
i = 0
j = 0
while (i <= len(waveform) - BIT_SAMPLES) & (j<17):
segment = waveform[i:i + BIT_SAMPLES]
bit,corr = match_bit_r_corr(segment, template_0, template_1) # match_bit时平局，视为错误，返回-1
if bit != -1:
bits.append(bit)
else:
bits.append(-1) # 标记错误位
i += BIT_SAMPLES # 移动一个位
j += 1 # 解码的数据个数加一
return bits
def decode_bits_from_waveform_r_peaks(waveform: np.ndarray,
template_0: np.ndarray,
template_1: np.ndarray):
"""从波形中解码出位流"""
bits = []
i = 0
corr = np.correlate(waveform, template_0, mode='valid')
res = find_alternating_peaks(corr) # 峰值数据
if res[0][1] == 'positive':
bits.append(0)
else:
bits.append(1)
flag = False
for i in range(1,len(res)):
if(flag):
if res[i][1] == 'positive':
bits.append(0)
else:
bits.append(1)
flag = not flag
if ((res[i][0] - res[i-1][0]) > 6+1):
flag = False
if res[i][1] == 'positive':
bits.append(0)
else:
bits.append(1)
else:
pass
return bits[0:17],res
# ========================
# 数据校验 如果符合偶数校验则返回True，否则返回False
# ========================
def even_parity_check(decoded_bits):
"""
对给定的二进制列表进行偶数奇偶校验。
参数:
decoded_bits (list): 二进制位列表（0或1），最后一位被认为是奇偶校验位
返回:
bool: 如果符合偶数奇偶校验则返回True，否则返回False
"""
# 计算除了校验位外的所有1的数目
count_ones = sum(decoded_bits[:-1])
# 判断包括校验位在内的所有1的数目是否为偶数
return (count_ones + decoded_bits[-1]) % 2 == 0
# ========================
# 绘制数据
# ========================
def plot_data(data):
plt.figure(figsize=(10, 5)) # 设置图表大小
plt.plot(data, marker='o') # 绘制数据点，添加标记使数据点更明显
plt.title('Data Plot') # 添加标题
plt.xlabel('Index') # X轴标签
plt.ylabel('Value') # Y轴标签
plt.grid(True) # 显示网格线
# 显示图表
plt.show()
print(data)
# 绘制同步头的位置
def plot_data_mark(data,sync_positions):
"""
绘制数据，并在 sync_positions 指定的位置添加垂直线标记（如同步头位置）
参数:
data (list or np.ndarray): 要绘制的数据
sync_positions (list of int): 需要标记的位置索引列表，例如 [150, 450]
"""
plt.figure(figsize=(12, 6)) # 设置图表大小
plt.plot(data, marker='o', markersize=3, linewidth=1, label='Signal') # 绘制数据
plt.axvline(x=sync_positions, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.axvline(x=sync_positions+30, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.title('Data Plot with Sync Markers') # 添加标题
plt.xlabel('Sample Index') # X轴标签
plt.ylabel('Amplitude') # Y轴标签
plt.grid(True, alpha=0.6) # 显示网格线
# 图例：只显示一次 "Sync" 标签
handles, labels = plt.gca().get_legend_handles_labels()
if handles:
plt.legend()
# 显示图表
plt.tight_layout()
plt.show()
# 输出数据（可选）
print("Data:", data)
print("Sync positions marked:", sync_positions)
# 绘制多个同步头的位置
def plot_data_mark_mut(data,sync_positions):
"""
绘制数据，并在 sync_positions 指定的位置添加垂直线标记（如同步头位置）
参数:
data (list or np.ndarray): 要绘制的数据
sync_positions (list of int): 需要标记的位置索引列表，例如 [150, 450]
"""
plt.figure(figsize=(12, 6)) # 设置图表大小
plt.plot(data, marker='o', markersize=3, linewidth=1, label='Signal') # 绘制数据
for pos in sync_positions:
plt.axvline(x=pos, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.axvline(x=pos+30, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.title('Data Plot with Sync Markers') # 添加标题
plt.xlabel('Sample Index') # X轴标签
plt.ylabel('Amplitude') # Y轴标签
plt.grid(True, alpha=0.6) # 显示网格线
# 图例：只显示一次 "Sync" 标签
handles, labels = plt.gca().get_legend_handles_labels()
if handles:
plt.legend()
# 显示图表
plt.tight_layout()
plt.show()
# 输出数据（可选）
print("Data:", data)
print("Sync positions marked:", sync_positions)
def plot_data_mark_mut_dot(data, sync_positions):
"""
绘制数据，并在 sync_positions 指定的位置添加垂直线标记（如同步头位置）
在第二个红线后以10为距离添加标记，直到下一个第一条红线为止，标记显示对应index
"""
plt.figure(figsize=(12, 6)) # 设置图表大小
plt.plot(data, marker='o', markersize=3, linewidth=1, label='Signal') # 绘制数据
# 先绘制所有红线
red_lines = []
for pos in sync_positions:
# 第一条红线位置
first_line = pos
# 第二条红线位置
second_line = pos + 30
plt.axvline(x=first_line, color='red', linestyle='--', linewidth=2, alpha=0.8)
plt.axvline(x=second_line, color='red', linestyle='--', linewidth=2, alpha=0.8)
red_lines.append((first_line, second_line))
# 为每个区间添加浮动文本标记
for i, (first_line, second_line) in enumerate(red_lines):
# 确定当前区间的结束位置（下一个第一条红线）
if i < len(red_lines) - 1:
end_pos = red_lines[i+1][0] # 下一个第一条红线位置
else:
end_pos = len(data) - 1 # 最后一个区间以数据长度为终点
# 从第二条红线后10个位置开始
start_pos = second_line + 10
# 确保起始位置不超过结束位置
if start_pos > end_pos:
continue
# 以10为间隔添加浮动标记（显示index）
current_pos = start_pos
while current_pos <= end_pos:
# 获取当前位置的数据值（用于定位文本）
data_value = data[current_pos]
# 添加浮动文本标记（显示当前index）
# 位置：在点的右下方（可通过xytext调整偏移）
# 仅第一个标记添加图例标签
label = 'Index Mark' if (i == 0 and current_pos == start_pos) else ""
if(current_pos in score_dict):
plt.text(
current_pos, data_value, # 文本锚点位置（对应数据点）
f'{score_dict[current_pos]}', # 显示的文本内容（index值）
fontsize=8, # 文本大小
color='blue', # 文本颜色
ha='left', va='bottom', # 水平/垂直对齐方式
bbox=dict(facecolor='white', alpha=0.3, boxstyle='round,pad=0.3'), # 白色背景框
label=label
)
# print("----",current_pos)
current_pos += 10
plt.title('Data Plot with Sync Markers') # 添加标题
plt.xlabel('Sample Index') # X轴标签
plt.ylabel('Amplitude') # Y轴标签
plt.grid(True, alpha=0.6) # 显示网格线
# 处理图例，确保只显示唯一标签
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys())
# 显示图表
plt.tight_layout()
plt.show()
# 输出数据（可选）
print("Data:", data)
print("Sync positions marked:", sync_positions)
# ========================
# 数据滤波
# ========================
def lowpass_filter_scipy(data, cutoff_freq, sampling_rate, filter_type='butter', order=4):
"""
使用SciPy对NumPy数组进行低通滤波
参数:
data: numpy数组，输入信号
cutoff_freq: float, 截止频率 (Hz)
sampling_rate: float, 采样频率 (Hz)
filter_type: str, 滤波器类型 ('butter', 'cheby1', 'bessel')
order: int, 滤波器阶数
返回:
filtered_data: 滤波后的numpy数组
"""
# 计算归一化截止频率 (Nyquist频率的倍数)
nyquist_freq = 0.5 * sampling_rate
normal_cutoff = cutoff_freq / nyquist_freq
if filter_type == 'butter':
b, a = signal.butter(order, normal_cutoff, btype='low', analog=False)
elif filter_type == 'cheby1':
b, a = signal.cheby1(order, 0.5, normal_cutoff, btype='low', analog=False)
elif filter_type == 'bessel':
b, a = signal.bessel(order, normal_cutoff, btype='low', analog=False)
else:
raise ValueError("不支持的滤波器类型")
# 应用滤波器
filtered_data = signal.filtfilt(b, a, data)
return filtered_data
def moving_average_filter(data, window_size):
"""
使用移动平均进行低通滤波
参数:
data: numpy数组
window_size: int, 窗口大小
返回:
filtered_data: 滤波后的数组
"""
if window_size <= 1:
return data.copy()
# 使用卷积实现移动平均
window = np.ones(window_size) / window_size
filtered_data = np.convolve(data, window, mode='same')
return filtered_data
res_dict = {}
score_dict = {}
if __name__ == "__main__":
sync_head_num = 0 # 查找到同步头的个数
paritiy_num = 0
paritiy_list = []
not_par_num = 0
not_par_list = []
data = np.loadtxt('***.dat', skiprows=1) # 加载数据文件，跳过第一行
# 预处理数据
data = -1*data
print(data)
startindex = 0
data = data[startindex:startindex+49000] # 100000
# 0.预处理
# data = moving_average_filter(data,2)
data = lowpass_filter_scipy(data,40000,100000)
# 1. 生成模板
template_0, template_1 = generate_manchester_templates()
sync_template = generate_sync_template(template_1, 3)
# 2. 波形数据
test_wave = data
# plot_data(test_wave)
data_len = data.size
index = 0
sync_list = [] # 标记同步头位置的列表
while(index<(data_len-1000)): # 留1000个缓冲空间
# 查找同步头
sync_positions = -1
window_size = 300
window_overlap = 100
while(sync_positions ==-1):
index = index + window_size - window_overlap # 窗口重叠
segment = test_wave[index:-1]
if(segment.size != 0):
sync_positions = find_sync_start_correlation(segment, window_size, sync_template, threshold= 55000.0)
else:
print("已遍历所有数据")
break
if(segment.size != 0):
sync_head_num = sync_head_num +1
print("第一个同步头的位置",index+sync_positions)
# plot_data_mark(test_wave[0:-1],index+sync_positions) # python 不支持函数重载
sync_list.append(index+sync_positions)
# 3. 解码位流
data = test_wave[index+sync_positions+30:index+sync_positions+30+170]
print("+++",index+sync_positions+30)
decoded_bits,r_list = decode_bits_from_waveform_r_peaks(data, template_0, template_1) # 解完17个就返回
res_dict[index+sync_positions] = decoded_bits
value = index+sync_positions # +30+10+0*10
score_dict[value] = r_list
# 4. 对list数据decoded_bits进行奇偶校验
is_valid = even_parity_check(decoded_bits)
print(f"奇偶校验{'通过' if is_valid else '未通过'}")
if is_valid:
paritiy_num = paritiy_num + 1
paritiy_list.append(index+sync_positions)
else:
not_par_num = not_par_num + 1
not_par_list.append(index+sync_positions)
print("-------------统计信息------------")
print("找到同步头",sync_head_num)
print("奇偶校验未通过个数",paritiy_num,paritiy_list)
print("奇偶校验通过个数",not_par_num,not_par_list)
print("---------------------------------")
print("-------------绘制同步头位置------------")
plot_data_mark_mut_dot(test_wave,sync_list)
print("---------------------------------")
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
