import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
fs = 100000         # サンプリング周波数 [Hz]
T = 0.1             # 信号長 [s]
t = np.arange(0, T, 1/fs)
fc = 10000          # 搬送波周波数 [Hz]
fm = 500            # 信号周波数 [Hz]
Ac = 1.0            # 搬送波振幅

# ベースバンド信号（正弦波）
x = np.cos(2 * np.pi * fm * t)

# 搬送波
carrier = np.cos(2 * np.pi * fc * t)

# 各変調度でのAM信号
m_list = [0.3, 0.7, 1.0, 1.5]

fig, axes = plt.subplots(len(m_list), 2, figsize=(16, 14))

for i, m in enumerate(m_list):
    # AM変調信号
    s_am = Ac * (1 + m * x) * carrier

    # 包絡線
    envelope_upper = Ac * (1 + m * x)
    envelope_lower = -Ac * (1 + m * x)

    # 時間波形
    axes[i, 0].plot(t * 1000, s_am, 'b-', linewidth=0.5, alpha=0.7)
    axes[i, 0].plot(t * 1000, envelope_upper, 'r--', linewidth=1.5,
                    label='Envelope')
    axes[i, 0].plot(t * 1000, envelope_lower, 'r--', linewidth=1.5)
    axes[i, 0].set_xlabel('Time [ms]')
    axes[i, 0].set_ylabel('Amplitude')
    axes[i, 0].set_title(f'AM Signal (m = {m})')
    axes[i, 0].set_xlim(0, 10)
    axes[i, 0].legend(loc='upper right')
    axes[i, 0].grid(True, alpha=0.3)

    if m > 1.0:
        axes[i, 0].set_title(f'AM Signal (m = {m}, Over-modulation!)',
                              color='red')

    # スペクトル
    N_fft = len(t)
    freq = np.fft.fftfreq(N_fft, 1/fs)
    S = np.abs(np.fft.fft(s_am)) / N_fft
    pos_mask = freq >= 0

    axes[i, 1].plot(freq[pos_mask] / 1000, 2 * S[pos_mask], 'b-',
                    linewidth=1)
    axes[i, 1].set_xlabel('Frequency [kHz]')
    axes[i, 1].set_ylabel('Magnitude')
    axes[i, 1].set_title(f'Spectrum (m = {m})')
    axes[i, 1].set_xlim(8, 12)
    axes[i, 1].grid(True, alpha=0.3)

    # 搬送波・USB・LSBの位置を注釈
    axes[i, 1].axvline(x=fc/1000, color='r', linestyle=':', alpha=0.5)
    axes[i, 1].axvline(x=(fc+fm)/1000, color='g', linestyle=':',
                        alpha=0.5)
    axes[i, 1].axvline(x=(fc-fm)/1000, color='g', linestyle=':',
                        alpha=0.5)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 変調度と電力効率の関係
m_range = np.linspace(0, 1, 100)
eta = m_range**2 / (2 + m_range**2)

fig, ax = plt.subplots(figsize=(10, 6))

ax.plot(m_range, eta * 100, 'b-', linewidth=2)
ax.set_xlabel('Modulation Index m')
ax.set_ylabel('Power Efficiency η [%]')
ax.set_title('AM Power Efficiency vs Modulation Index')
ax.grid(True, alpha=0.3)

# 特定の点を注釈
for m_val in [0.3, 0.5, 0.7, 1.0]:
    eta_val = m_val**2 / (2 + m_val**2) * 100
    ax.plot(m_val, eta_val, 'ro', markersize=8)
    ax.annotate(f'm={m_val}: η={eta_val:.1f}%',
                xy=(m_val, eta_val),
                xytext=(m_val + 0.05, eta_val + 2),
                fontsize=10, color='red')

ax.set_xlim(0, 1.05)
ax.set_ylim(0, 40)
plt.tight_layout()
plt.show()

# 数値確認
print("=== AM変調の電力効率 ===")
for m_val in [0.3, 0.5, 0.7, 1.0]:
    eta_val = m_val**2 / (2 + m_val**2) * 100
    print(f"  m = {m_val:.1f}: η = {eta_val:.2f}%")

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
fs = 200000
T = 0.05
t = np.arange(0, T, 1/fs)
fc = 20000
Ac = 1.0
m = 0.8

# 複合ベースバンド信号
x = 0.6 * np.cos(2 * np.pi * 300 * t) + 0.4 * np.cos(2 * np.pi * 800 * t)
x = x / np.max(np.abs(x))  # 正規化

# AM変調
s_am = Ac * (1 + m * x) * np.cos(2 * np.pi * fc * t)

# 雑音付加（SNR = 20 dB）
snr_dB = 20
signal_power = np.mean(s_am**2)
noise_power = signal_power / (10**(snr_dB/10))
noise = np.sqrt(noise_power) * np.random.randn(len(t))
s_noisy = s_am + noise

# --- 包絡線検波（ヒルベルト変換を使用） ---
from scipy.signal import hilbert
analytic_signal = hilbert(s_noisy)
envelope = np.abs(analytic_signal)

