"""
干涉雷达多特征生命体征签名识别
论文: MDPI Sensors 2026, 26(18), 5724

核心: 四维特征签名 + ML验证
"""

import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F
from scipy import signal as sp_signal
from typing import Dict, List, Tuple

class RespiratoryHarmonicFeatures:
    """
    特征1: 呼吸谐波关系
    
    提取呼吸基频与各次谐波的幅度比、相位差
    
    论文: "respiratory harmonic relationships"
    """
    
    def __init__(self, num_harmonics: int = 5):
        self.num_harmonics = num_harmonics
    
    def extract(self, chest_signal: np.ndarray, 
                sample_rate: float = 100.0) -> np.ndarray:
        """
        提取呼吸谐波特征
        
        Args:
            chest_signal: 胸壁位移信号
            sample_rate: 采样率
            
        Returns:
            features: (num_harmonics * 2,) 谐波幅度比+相位差
        """
        # FFT
        freqs = np.fft.rfftfreq(len(chest_signal), 1/sample_rate)
        spectrum = np.fft.rfft(chest_signal)
        magnitude = np.abs(spectrum)
        phase = np.angle(spectrum)
        
        # 找呼吸基频
        mask = (freqs >= 0.1) & (freqs <= 0.5)
        resp_idx = np.argmax(magnitude[mask])
        resp_freq = freqs[mask][resp_idx]
        
        # 提取各次谐波
        features = []
        for h in range(1, self.num_harmonics + 1):
            h_freq = resp_freq * h
            h_idx = np.argmin(np.abs(freqs - h_freq))
            
            # 幅度比（相对基频）
            amp_ratio = magnitude[h_idx] / (magnitude[np.argmin(
                np.abs(freqs - resp_freq))] + 1e-8)
            
            # 相位差（相对基频）
            phase_diff = phase[h_idx] - phase[
                np.argmin(np.abs(freqs - resp_freq))]
            
            features.extend([amp_ratio, phase_diff])
        
        return np.array(features)


class InterHarmonicConsistency:
    """
    特征2: 谐波间一致性
    
    评估各次谐波之间的频率关系是否自洽
    
    论文: "inter-harmonic consistency"
    """
    
    def __init__(self, num_harmonics: int = 5):
        self.num_harmonics = num_harmonics
    
    def extract(self, chest_signal: np.ndarray,
                sample_rate: float = 100.0) -> np.ndarray:
        """
        提取谐波间一致性特征
        
        Returns:
            features: (num_harmonics - 1,) 一致性分数
        """
        freqs = np.fft.rfftfreq(len(chest_signal), 1/sample_rate)
        spectrum = np.abs(np.fft.rfft(chest_signal))
        
        # 找基频
        mask = (freqs >= 0.1) & (freqs <= 0.5)
        resp_idx = np.argmax(spectrum[mask])
        resp_freq = freqs[mask][resp_idx]
        
        features = []
        for h in range(2, self.num_harmonics + 1):
            h_freq = resp_freq * h
            h_idx = np.argmin(np.abs(freqs - h_freq))
            
            # 检查峰值附近是否有显著能量
            window = max(2, int(0.05 / (freqs[1] - freqs[0])))
            local_energy = spectrum[max(0, h_idx-window):h_idx+window+1]
            peak_prominence = spectrum[h_idx] / (local_energy.mean() + 1e-8)
            
            # 频率精确度：实际峰频 vs 预期谐波频率
            actual_peak = freqs[np.argmax(local_energy)]
            freq_accuracy = 1 - abs(actual_peak - h_freq) / h_freq
            
            features.append(peak_prominence * freq_accuracy)
        
        return np.array(features)


