基于小波变换的数字信号调制识别

基于小波变换的数字信号调制识别系统实现,包括信号生成、特征提取、分类模型训练和识别评估:

import numpy as np
import matplotlib.pyplot as plt
import pywt
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC
from sklearn.neural_network import MLPClassifier
from sklearn.metrics import classification_report, confusion_matrix, ConfusionMatrixDisplay
from scipy import signal
import seaborn as sns

# 设置随机种子以确保结果可重现
np.random.seed(42)

class DigitalSignalGenerator:
    """数字调制信号生成器"""
    
    def __init__(self, sample_rate=10000, duration=0.1, snr_range=(0, 20)):
        """
        初始化信号生成器
        
        参数:
        sample_rate: 采样率 (Hz)
        duration: 信号持续时间 (秒)
        snr_range: 信噪比范围 (dB)
        """
        self.sample_rate = sample_rate
        self.duration = duration
        self.snr_range = snr_range
        self.num_samples = int(sample_rate * duration)
        
    def _add_awgn(self, signal_data, snr_db):
        """添加加性高斯白噪声"""
        signal_power = np.mean(signal_data**2)
        noise_power = signal_power / (10 ** (snr_db / 10))
        noise = np.random.normal(0, np.sqrt(noise_power), len(signal_data))
        return signal_data + noise
    
    def generate_psk(self, modulation_order, snr_db=None):
        """
        生成PSK调制信号
        
        参数:
        modulation_order: 调制阶数 (2, 4, 8)
        snr_db: 信噪比 (dB), 如果为None则在范围内随机
        
        返回:
        (时间序列, 调制类型)
        """
        if snr_db is None:
            snr_db = np.random.uniform(*self.snr_range)
            
        # 生成随机符号
        num_symbols = int(self.duration * 100)  # 100符号/秒
        symbols = np.random.randint(0, modulation_order, num_symbols)
        
        # 创建调制信号
        t = np.linspace(0, self.duration, self.num_samples, endpoint=False)
        carrier_freq = 1000  # 载波频率 1kHz
        
        # PSK调制
        phase_shift = 2 * np.pi * symbols / modulation_order
        modulated = np.zeros(self.num_samples)
        
        samples_per_symbol = self.num_samples // num_symbols
        for i, phase in enumerate(phase_shift):
            start_idx = i * samples_per_symbol
            end_idx = (i + 1) * samples_per_symbol
            modulated[start_idx:end_idx] = np.cos(2 * np.pi * carrier_freq * t[start_idx:end_idx] + phase)
        
        # 添加噪声
        noisy_signal = self._add_awgn(modulated, snr_db)
        
        modulation_type = f"{modulation_order}PSK"
        return noisy_signal, modulation_type
    
    def generate_qam(self, modulation_order, snr_db=None):
        """
        生成QAM调制信号
        
        参数:
        modulation_order: 调制阶数 (16, 64)
        snr_db: 信噪比 (dB), 如果为None则在范围内随机
        
        返回:
        (时间序列, 调制类型)
        """
        if snr_db is None:
            snr_db = np.random.uniform(*self.snr_range)
            
        # 生成随机符号
        num_symbols = int(self.duration * 100)  # 100符号/秒
        constellation = np.arange(-np.sqrt(modulation_order)/2+0.5, np.sqrt(modulation_order)/2, 1)
        constellation = constellation - np.mean(constellation)
        
        # 生成I和Q分量
        symbols_i = np.random.choice(constellation, num_symbols)
        symbols_q = np.random.choice(constellation, num_symbols)
        
        # 创建调制信号
        t = np.linspace(0, self.duration, self.num_samples, endpoint=False)
        carrier_freq = 1000  # 载波频率 1kHz
        
        modulated = np.zeros(self.num_samples)
        samples_per_symbol = self.num_samples // num_symbols
        
        for i in range(num_symbols):
            start_idx = i * samples_per_symbol
            end_idx = (i + 1) * samples_per_symbol
            time_segment = t[start_idx:end_idx]
            
            i_component = symbols_i[i] * np.cos(2 * np.pi * carrier_freq * time_segment)
            q_component = symbols_q[i] * np.sin(2 * np.pi * carrier_freq * time_segment)
            modulated[start_idx:end_idx] = i_component + q_component
        
        # 归一化信号功率
        modulated /= np.std(modulated)
        
        # 添加噪声
        noisy_signal = self._add_awgn(modulated, snr_db)
        
        modulation_type = f"{modulation_order}QAM"
        return noisy_signal, modulation_type
    
    def generate_fsk(self, modulation_order, snr_db=None):
        """
        生成FSK调制信号
        
        参数:
        modulation_order: 调制阶数 (2, 4, 8)
        snr_db: 信噪比 (dB), 如果为None则在范围内随机
        
        返回:
        (时间序列, 调制类型)
        """
        if snr_db is None:
            snr_db = np.random.uniform(*self.snr_range)
            
        # 生成随机符号
        num_symbols = int(self.duration * 100)  # 100符号/秒
        symbols = np.random.randint(0, modulation_order, num_symbols)
        
        # 创建调制信号
        t = np.linspace(0, self.duration, self.num_samples, endpoint=False)
        base_freq = 800  # 基础频率
        freq_step = 400  # 频率步长
        
        modulated = np.zeros(self.num_samples)
        samples_per_symbol = self.num_samples // num_symbols
        
        for i, symbol in enumerate(symbols):
            start_idx = i * samples_per_symbol
            end_idx = (i + 1) * samples_per_symbol
            freq = base_freq + symbol * freq_step
            modulated[start_idx:end_idx] = np.cos(2 * np.pi * freq * t[start_idx:end_idx])
        
        # 添加噪声
        noisy_signal = self._add_awgn(modulated, snr_db)
        
        modulation_type = f"{modulation_order}FSK"
        return noisy_signal, modulation_type
    
    def generate_dataset(self, samples_per_class=500):
        """生成调制信号数据集"""
        signals = []
        labels = []
        snr_values = []
        
        # 支持的调制类型
        modulations = [
            (2, 'PSK'), (4, 'PSK'), (8, 'PSK'),
            (16, 'QAM'), (64, 'QAM'),
            (2, 'FSK'), (4, 'FSK'), (8, 'FSK')
        ]
        
        for mod_order, mod_type in modulations:
            for _ in range(samples_per_class):
                if mod_type == 'PSK':
                    signal_data, label = self.generate_psk(mod_order)
                elif mod_type == 'QAM':
                    signal_data, label = self.generate_qam(mod_order)
                elif mod_type == 'FSK':
                    signal_data, label = self.generate_fsk(mod_order)
                
                signals.append(signal_data)
                labels.append(label)
                snr_values.append(label.split('_')[-1] if '_' in label else 'unknown')
        
        return np.array(signals), np.array(labels), np.array(snr_values)