# DCブロック（ハイパスフィルタ代用: 平均値を引く）
envelope_ac = envelope - np.mean(envelope)

# --- 同期検波 ---
# 局部発振器
lo = np.cos(2 * np.pi * fc * t)
demod_sync = s_noisy * lo * 2  # 2倍は振幅補正

# ローパスフィルタ
from scipy.signal import butter, filtfilt
b, a = butter(5, 2000 / (fs/2), btype='low')
demod_sync_filtered = filtfilt(b, a, demod_sync)
demod_sync_ac = demod_sync_filtered - np.mean(demod_sync_filtered)

# 可視化
fig, axes = plt.subplots(4, 1, figsize=(14, 14))

# 元の信号
axes[0].plot(t * 1000, x, 'b-', linewidth=1.5)
axes[0].set_xlabel('Time [ms]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Original Baseband Signal x(t)')
axes[0].set_xlim(0, 20)
axes[0].grid(True, alpha=0.3)

# AM変調信号（雑音付き）
axes[1].plot(t * 1000, s_noisy, 'b-', linewidth=0.3, alpha=0.5)
axes[1].plot(t * 1000, Ac * (1 + m * x), 'r-', linewidth=1.5,
             label='True Envelope')
axes[1].set_xlabel('Time [ms]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title(f'Noisy AM Signal (SNR = {snr_dB} dB)')
axes[1].set_xlim(0, 20)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# 包絡線検波結果
scale_env = m * Ac
axes[2].plot(t * 1000, m * x, 'b-', linewidth=1.5, label='Original (scaled)')
axes[2].plot(t * 1000, envelope_ac / Ac, 'r-', linewidth=1.5, alpha=0.7,
             label='Envelope Detection')
axes[2].set_xlabel('Time [ms]')
axes[2].set_ylabel('Amplitude')
axes[2].set_title('Envelope Detection Result')
axes[2].set_xlim(0, 20)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

# 同期検波結果
axes[3].plot(t * 1000, m * x, 'b-', linewidth=1.5, label='Original (scaled)')
axes[3].plot(t * 1000, demod_sync_ac / Ac, 'g-', linewidth=1.5, alpha=0.7,
             label='Synchronous Detection')
axes[3].set_xlabel('Time [ms]')
axes[3].set_ylabel('Amplitude')
axes[3].set_title('Synchronous Detection Result')
axes[3].set_xlim(0, 20)
axes[3].legend()
axes[3].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert

# パラメータ
fs = 100000
T = 0.05
t = np.arange(0, T, 1/fs)
fc = 10000
Ac = 1.0

# ベースバンド信号
fm1, fm2 = 400, 900
x = 0.7 * np.cos(2 * np.pi * fm1 * t) + 0.3 * np.cos(2 * np.pi * fm2 * t)

# AM (DSB-FC)
m = 0.8
s_am = Ac * (1 + m * x) * np.cos(2 * np.pi * fc * t)

# DSB-SC
s_dsb = Ac * x * np.cos(2 * np.pi * fc * t)

# SSB (USB) — ヒルベルト変換を使用
x_hat = np.imag(hilbert(x))  # x(t)のヒルベルト変換
s_ssb = 0.5 * Ac * (x * np.cos(2 * np.pi * fc * t)
                     - x_hat * np.sin(2 * np.pi * fc * t))

# スペクトル計算
def compute_spectrum(signal, fs):
    N_fft = len(signal)
    freq = np.fft.fftfreq(N_fft, 1/fs)
    S = np.abs(np.fft.fft(signal)) / N_fft
    pos = freq >= 0
    return freq[pos], 2 * S[pos]

fig, axes = plt.subplots(3, 2, figsize=(16, 12))

signals = [('AM (DSB-FC)', s_am), ('DSB-SC', s_dsb), ('SSB (USB)', s_ssb)]

for i, (name, s) in enumerate(signals):
    # 時間波形
    axes[i, 0].plot(t * 1000, s, 'b-', linewidth=0.5)
    axes[i, 0].set_xlabel('Time [ms]')
    axes[i, 0].set_ylabel('Amplitude')
    axes[i, 0].set_title(f'{name} - Time Domain')
    axes[i, 0].set_xlim(0, 10)
    axes[i, 0].grid(True, alpha=0.3)

    # スペクトル
    f, S_mag = compute_spectrum(s, fs)
    axes[i, 1].plot(f / 1000, S_mag, 'b-', linewidth=1)
    axes[i, 1].set_xlabel('Frequency [kHz]')
    axes[i, 1].set_ylabel('Magnitude')
    axes[i, 1].set_title(f'{name} - Spectrum')
    axes[i, 1].set_xlim(8, 12)
    axes[i, 1].grid(True, alpha=0.3)

    # 搬送波位置
    axes[i, 1].axvline(x=fc/1000, color='r', linestyle=':', alpha=0.5)

plt.tight_layout()
plt.show()

# 電力比較
powers = {}
for name, s in signals:
    powers[name] = np.mean(s**2)

print("=== 各方式の平均電力 ===")
for name, p in powers.items():
    print(f"  {name}: {p:.4f}")