class ChestDisplacementSpectral:
    """
    特征3: 胸壁位移频谱特征
    
    频谱形状、能量分布、峭度等统计量
    
    论文: "chest-displacement spectral characteristics"
    """
    
    def extract(self, chest_signal: np.ndarray,
                sample_rate: float = 100.0) -> np.ndarray:
        """
        提取频谱统计特征
        
        Returns:
            features: (8,) 频谱统计量
        """
        spectrum = np.abs(np.fft.rfft(chest_signal))
        spectrum_norm = spectrum / (spectrum.sum() + 1e-8)
        
        freqs = np.fft.rfftfreq(len(chest_signal), 1/sample_rate)
        
        # 频谱质心
        centroid = np.sum(freqs * spectrum_norm)
        
        # 频谱带宽
        bandwidth = np.sqrt(np.sum(((freqs - centroid) ** 2) * spectrum_norm))
        
        # 频谱平坦度
        geometric_mean = np.exp(np.mean(np.log(spectrum + 1e-8)))
        arithmetic_mean = spectrum.mean()
        flatness = geometric_mean / (arithmetic_mean + 1e-8)
        
        # 频谱峭度
        kurtosis = np.sum(((freqs - centroid) ** 4) * spectrum_norm) / (
            bandwidth ** 4 + 1e-8)
        
        # 低频能量比 (0-1 Hz)
        low_freq_energy = spectrum[freqs < 1.0].sum() / (spectrum.sum() + 1e-8)
        
        # 高频能量比 (>1 Hz)
        high_freq_energy = spectrum[freqs > 1.0].sum() / (spectrum.sum() + 1e-8)
        
        # 总能量
        total_energy = np.sum(chest_signal ** 2) / len(chest_signal)
        
        # 峰值频率
        peak_freq = freqs[np.argmax(spectrum)]
        
        return np.array([
            centroid, bandwidth, flatness, kurtosis,
            low_freq_energy, high_freq_energy, total_energy, peak_freq
        ])


class CardiacMechanicalSignature:
    """
    特征4: 心机机械信号
    
    从胸壁位移中提取心脏机械活动引起的微弱信号
    
    论文: "radar-derived cardiac mechanical signatures"
    """
    
    def __init__(self, sample_rate: float = 100.0):
        self.sample_rate = sample_rate
    
    def extract(self, chest_signal: np.ndarray) -> np.ndarray:
        """
        提取心机机械特征
        
        Returns:
            features: (6,) 心机机械特征
        """
        # 高通滤波分离心跳
        b, a = sp_signal.butter(4, 0.8, btype='high', 
                                 fs=self.sample_rate)
        cardiac = sp_signal.filtfilt(b, a, chest_signal)
        
        # 心跳峰值检测
        peaks, _ = sp_signal.find_peaks(
            cardiac,
            distance=int(0.4 * self.sample_rate),
            prominence=0.5 * np.std(cardiac)
        )
        
        if len(peaks) < 3:
            return np.zeros(6)
        
        # 心跳间隔 (IBI)
        ibi = np.diff(peaks) / self.sample_rate * 1000  # ms
        
        # 特征
        mean_hr = 60 / (np.mean(ibi) / 1000) if np.mean(ibi) > 0 else 0
        std_hr = np.std(ibi)
        hrv_rmssd = np.sqrt(np.mean(np.diff(ibi) ** 2))  # RMSSD
        hrv_sdnn = np.std(ibi)  # SDNN
        
        # 心跳幅度
        amplitudes = cardiac[peaks]
        mean_amp = np.mean(amplitudes)
        std_amp = np.std(amplitudes)
        
        return np.array([
            mean_hr, std_hr, hrv_rmssd, hrv_sdnn,
            mean_amp, std_amp
        ])