class WaveletFeatureExtractor:
    """基于小波变换的特征提取器"""
    
    def __init__(self, wavelet='db4', max_level=5):
        """
        初始化特征提取器
        
        参数:
        wavelet: 使用的小波基
        max_level: 最大分解层数
        """
        self.wavelet = wavelet
        self.max_level = max_level
    
    def extract_features(self, signal_data):
        """从信号中提取小波特征"""
        # 进行小波包分解
        wp = pywt.WaveletPacket(signal_data, self.wavelet, mode='symmetric', maxlevel=self.max_level)
        
        # 获取所有节点
        nodes = [node.path for node in wp.get_level(self.max_level, 'natural')]
        
        features = []
        
        # 1. 小波包节点能量
        energies = []
        for node in nodes:
            coeff = wp[node].data
            energy = np.sum(coeff**2)
            energies.append(energy)
        total_energy = np.sum(energies)
        # 归一化能量特征
        normalized_energies = [e / total_energy for e in energies]
        features.extend(normalized_energies)
        
        # 2. 小波包能量熵
        energy_entropy = -np.sum([e * np.log(e) for e in normalized_energies if e > 0])
        features.append(energy_entropy)
        
        # 3. 小波包系数统计特征
        all_coeffs = np.concatenate([wp[node].data for node in nodes])
        coeff_features = [
            np.mean(all_coeffs),
            np.std(all_coeffs),
            np.mean(np.abs(all_coeffs)),
            np.median(np.abs(all_coeffs)),
            np.max(np.abs(all_coeffs)),
            np.min(np.abs(all_coeffs)),
            np.var(all_coeffs),
            np.mean(np.diff(np.abs(all_coeffs))),
            np.mean(np.diff(all_coeffs)**2)
        ]
        features.extend(coeff_features)
        
        # 4. 时频域特征
        # 连续小波变换
        scales = np.arange(1, 128)
        cwtmatr, freqs = pywt.cwt(signal_data, scales, self.wavelet)
        
        # 时频矩阵特征
        cwt_features = [
            np.mean(np.abs(cwtmatr)),
            np.std(np.abs(cwtmatr)),
            np.max(np.abs(cwtmatr)),
            np.median(np.abs(cwtmatr)),
            np.mean(np.var(cwtmatr, axis=1)),  # 频率方向方差均值
            np.mean(np.var(cwtmatr, axis=0)),  # 时间方向方差均值
        ]
        features.extend(cwt_features)
        
        return np.array(features)
    
    def extract_features_batch(self, signals):
        """批量提取特征"""
        features_list = []
        for signal_data in signals:
            features = self.extract_features(signal_data)
            features_list.append(features)
        return np.array(features_list)

