"""
goVISCAN 多光谱 rPPG（远程光电容积描记）原理复现

核心原理：
1. 血液脉动导致皮肤亮度微小变化
2. 不同波长对血液的吸收/反射特性不同
3. 多通道同时采集 → 信噪比提升 → 运动鲁棒性

对比传统单光谱 rPPG：
- 传统：RGB 三通道，受环境光照影响大
- goVISCAN：9 通道（含近红外），光照鲁棒性强
"""

import numpy as np
from typing import Tuple


class MultispectralRPPG:
    """
    多光谱 rPPG 信号提取
    
    模拟 goVISCAN 的 9 通道多光谱采集
    """
    
    # 9 个通道的中心波长 (nm)
    WAVELENGTHS = {
        'ch1': 450,   # 蓝
        'ch2': 470,   # 蓝
        'ch3': 505,   # 青绿
        'ch4': 530,   # 绿
        'ch5': 565,   # 黄绿
        'ch6': 585,   # 橙
        'ch7': 625,   # 红
        'ch8': 660,   # 深红
        'ch9': 850,   # 近红外
    }
    
    # 血红蛋白吸收系数（近似值，单位 cm^-1）
    # HbO2（氧合血红蛋白）和 Hb（脱氧血红蛋白）
    HEMOGLOBIN_ABSORPTION = {
        'ch1': (8.0, 8.5),   # (HbO2, Hb)
        'ch2': (7.0, 7.5),
        'ch3': (4.5, 5.0),
        'ch4': (3.0, 3.5),
        'ch5': (2.0, 2.5),
        'ch6': (1.5, 2.0),
        'ch7': (0.8, 1.2),
        'ch8': (0.5, 1.0),
        'ch9': (0.3, 0.2),  # 近红外：HbO2 > Hb（反转）
    }
    
    def __init__(self, fps=30, num_channels=9):
        self.fps = fps
        self.num_channels = num_channels
    
    def simulate_ppg_signal(self, heart_rate=72, 
                             breathing_rate=16,
                             spo2=98, duration=10):
        """
        模拟多光谱 PPG 信号
        
        Args:
            heart_rate: 心率 (bpm)
            breathing_rate: 呼吸率 (次/分)
            spo2: 血氧饱和度 (%)
            duration: 持续时间 (秒)
            
        Returns:
            signals: shape=(num_channels, num_samples)
        """
        num_samples = int(self.fps * duration)
        t = np.linspace(0, duration, num_samples)
        
        # 心率基频
        hr_freq = heart_rate / 60.0  # Hz
        # 呼吸基频
        br_freq = breathing_rate / 60.0  # Hz
        
        signals = np.zeros((self.num_channels, num_samples))
        
        for i, (ch, wl) in enumerate(self.WAVELENGTHS.items()):
            # 血红蛋白吸收
            hbo2_abs, hb_abs = self.HEMOGLOBIN_ABSORPTION[ch]
            
            # 氧合度对吸收的影响
            actual_hbo2_abs = hbo2_abs * (spo2 / 100.0)
            actual_hb_abs = hb_abs * ((100 - spo2) / 100.0)
            total_abs = actual_hbo2_abs + actual_hb_abs
            
            # PPG 信号 = 心率脉动 + 呼吸调制 + 噪声
            cardiac = np.sin(2 * np.pi * hr_freq * t)
            respiratory = 0.3 * np.sin(2 * np.pi * br_freq * t)
            noise = 0.1 * np.random.randn(num_samples)
            
            # 通道间差异：不同波长对血液体积变化敏感度不同
            sensitivity = total_abs / max(
                self.HEMOGLOBIN_ABSORPTION['ch4']
            )
            
            signals[i, :] = (
                sensitivity * (cardiac + respiratory) + noise
            )
        
        return signals
    
    def extract_heart_rate(self, signals):
        """
        从多光谱信号提取心率
        
        使用多通道融合 + FFT
        
        Args:
            signals: shape=(num_channels, num_samples)
            
        Returns:
            hr: 心率 (bpm)
        """
        # 方法：对每个通道做 FFT，然后融合
        freqs = np.fft.rfftfreq(signals.shape[1], 1/self.fps)
        
        # 去趋势
        signals_detrended = signals - np.mean(
            signals, axis=1, keepdims=True
        )
        
        # 每个通道的频谱
        spectra = np.abs(
            np.fft.rfft(signals_detrended, axis=1)
        )
        
        # 多通道融合：取加权平均频谱
        # 权重：绿光通道（对血液最敏感）
        weights = np.ones(self.num_channels)
        weights[3] = 2.0  # ch4 (530nm 绿)
        weights[4] = 1.5  # ch5
        weights /= weights.sum()
        
        fused_spectrum = np.average(spectra, 
                                    axis=0, 
                                    weights=weights)
        
        # 在心率范围 (0.7-3.5 Hz = 42-210 bpm) 扻最大值
        hr_mask = (freqs > 0.7) & (freqs < 3.5)
        hr_indices = np.where(hr_mask)[0]
        
        if len(hr_indices) > 0:
            peak_idx = hr_indices[
                np.argmax(fused_spectrum[hr_indices])
            ]
            hr = freqs[peak_idx] * 60
        else:
            hr = 0
        
        return hr
    
    def extract_breathing_rate(self, signals):
        """
        从多光谱信号提取呼吸率
        
        呼吸信号通过低频幅度调制传递
        
        Returns:
            br: 呼吸率 (次/分)
        """
        # 提取低频包络
        from scipy.signal import hilbert
        
        # 取绿光通道（对血液最敏感）
        green_signal = signals[3, :]
        
        # Hilbert 变换取包络
        analytic = hilbert(green_signal - np.mean(green_signal))
        envelope = np.abs(analytic)
        
        # FFT 找呼吸频率
        freqs = np.fft.rfftfreq(len(envelope), 1/self.fps)
        spectrum = np.abs(np.fft.rfft(envelope - np.mean(envelope)))
        
        # 呼吸范围 (0.1-0.5 Hz = 6-30 次/分)
        br_mask = (freqs > 0.1) & (freqs < 0.5)
        br_indices = np.where(br_mask)[0]
        
        if len(br_indices) > 0:
            peak_idx = br_indices[
                np.argmax(spectrum[br_indices])
            ]
            br = freqs[peak_idx] * 60
        else:
            br = 0
        
        return br
    
    def estimate_spo2(self, signals):
        """
        估算血氧饱和度
        
        利用红光和近红外通道的比值
        SpO2 = f(AC_red/DC_red) / (AC_nir/DC_nir)
        
        Returns:
            spo2: 血氧饱和度 (%)
        """
        # 红光通道 (660nm)
        red_signal = signals[7, :]
        # 近红外通道 (850nm)
        nir_signal = signals[8, :]
        
        # AC 分量（脉动）
        red_ac = np.std(red_signal)
        nir_ac = np.std(nir_signal)
        
        # DC 分量（直流）
        red_dc = np.mean(red_signal)
        nir_dc = np.mean(nir_signal)
        
        # 比值
        if red_dc > 0 and nir_dc > 0 and nir_ac > 0:
            ratio = (red_ac / red_dc) / (nir_ac / nir_dc)
            # 经验公式（简化版）
            spo2 = 110 - 25 * ratio
            spo2 = np.clip(spo2, 70, 100)
        else:
            spo2 = 0
        
        return spo2


