ITU-T勧告G.227特性フィルタを通した白色雑音の擬似音声発生器のPythonプログラムとWAV出力

import numpy as np
import scipy.signal as signal
import scipy.io.wavfile as wavfile
import matplotlib.pyplot as plt
import os # ファイルパス操作のため

# --- パラメータ設定 ---
fs = 48000 # Hz (サンプリング周波数) - 一般的なオーディオ用
duration = 10 # seconds (生成する信号の長さ)

# 目標電力のデフォルト設定 (Sメーター調整目的の時の値を維持)
# スプリアス測定では、この絶対レベルよりも基本波に対する相対レベルが重要になりますが、
# 出力されるWAVファイルのRMSレベルの目安として残しておきます。
target_power_dbm = -73 # dBm (例えば、送信機入力レベルを合わせる際の目安に)
R = 50 # Ω (負荷抵抗) - 50Ωシステムを想定

# --- G.227フィルタの設計関数 (共通部分として定義) ---
def design_g227_filter(fs_sampling):
    """
    ITU-T G.227フィルタのデジタル係数を設計する関数。
    Args:
        fs_sampling (int): サンプリング周波数 (Hz)
    Returns:
        tuple: (b, a) デジタルフィルタの分子・分母係数
    """
    # G.227の式をフィルタのゲインとして解釈するための分子・分母係数
    num_g227_corrected = [1, 130, 4001, 36040, 400] # G.227元の式の分母
    den_g227_corrected = [11638, 54050, 91238, 67280, 18400] # G.227元の式の分子

    # p = s / (2 * pi * 1000) という関係に基づき、s領域の係数にスケーリング
    K = 2 * np.pi * 1000 # アナログ周波数スケーリング定数 (rad/s)

    num_s = [
        num_g227_corrected[0] / (K**4),
        num_g227_corrected[1] / (K**3),
        num_g227_corrected[2] / (K**2),
        num_g227_corrected[3] / (K**1),
        num_g227_corrected[4] / (K**0)
    ]
    den_s = [
        den_g227_corrected[0] / (K**4),
        den_g227_corrected[1] / (K**3),
        den_g227_corrected[2] / (K**2),
        den_g227_corrected[3] / (K**1),
        den_g227_corrected[4] / (K**0)
    ]
    b, a = signal.bilinear(num_s, den_s, fs=fs_sampling)
    return b, a

# --- 信号生成関数 ---
def generate_signal(signal_type, fs, duration, b_g227=None, a_g227=None):
    num_samples = int(fs * duration)
    t = np.linspace(0, duration, num_samples, endpoint=False)
    signal_data = np.zeros(num_samples)
    signal_name = ""

    if signal_type == "g227_noise":
        if b_g227 is None or a_g227 is None:
            print("エラー: G.227フィルタ係数が提供されていません。")
            return None, ""
        white_noise = np.random.randn(num_samples)
        signal_data = signal.lfilter(b_g227, a_g227, white_noise)
        signal_name = "ITU-T_G227_Filtered_White_Noise"
    elif signal_type == "single_tone":
        freq_hz = 1000 # Hz (例: 1kHzトーン)
        signal_data = np.sin(2 * np.pi * freq_hz * t)
        signal_name = f"Single_Tone_{freq_hz}Hz"
    elif signal_type == "two_tone":
        freq1_hz = 700  # Hz
        freq2_hz = 1900 # Hz
        # 各トーンの振幅は0.5にし、合成波のピークが1.0を超えないようにする
        signal_data = 0.5 * np.sin(2 * np.pi * freq1_hz * t) + \
                      0.5 * np.sin(2 * np.pi * freq2_hz * t)
        signal_name = f"Two_Tone_{freq1_hz}Hz_{freq2_hz}Hz"
    else:
        print("無効な信号タイプです。'g227_noise', 'single_tone', 'two_tone' のいずれかを指定してください。")
        return None, ""

    return signal_data, signal_name