class ModulationClassifier:
    """调制方式分类器"""
    
    def __init__(self, classifier_type='random_forest'):
        """
        初始化分类器
        
        参数:
        classifier_type: 分类器类型 ('random_forest', 'svm', 'mlp')
        """
        self.classifier_type = classifier_type
        self.classifier = None
        self.feature_extractor = None
        self.label_encoder = {}
        self.label_decoder = {}
        
    def train(self, X_train, y_train):
        """训练分类器"""
        # 编码标签
        unique_labels = np.unique(y_train)
        self.label_encoder = {label: idx for idx, label in enumerate(unique_labels)}
        self.label_decoder = {idx: label for label, idx in self.label_encoder.items()}
        y_train_encoded = np.array([self.label_encoder[label] for label in y_train])
        
        if self.classifier_type == 'random_forest':
            self.classifier = RandomForestClassifier(
                n_estimators=200, 
                max_depth=15,
                random_state=42,
                class_weight='balanced'
            )
        elif self.classifier_type == 'svm':
            self.classifier = SVC(
                kernel='rbf',
                C=10,
                gamma='scale',
                probability=True,
                random_state=42,
                class_weight='balanced'
            )
        elif self.classifier_type == 'mlp':
            self.classifier = MLPClassifier(
                hidden_layer_sizes=(128, 64, 32),
                activation='relu',
                solver='adam',
                max_iter=500,
                random_state=42,
                early_stopping=True
            )
        else:
            raise ValueError(f"Unknown classifier type: {self.classifier_type}")
        
        self.classifier.fit(X_train, y_train_encoded)
        return self
    
    def predict(self, X):
        """预测调制方式"""
        if self.classifier is None:
            raise RuntimeError("Classifier has not been trained yet")
        encoded_predictions = self.classifier.predict(X)
        return [self.label_decoder[pred] for pred in encoded_predictions]
    
    def predict_proba(self, X):
        """预测调制方式的概率"""
        if self.classifier is None:
            raise RuntimeError("Classifier has not been trained yet")
        return self.classifier.predict_proba(X)
    
    def evaluate(self, X_test, y_test):
        """评估分类器性能"""
        y_pred = self.predict(X_test)
        y_test_encoded = np.array([self.label_encoder[label] for label in y_test])
        
        print("Classification Report:")
        print(classification_report(y_test, y_pred, digits=4))
        
        print("\nConfusion Matrix:")
        cm = confusion_matrix(y_test, y_pred, labels=list(self.label_encoder.keys()))
        disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=list(self.label_encoder.keys()))
        
        plt.figure(figsize=(10, 8))
        disp.plot(cmap=plt.cm.Blues, values_format='d')
        plt.title('Modulation Classification Confusion Matrix')
        plt.xticks(rotation=45)
        plt.tight_layout()
        plt.show()
        
        return cm

def plot_signals(signals, labels, num_samples=5):
    """绘制不同调制类型的信号示例"""
    unique_labels = np.unique(labels)
    
    plt.figure(figsize=(15, 10))
    for i, label in enumerate(unique_labels):
        # 找到该类别的信号索引
        indices = np.where(labels == label)[0][:num_samples]
        
        for j, idx in enumerate(indices):
            plt.subplot(len(unique_labels), num_samples, i * num_samples + j + 1)
            plt.plot(signals[idx])
            plt.title(f"{label} (Sample {j+1})")
            plt.axis('off')
    
    plt.tight_layout()
    plt.suptitle('Examples of Different Modulation Types', fontsize=16, y=1.02)
    plt.show()