# ========== 测试 ==========

if __name__ == "__main__":
    rppg = MultispectralRPPG(fps=30, num_channels=9)
    
    # 模拟信号
    np.random.seed(42)
    true_hr = 72
    true_br = 16
    true_spo2 = 98
    
    signals = rppg.simulate_ppg_signal(
        heart_rate=true_hr,
        breathing_rate=true_br,
        spo2=true_spo2,
        duration=10
    )
    
    # 提取参数
    est_hr = rppg.extract_heart_rate(signals)
    est_br = rppg.extract_breathing_rate(signals)
    est_spo2 = rppg.estimate_spo2(signals)
    
    print(f"=== goVISCAN 多光谱 rPPG 模拟结果 ===")
    print(f"  设定心率: {true_hr} bpm | 检测: {est_hr:.1f} bpm | "
          f"误差: {abs(est_hr-true_hr):.1f}")
    print(f"  设定呼吸: {true_br}/min | 检测: {est_br:.1f}/min | "
          f"误差: {abs(est_br-true_br):.1f}")
    print(f"  设定SpO2: {true_spo2}% | 检测: {est_spo2:.1f}% | "
          f"误差: {abs(est_spo2-true_spo2):.1f}")
    print(f"\n  goVISCAN 论文报告误差:")
    print(f"    呼吸: 1.3/min | 心率: 2.8 bpm | SpO2: 3.8%")