class VitalSignMLValidator(nn.Module):
    """
    机器学习验证器
    
    将四维特征签名输入ML模型进行人体检测验证
    
    论文: "ML Validation for Privacy-Preserving Human Detection"
    """
    
    def __init__(self, 
                 harmonic_dim: int = 10,
                 consistency_dim: int = 4,
                 spectral_dim: int = 8,
                 cardiac_dim: int = 6,
                 hidden_dim: int = 128):
        super().__init__()
        input_dim = harmonic_dim + consistency_dim + spectral_dim + cardiac_dim
        
        self.classifier = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.BatchNorm1d(hidden_dim),
            nn.ReLU(),
            nn.Dropout(0.3),
            nn.Linear(hidden_dim, hidden_dim // 2),
            nn.BatchNorm1d(hidden_dim // 2),
            nn.ReLU(),
            nn.Dropout(0.3),
            nn.Linear(hidden_dim // 2, 3),  # 3类:无人/有人/不确定
        )
    
    def forward(self, features: torch.Tensor) -> Dict:
        """
        Args:
            features: (B, feature_dim) 四维特征签名
            
        Returns:
            prediction: 人体检测结果
        """
        logits = self.classifier(features)
        probs = F.softmax(logits, dim=-1)
        
        idx = probs.argmax(dim=-1)
        labels = ['no_person', 'person_present', 'uncertain']
        
        results = []
        for i in range(features.shape[0]):
            results.append({
                'detection': labels[idx[i].item()],
                'confidence': probs[i, idx[i]].item(),
                'probabilities': {
                    'no_person': probs[i, 0].item(),
                    'person_present': probs[i, 1].item(),
                    'uncertain': probs[i, 2].item(),
                }
            })
        
        return {'results': results}


# ========== 完整系统 ==========
class InterferometricVitalSignSystem:
    """
    干涉雷达多特征生命体征签名识别系统
    
    论文: MDPI Sensors 2026, 26(18), 5724
    
    四维特征:
    1. 呼吸谐波关系
    2. 谐波间一致性
    3. 胸壁位移频谱特征
    4. 心机机械信号
    
    ML验证: 人体检测+隐私保护
    """
    
    def __init__(self):
        self.harmonic = RespiratoryHarmonicFeatures(num_harmonics=5)
        self.consistency = InterHarmonicConsistency(num_harmonics=5)
        self.spectral = ChestDisplacementSpectral()
        self.cardiac = CardiacMechanicalSignature()
        self.validator = VitalSignMLValidator()
        
    def process(self, radar_signal: np.ndarray,
                sample_rate: float = 100.0) -> Dict:
        """
        完整处理流程
        
        Args:
            radar_signal: 雷达解调后的胸壁位移信号
            sample_rate: 采样率
            
        Returns:
            result: 检测结果+特征签名
        """
        # 提取四维特征
        f1 = self.harmonic.extract(radar_signal, sample_rate)
        f2 = self.consistency.extract(radar_signal, sample_rate)
        f3 = self.spectral.extract(radar_signal, sample_rate)
        f4 = self.cardiac.extract(radar_signal)
        
        # 拼接特征
        feature_vector = np.concatenate([f1, f2, f3, f4])
        
        # ML验证
        features_tensor = torch.from_numpy(feature_vector).float().unsqueeze(0)
        with torch.no_grad():
            result = self.validator(features_tensor)
        
        return {
            'detection': result['results'][0],
            'feature_signature': {
                'harmonic_relations': f1.tolist(),
                'inter_harmonic_consistency': f2.tolist(),
                'spectral_features': f3.tolist(),
                'cardiac_mechanical': f4.tolist(),
            },
            'feature_dim': len(feature_vector),
        }


# ========== 测试 ==========
if __name__ == "__main__":
    np.random.seed(42)
    torch.manual_seed(42)
    
    system = InterferometricVitalSignSystem()
    
    # 模拟有人信号
    sample_rate = 100.0
    duration = 30.0
    N = int(duration * sample_rate)
    t = np.arange(N) / sample_rate
    
    # 呼吸 + 心跳 + 噪声
    resp = 6.0 * np.sin(2 * np.pi * 0.25 * t)
    cardiac = 1.2 * np.sin(2 * np.pi * 1.2 * t)
    noise = np.random.normal(0, 0.2, N)
    person_signal = resp + cardiac + noise
    
    # 模拟无人信号（只有噪声）
    empty_signal = np.random.normal(0, 0.1, N)
    
    print("=" * 60)
    print("干涉雷达多特征生命体征签名识别测试")
    
    # 测试1: 有人
    print("\n场景1: 有生命体征存在")
    result1 = system.process(person_signal, sample_rate)
    print(f"  检测: {result1['detection']['detection']}")
    print(f"  置信度: {result1['detection']['confidence']:.3f}")
    print(f"  特征维度: {result1['feature_dim']}")
    sig = result1['feature_signature']
    print(f"  呼吸谐波: {len(sig['harmonic_relations'])}维")
    print(f"  谐波一致性: {len(sig['inter_harmonic_consistency'])}维")
    print(f"  频谱特征: {len(sig['spectral_features'])}维")
    print(f"  心机机械: {len(sig['cardiac_mechanical'])}维")
    
    # 测试2: 无人
    print("\n场景2: 无生命体征")
    result2 = system.process(empty_signal, sample_rate)
    print(f"  检测: {result2['detection']['detection']}")
    print(f"  置信度: {result2['detection']['confidence']:.3f}")
    
    print("\n✅ 测试通过")