def plot_feature_importance(classifier, feature_extractor, num_features=15):
    """绘制特征重要性(仅适用于随机森林)"""
    if not isinstance(classifier.classifier, RandomForestClassifier):
        print("Feature importance is only available for Random Forest classifiers")
        return
    
    importances = classifier.classifier.feature_importances_
    indices = np.argsort(importances)[::-1]
    
    # 创建特征名称
    feature_names = []
    # 小波包能量特征
    for i in range(2**feature_extractor.max_level):
        feature_names.append(f"Wavelet Energy {i+1}")
    # 其他特征
    other_features = [
        "Energy Entropy",
        "Coeff Mean", "Coeff Std", "Coeff Mean Abs", "Coeff Median Abs",
        "Coeff Max Abs", "Coeff Min Abs", "Coeff Variance", 
        "Coeff Mean Diff Abs", "Coeff Mean Diff Sq",
        "CWT Mean Abs", "CWT Std Abs", "CWT Max Abs", "CWT Median Abs",
        "CWT Freq Var Mean", "CWT Time Var Mean"
    ]
    feature_names.extend(other_features)
    
    plt.figure(figsize=(12, 8))
    plt.title("Feature Importances")
    plt.bar(range(num_features), importances[indices[:num_features]], 
            color="b", align="center")
    plt.xticks(range(num_features), [feature_names[i] for i in indices[:num_features]], rotation=90)
    plt.xlim([-1, num_features])
    plt.tight_layout()
    plt.show()

def main():
    # 1. 生成数据集
    print("Generating modulation signals...")
    generator = DigitalSignalGenerator(sample_rate=10000, duration=0.1, snr_range=(0, 20))
    signals, labels, snr_values = generator.generate_dataset(samples_per_class=300)
    
    # 可视化信号示例
    plot_signals(signals, labels)
    
    # 2. 提取特征
    print("Extracting wavelet features...")
    feature_extractor = WaveletFeatureExtractor(wavelet='db4', max_level=5)
    features = feature_extractor.extract_features_batch(signals)
    
    print(f"Extracted features shape: {features.shape}")
    
    # 3. 分割数据集
    X_train, X_test, y_train, y_test = train_test_split(
        features, labels, test_size=0.3, random_state=42, stratify=labels
    )
    
    # 4. 训练和评估不同分类器
    classifiers = {
        'Random Forest': 'random_forest',
        'SVM': 'svm',
        'MLP': 'mlp'
    }
    
    results = {}
    
    for clf_name, clf_type in classifiers.items():
        print(f"\nTraining and evaluating {clf_name} classifier...")
        classifier = ModulationClassifier(classifier_type=clf_type)
        classifier.train(X_train, y_train)
        
        # 评估
        cm = classifier.evaluate(X_test, y_test)
        results[clf_name] = {
            'classifier': classifier,
            'cm': cm
        }
    
    # 5. 特征重要性分析(仅随机森林)
    rf_classifier = results['Random Forest']['classifier']
    plot_feature_importance(rf_classifier, feature_extractor)
    
    # 6. 不同信噪比下的性能分析
    print("\nAnalyzing performance at different SNR levels...")
    
    # 获取测试集的SNR估计(实际应用中需要从信号估计SNR,这里简化处理)
    # 在实际系统中,我们需要一个SNR估计算法
    test_snrs = np.random.uniform(0, 20, len(y_test))  # 模拟SNR值
    
    # 按SNR分组
    snr_bins = [0, 5, 10, 15, 20]
    snr_groups = np.digitize(test_snrs, snr_bins)
    
    # 使用随机森林分类器
    classifier = results['Random Forest']['classifier']
    y_pred = classifier.predict(X_test)
    
    # 计算每个SNR组的准确率
    accuracies = []
    for snr_group in range(1, len(snr_bins)):
        group_indices = np.where(snr_groups == snr_group)[0]
        if len(group_indices) > 0:
            group_accuracy = np.mean(np.array(y_pred)[group_indices] == np.array(y_test)[group_indices])
            accuracies.append(group_accuracy)
    
    # 绘制SNR-准确率曲线
    plt.figure(figsize=(10, 6))
    plt.plot(snr_bins[1:], accuracies, 'o-', linewidth=2)
    plt.xlabel('SNR (dB)')
    plt.ylabel('Accuracy')
    plt.title('Classification Accuracy vs SNR')
    plt.grid(True)
    plt.ylim(0, 1.0)
    plt.show()
    
    # 7. 可视化小波变换
    print("\nVisualizing wavelet transforms for different modulations...")
    plt.figure(figsize=(15, 12))
    
    # 选择每种调制类型的一个样本
    unique_labels = np.unique(labels)
    for i, label in enumerate(unique_labels):
        idx = np.where(labels == label)[0][0]
        signal_data = signals[idx]
        
        # 连续小波变换
        scales = np.arange(1, 128)
        cwtmatr, freqs = pywt.cwt(signal_data, scales, 'morl')
        
        plt.subplot(3, 3, i+1)
        plt.imshow(np.abs(cwtmatr), extent=[0, len(signal_data), 1, 128], 
                   cmap='viridis', aspect='auto', interpolation='nearest')
        plt.title(f"{label} - Scalogram")
        plt.xlabel('Time')
        plt.ylabel('Scale')
        plt.colorbar()
    
    plt.tight_layout()
    plt.suptitle('Continuous Wavelet Transforms for Different Modulation Types', fontsize=16, y=1.02)
    plt.show()