# --- メイン処理 ---
if __name__ == "__main__":
    # G.227フィルタを事前に設計
    b_g227, a_g227 = design_g227_filter(fs)
    print("G.227フィルタ設計完了。")

    # フィルタの周波数応答をプロット (確認用)
    w_freqz, h_freqz = signal.freqz(b_g227, a_g227, worN=8000, fs=fs)
    plt.figure(figsize=(10, 6))
    plt.plot(w_freqz, 20 * np.log10(abs(h_freqz)))
    plt.title('Calculated G.227 Filter Frequency Response (dB)')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dB]')
    plt.grid(True, which="both", ls="-")
    plt.xscale('log')
    plt.xlim([20, fs/2])
    plt.ylim([-50, 5])
    plt.axvline(x=500, color='r', linestyle='--', label='500 Hz (Min. Loss)')
    plt.legend()
    plt.show()

    # 生成する信号タイプの選択
    print("\n生成する信号タイプを選択してください:")
    print("1: ITU-T G.227フィルタ付き白色雑音 (音声帯域ノイズ)")
    print("2: 単一トーン信号 (例: 1kHz)")
    print("3: 二重トーン信号 (SSB線形性評価用)")
    choice = input("選択 (1/2/3): ")

    signal_type_map = {
        "1": "g227_noise",
        "2": "single_tone",
        "3": "two_tone"
    }
    selected_signal_type = signal_type_map.get(choice)

    if not selected_signal_type:
        print("無効な選択です。プログラムを終了します。")
        exit()

    print(f"\n'{selected_signal_type}' 信号を生成します。")

    # 信号生成
    generated_signal, signal_base_name = generate_signal(selected_signal_type, fs, duration, b_g227, a_g227)

    if generated_signal is None:
        exit()

    # S9 = -73 dBm への調整（目安として）
    target_power_watts = 10**(target_power_dbm / 10) * 1e-3 # W
    target_Vrms = np.sqrt(target_power_watts * R)

    current_Vrms = np.sqrt(np.mean(generated_signal**2))
    print(f"RAW信号のRMS電圧: {current_Vrms:.4e} V")

    # WAVファイルへの出力と正規化 (ピークを-1.0から1.0に)
    # スプリアス測定では、WAVファイル自体の絶対レベルよりも、
    # 送信機への入力レベルと、スペアナで観測される相対レベルが重要です。
    # ここでは、WAVファイルを一般的なオーディオとして適切な音量で保存するため、ピーク正規化を行います。
    max_amplitude = np.max(np.abs(generated_signal))

    if max_amplitude > 1e-18:
        normalized_signal = generated_signal / max_amplitude
    else:
        normalized_signal = np.zeros_like(generated_signal)
        print("警告: 信号の最大振幅がゼロのため、正規化できませんでした。")

    output_filename = f"SSB_Test_Signal_{signal_base_name}_{int(fs/1000)}kHz_{duration}s.wav"
    wavfile.write(output_filename, fs, normalized_signal.astype(np.float32))

    print(f"\nテスト信号を '{output_filename}' として保存しました。")
    print(f"このWAVファイルは {fs} Hz でサンプリングされ、ピーク値が正規化されています。")
    print(f"元のRAW信号のRMS電圧は {current_Vrms:.4e} V です。（これはスケーリング前の値です）")

    print("\n--- SSB送信機スプリアス測定の注意点 ---")
    print("1. 送信機への入力レベル調整:")
    print(f"   生成されたWAVファイルはピーク正規化されています。オーディオインターフェースの出力レベルと、")
    print(f"   必要に応じて外部アッテネータやゲインアンプを調整し、**送信機の規定入力レベル**に合わせる必要があります。")
    print("   例えば、送信機がマイク入力に-30dBmを要求する場合、それに合わせた電圧に調整します。")
    print("2. 測定器の準備:")
    print("   送信機出力とスペクトラムアナライザの間に、**適切な電力のアッテネータやダミーロード、LPFなど**を接続し、")
    print("   スペアナの入力が飽和しないように注意してください。")
    print("3. スプリアスの種類と測定:")
    print("   - **無変調/単一トーン**時: 搬送波漏洩、高調波、広帯域ノイズを確認。")
    print("   - **二重トーン**時: 3次・5次相互変調積 (IMD) のレベルを測定し、メイン信号との差 (IMD3, IMD5など) を評価。")
    print("   - **G.227ノイズ**時: 隣接チャネル漏洩電力 (ACPR) や広帯域ノイズフロアを評価。")
    print("4. 測定周波数範囲:")
    print("   スペアナで、送信機の基本波から離れた広範囲の周波数を掃引し、スプリアスを確認してください。")