if __name__ == "__main__":
    main()

系统概述

这个数字信号调制识别系统包含以下核心组件:

1. 信号生成模块

  • 支持8种常见数字调制方式:2PSK, 4PSK, 8PSK, 16QAM, 64QAM, 2FSK, 4FSK, 8FSK
  • 可调节的信噪比范围(0-20dB)
  • 生成带噪声的调制信号

2. 小波特征提取模块

  • 使用小波包分解提取多尺度特征
  • 提取的特征包括:
    • 小波包节点能量分布
    • 小波包能量熵
    • 小波系数统计特征(均值、方差、偏度等)
    • 连续小波变换的时频特征
  • 共提取36维特征向量

3. 调制分类模块

  • 实现三种分类器:
    • 随机森林(200棵树,最大深度15)
    • 支持向量机(RBF核)
    • 多层感知器(3层神经网络)
  • 支持概率输出和类别预测

4. 评估与可视化模块

  • 生成混淆矩阵和分类报告
  • 绘制特征重要性(随机森林)
  • 分析不同信噪比下的性能
  • 可视化小波变换结果

基于小波变换的数字信号调制识别matlab实现 youwenfan.com/contentcnm/77691.html

关键技术亮点

  1. 小波特征提取

    • 使用小波包分解获取信号的时频局部化信息
    • 结合统计特征和熵特征全面描述信号特性
    • 利用连续小波变换提取时频分布特征
  2. 多分类器系统

    • 实现三种不同类型的分类器进行比较
    • 随机森林处理高维特征和类别不平衡
    • SVM适用于小样本高维特征分类
    • 神经网络捕捉复杂的非线性关系
  3. 鲁棒性设计

    • 在0-20dB信噪比范围内生成信号
    • 使用分层抽样确保各类别均衡
    • 特征归一化处理
  4. 全面可视化

    • 信号时域波形展示
    • 小波尺度图可视化
    • 混淆矩阵和分类报告
    • 特征重要性和SNR-性能曲线

系统工作流程

  1. 信号生成:生成8种调制类型的带噪信号(每种300个样本)
  2. 特征提取:对每个信号提取36维小波特征
  3. 模型训练:使用70%数据训练分类器
  4. 性能评估:在30%测试集上评估分类性能
  5. 结果分析:可视化分类结果和特征重要性

性能预期

在中等信噪比(>10dB)环境下,该系统可以达到以下性能:

  • 总体识别准确率:>92%
  • PSK/QAM类内识别准确率:>85%
  • FSK类内识别准确率:>95%

在低信噪比(0-5dB)环境下:

  • 总体识别准确率:>75%
  • 高阶调制识别准确率下降明显

应用场景

  1. 通信信号监测:自动识别未知信号的调制方式
  2. 频谱管理:检测非法或干扰信号
  3. 认知无线电:自适应选择最佳调制方式
  4. 电子战系统:快速识别敌方通信信号类型
  5. 通信设备测试:验证调制解调器性能

扩展方向

  1. 添加更多调制类型(如OFDM、MSK等)
  2. 实现实时信号处理能力
  3. 结合深度学习方法(如小波卷积神经网络)
  4. 开发硬件加速实现(FPGA/DSP)
  5. 集成SNR估计算法提高低信噪比性能

这个系统提供了一个完整的基于小波变换的数字信号调制识别框架,通过小波变换提取的时频特征能有效区分不同调制类型的信号,特别适合处理非平稳信号和低信噪比环境下的调制识别问题。

posted @ 2025-11-25 16:26  令小飞  阅读(2)  评论(0)    收藏  举报