Python のまとめ

set MSVC_VERSION=16
set CMAKE_GENERATOR=Visual Studio 16 2019
set PREBUILD_COMMAND=vcvarsall.bat
set PREBUILD_COMMAND_ARGS=x64
# call "vcvarsall.bat" x64 -vcvars_ver=14.11

python %LOCALAPPDATA%\Programs\Python\Python310\Tools\scripts\2to3.py -w .

pyenv shell 3.8.3
python -m pip install jupyterlab jupyter jupyter-console jupytext

sudo apt install python3-spyder
sudo apt install spyder3

import hoge
print(hoge.tax(10))

import numpy as np
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
C = A @ B  # または C = np.dot(A, B)

import numpy as np
from scipy import linalg
from scipy import signal

def matrix_operations_demo():
    """行列計算の基本操作"""
    # 1. 行列の積
    print("1. 行列の積の計算")
    A = np.array([[1, 2], [3, 4]])
    B = np.array([[5, 6], [7, 8]])
    C = A @ B  # または np.dot(A, B)
    print("行列A:\n", A)
    print("行列B:\n", B)
    print("A×B:\n", C)

    # 2. LU分解
    print("\n2. LU分解")
    # 3×3の行列を作成
    A = np.array([[2, -1, 1],
                  [4, 1, -1],
                  [1, 1, 1]])
    P, L, U = linalg.lu(A)
    print("元の行列A:\n", A)
    print("置換行列P:\n", P)
    print("下三角行列L:\n", L)
    print("上三角行列U:\n", U)
    # LU分解の検証
    print("P^T×L×U:\n", P.T @ (L @ U))
    print("元の行列との誤差:", np.max(np.abs(A - P.T @ (L @ U))))

    # LU分解を使用した連立方程式の解法
    b = np.array([1, 2, 3])
    x = linalg.lu_solve(linalg.lu_factor(A), b)
    print("連立方程式Ax = bの解:", x)
    print("残差 ||Ax - b||:", np.linalg.norm(A @ x - b))

    # 3. 逆行列
    print("\n3. 逆行列")
    A = np.array([[4, 7], [2, 6]])
    try:
        A_inv = np.linalg.inv(A)
        print("行列A:\n", A)
        print("逆行列:\n", A_inv)
        # 検証．A×A^(-1) = I
        print("A×A^(-1):\n", A @ A_inv)
        print("逆行列の精度:", np.max(np.abs(A @ A_inv - np.eye(2))))
    except np.linalg.LinAlgError:
        print("逆行列が存在しません")

    # 4. 行列式
    print("\n4. 行列式")
    det_A = np.linalg.det(A)
    print("行列A:\n", A)
    print("行列式:", det_A)

    # 5. 畳み込み積分
    print("\n5. 畳み込み積分")
    # 2D配列（画像を想定）
    image = np.array([[1, 2, 3, 4],
                      [5, 6, 7, 8],
                      [9, 10, 11, 12],
                      [13, 14, 15, 16]])
    # カーネル（エッジ検出フィルタ）
    kernel = np.array([[-1, -1, -1],
                      [-1,  8, -1],
                      [-1, -1, -1]])
    # 異なるモードでの畳み込み演算
    modes = ['full', 'same', 'valid']
    for mode in modes:
        conv_result = signal.convolve2d(image, kernel, mode=mode)
        print(f"\n畳み込み結果（mode='{mode}'）:\n", conv_result)
        print(f"出力サイズ: {conv_result.shape}")

if __name__ == "__main__":
    matrix_operations_demo()

pip install numpy scipy

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

def pca_analysis():
    """主成分分析のデモンストレーション"""

    # サンプルデータの生成
    np.random.seed(42)
    n_samples = 1000
    n_features = 5

    # 相関のあるデータを生成
    cov = np.array([[1, .7, .5, .3, .1],
                    [.7, 1, .4, .2, .1],
                    [.5, .4, 1, .4, .2],
                    [.3, .2, .4, 1, .3],
                    [.1, .1, .2, .3, 1]])
    X = np.random.multivariate_normal(mean=np.zeros(n_features), cov=cov, size=n_samples)

    # データの標準化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # PCAの実行
    pca = PCA()
    X_pca = pca.fit_transform(X_scaled)

    # 結果の表示
    print("各主成分の寄与率:")
    for i, ratio in enumerate(pca.explained_variance_ratio_, 1):
        print(f"第{i}主成分: {ratio:.4f}")

    print("\n累積寄与率:")
    cumsum = np.cumsum(pca.explained_variance_ratio_)
    for i, ratio in enumerate(cumsum, 1):
        print(f"第{i}主成分まで: {ratio:.4f}")

    print("\n主成分負荷量:")
    for i, loadings in enumerate(pca.components_, 1):
        print(f"\n第{i}主成分:")
        for j, loading in enumerate(loadings):
            print(f"変数{j+1}: {loading:.4f}")

    # 可視化
    plt.figure(figsize=(10, 4))

    # スクリープロット
    plt.subplot(121)
    plt.plot(range(1, len(pca.explained_variance_ratio_) + 1),
             pca.explained_variance_ratio_, 'o-')
    plt.xlabel('主成分番号')
    plt.ylabel('寄与率')
    plt.title('スクリープロット')

    # 第1・第2主成分でのデータ分布
    plt.subplot(122)
    plt.scatter(X_pca[:, 0], X_pca[:, 1], alpha=0.5)
    plt.xlabel('第1主成分')
    plt.ylabel('第2主成分')
    plt.title('主成分空間でのデータ分布')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    pca_analysis()

import numpy as np
from sklearn.decomposition import FactorAnalysis
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

def factor_analysis():
    """因子分析のデモンストレーション"""

    # サンプルデータの生成（質問紙調査を想定）
    np.random.seed(42)
    n_samples = 500

    # 2つの潜在因子を仮定
    true_factors = np.random.normal(0, 1, (n_samples, 2))

    # 観測変数を生成（6つの質問項目を想定）
    loadings = np.array([[0.8, 0.2],
                        [0.7, 0.1],
                        [0.6, 0.3],
                        [0.2, 0.8],
                        [0.1, 0.7],
                        [0.3, 0.6]])

    X = true_factors @ loadings.T + np.random.normal(0, 0.1, (n_samples, 6))

    # データの標準化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # 因子分析の実行
    fa = FactorAnalysis(n_components=2, random_state=42)
    X_fa = fa.fit_transform(X_scaled)

    # 結果の表示
    print("因子負荷量:")
    for i, loadings in enumerate(fa.components_.T, 1):
        print(f"\n観測変数{i}:")
        for j, loading in enumerate(loadings, 1):
            print(f"因子{j}: {loading:.4f}")

    # 共通性の計算
    communalities = np.sum(fa.components_**2, axis=0)
    print("\n共通性:")
    for i, comm in enumerate(communalities, 1):
        print(f"観測変数{i}: {comm:.4f}")

    # 因子間相関の計算
    factor_corr = np.corrcoef(X_fa.T)
    print("\n因子間相関:")
    print(factor_corr)

    # 可視化
    plt.figure(figsize=(10, 4))

    # 因子負荷量のプロット
    plt.subplot(121)
    plt.scatter(fa.components_[0], fa.components_[1])
    for i, (x, y) in enumerate(zip(fa.components_[0], fa.components_[1]), 1):
        plt.annotate(f'Var{i}', (x, y))
    plt.xlabel('第1因子')
    plt.ylabel('第2因子')
    plt.title('因子負荷量プロット')
    plt.grid(True)

    # 因子得点の分布
    plt.subplot(122)
    plt.scatter(X_fa[:, 0], X_fa[:, 1], alpha=0.5)
    plt.xlabel('第1因子得点')
    plt.ylabel('第2因子得点')
    plt.title('因子得点の分布')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    factor_analysis()

import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

def lda_analysis():
    """線形判別分析のデモンストレーション"""

    # サンプルデータの生成（3クラス）
    np.random.seed(42)
    n_samples = 300
    n_features = 4

    # 各クラスのデータを生成
    class1 = np.random.multivariate_normal(
        mean=[-2, -2, -2, -2],
        cov=np.eye(4),
        size=n_samples
    )
    class2 = np.random.multivariate_normal(
        mean=[0, 0, 0, 0],
        cov=np.eye(4),
        size=n_samples
    )
    class3 = np.random.multivariate_normal(
        mean=[2, 2, 2, 2],
        cov=np.eye(4),
        size=n_samples
    )

    # データの結合
    X = np.vstack([class1, class2, class3])
    y = np.hstack([np.zeros(n_samples),
                   np.ones(n_samples),
                   np.full(n_samples, 2)])

    # データの標準化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # LDAの実行
    lda = LinearDiscriminantAnalysis()
    X_lda = lda.fit_transform(X_scaled, y)

    # 結果の表示
    print("判別係数:")
    for i, coef in enumerate(lda.coef_, 1):
        print(f"\n判別関数{i}:")
        for j, c in enumerate(coef, 1):
            print(f"変数{j}: {c:.4f}")

    print("\nクラス間正答率:", lda.score(X_scaled, y))

    print("\n各判別関数の寄与率:")
    explained_ratios = lda.explained_variance_ratio_
    for i, ratio in enumerate(explained_ratios, 1):
        print(f"判別関数{i}: {ratio:.4f}")

    # 可視化
    plt.figure(figsize=(12, 4))

    # 第1・第2判別関数でのデータ分布
    plt.subplot(121)
    scatter = plt.scatter(X_lda[:, 0], X_lda[:, 1], c=y, cmap='viridis')
    plt.xlabel('第1判別関数')
    plt.ylabel('第2判別関数')
    plt.title('判別空間でのデータ分布')
    plt.colorbar(scatter)

    # 判別関数の寄与率
    plt.subplot(122)
    plt.bar(range(1, len(explained_ratios) + 1), explained_ratios)
    plt.xlabel('判別関数')
    plt.ylabel('寄与率')
    plt.title('判別関数の寄与率')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    lda_analysis()

import numpy as np
from sklearn.cluster import KMeans, AgglomerativeClustering
from sklearn.preprocessing import StandardScaler
from scipy.cluster.hierarchy import dendrogram, linkage
import matplotlib.pyplot as plt

def clustering_analysis():
    """クラスタリング分析のデモンストレーション"""

    # サンプルデータの生成（3つの自然なクラスタ）
    np.random.seed(42)
    n_samples_per_cluster = 100

    # クラスタ1: 左下のグループ
    cluster1 = np.random.normal(loc=[-2, -2], scale=0.3,
                               size=(n_samples_per_cluster, 2))

    # クラスタ2: 右下のグループ
    cluster2 = np.random.normal(loc=[2, -2], scale=0.3,
                               size=(n_samples_per_cluster, 2))

    # クラスタ3: 上部中央のグループ
    cluster3 = np.random.normal(loc=[0, 2], scale=0.3,
                               size=(n_samples_per_cluster, 2))

    # データの結合
    X = np.vstack([cluster1, cluster2, cluster3])

    # データの標準化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # k-meansクラスタリング
    kmeans = KMeans(n_clusters=3, random_state=42)
    kmeans_labels = kmeans.fit_predict(X_scaled)

    # 階層的クラスタリング
    hc = AgglomerativeClustering(n_clusters=3)
    hc_labels = hc.fit_predict(X_scaled)

    # 結果の表示
    print("K-means クラスタリング結果:")
    print("クラスタ中心:\n", scaler.inverse_transform(kmeans.cluster_centers_))
    print("各クラスタのサンプル数:", np.bincount(kmeans_labels))

    print("\n階層的クラスタリング結果:")
    print("各クラスタのサンプル数:", np.bincount(hc_labels))

    # 可視化
    plt.figure(figsize=(15, 5))

    # オリジナルデータ
    plt.subplot(131)
    X_orig = scaler.inverse_transform(X_scaled)
    plt.scatter(X_orig[:, 0], X_orig[:, 1], c='gray', alpha=0.5)
    plt.title('オリジナルデータ')
    plt.xlabel('Feature 1')
    plt.ylabel('Feature 2')

    # k-meansの結果
    plt.subplot(132)
    scatter = plt.scatter(X_orig[:, 0], X_orig[:, 1],
                         c=kmeans_labels, cmap='viridis')
    centers = scaler.inverse_transform(kmeans.cluster_centers_)
    plt.scatter(centers[:, 0], centers[:, 1], c='red',
                marker='x', s=200, linewidths=3)
    plt.title('K-means クラスタリング結果')
    plt.xlabel('Feature 1')
    plt.ylabel('Feature 2')
    plt.colorbar(scatter)

    # デンドログラム
    plt.subplot(133)
    linkage_matrix = linkage(X_scaled, method='ward')
    dendrogram(linkage_matrix)
    plt.title('階層的クラスタリング デンドログラム')
    plt.xlabel('サンプル番号')
    plt.ylabel('距離')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    clustering_analysis()



5. 特異値分解（SVD）



特異値分解は，行列を3つの行列（U，Σ，V^T）の積に分解する手法である．この分解により，データ圧縮や潜在的意味解析（LSA）などの分析が可能となる．



特異値の大きさは，各成分の重要度を示す指標である．小さい特異値を無視することで，データの次元削減が実現できる．



import numpy as np
import matplotlib.pyplot as plt

def svd_analysis():
    """特異値分解のデモンストレーション"""

    # サンプルデータの生成（ランク3の行列 + ノイズ）
    np.random.seed(42)
    n_rows, n_cols = 100, 50

    # ランク3の行列を生成
    U_true = np.random.normal(0, 1, (n_rows, 3))
    V_true = np.random.normal(0, 1, (n_cols, 3))
    singular_values_true = np.array([10, 5, 2])

    X = U_true @ np.diag(singular_values_true) @ V_true.T

    # ノイズを追加
    X_noisy = X + np.random.normal(0, 0.1, (n_rows, n_cols))

    # SVDの実行
    U, s, Vt = np.linalg.svd(X_noisy, full_matrices=False)

    # 結果の表示
    print("特異値:")
    for i, value in enumerate(s, 1):
        print(f"σ_{i}: {value:.4f}")

    # 低ランク近似の誤差を計算
    errors = []
    ranks = range(1, min(n_rows, n_cols) + 1)

    for r in ranks:
        # ランクrでの近似
        X_approx = U[:, :r] @ np.diag(s[:r]) @ Vt[:r, :]
        # フロベニウスノルムでの誤差
        error = np.linalg.norm(X_noisy - X_approx, 'fro')
        errors.append(error)

    # 可視化
    plt.figure(figsize=(12, 4))

    # 特異値のプロット
    plt.subplot(131)
    plt.plot(range(1, len(s) + 1), s, 'o-')
    plt.yscale('log')
    plt.xlabel('成分番号')
    plt.ylabel('特異値')
    plt.title('特異値のスケープロット')

    # 累積寄与率
    plt.subplot(132)
    cumulative_variance = np.cumsum(s**2) / np.sum(s**2)
    plt.plot(range(1, len(s) + 1), cumulative_variance, 'o-')
    plt.xlabel('成分数')
    plt.ylabel('累積寄与率')
    plt.title('累積寄与率')

    # 近似誤差
    plt.subplot(133)
    plt.plot(ranks, errors, 'o-')
    plt.xlabel('使用する特異値の数')
    plt.ylabel('近似誤差')
    plt.title('低ランク近似の誤差')

    plt.tight_layout()
    plt.show()

    # 画像の再構成例
    plt.figure(figsize=(15, 5))
    ranks_demo = [1, 5, 20]  # デモ用のランク

    plt.subplot(141)
    plt.imshow(X_noisy, cmap='coolwarm')
    plt.title('オリジナルデータ')
    plt.colorbar()

    for i, r in enumerate(ranks_demo, 1):
        X_approx = U[:, :r] @ np.diag(s[:r]) @ Vt[:r, :]
        plt.subplot(1, 4, i+1)
        plt.imshow(X_approx, cmap='coolwarm')
        plt.title(f'ランク{r}での近似')
        plt.colorbar()

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    svd_analysis()




6. QR分解


QR分解は，行列を直交行列Q（Q^T Q = I）と上三角行列Rの積に分解する手法である．この手法は，最小二乗法や固有値計算において重要な役割を持つ．


QR分解は，Gram-Schmidt直交化法を実装したものとして解釈できるという考え方がある．連立方程式の解法や固有値計算のQR法における基礎となる．

import numpy as np
import matplotlib.pyplot as plt
def qr_analysis():
    """QR分解のデモンストレーション"""
    # サンプルデータの生成
    np.random.seed(42)
    n_rows, n_cols = 5, 3  # 長方形行列
    A = np.random.randn(n_rows, n_cols)
    # QR分解の実行（numpyの実装）
    Q, R = np.linalg.qr(A)
    # Gram-Schmidt過程の実装（教育用）
    def gram_schmidt(A):
        m, n = A.shape
        Q = np.zeros((m, n))
        R = np.zeros((n, n))
        for j in range(n):
            v = A[:, j]
            for i in range(j):
                R[i, j] = Q[:, i].T @ A[:, j]
                v = v - R[i, j] * Q[:, i]
            R[j, j] = np.linalg.norm(v)
            Q[:, j] = v / R[j, j]
        return Q, R
    Q_gs, R_gs = gram_schmidt(A)
    # 結果の表示
    print("元の行列 A:")
    print(A)
    print("\nNumPyによるQR分解")
    print("Q行列（直交行列）:")
    print(Q)
    print("\nR行列（上三角行列）:")
    print(R)
    # 検証
    print("\n検証結果:")
    print("Q^T Q ≈ I:")
    print(np.round(Q.T @ Q, decimals=10))
    print("\nQR ≈ A（再構成誤差）:")
    print(np.linalg.norm(A - Q @ R))
    # Gram-Schmidt実装との比較
    print("\nGram-Schmidt実装との差:")
    print("Q の差のノルム:", np.linalg.norm(abs(Q) - abs(Q_gs)))
    print("R の差のノルム:", np.linalg.norm(abs(R) - abs(R_gs)))
    # 可視化
    plt.figure(figsize=(15, 5))
    # 元の行列の可視化
    plt.subplot(131)
    plt.imshow(A, cmap='coolwarm')
    plt.title('元の行列 A')
    plt.colorbar()
    # Q行列の可視化
    plt.subplot(132)
    plt.imshow(Q, cmap='coolwarm')
    plt.title('直交行列 Q')
    plt.colorbar()
    # R行列の可視化
    plt.subplot(133)
    plt.imshow(R, cmap='coolwarm')
    plt.title('上三角行列 R')
    plt.colorbar()
    plt.tight_layout()
    plt.show()
    # Q の列ベクトルの直交性を示す散布図
    if n_cols >= 2:  # 2列以上ある場合のみ
        plt.figure(figsize=(6, 6))
        plt.scatter(Q[:, 0], Q[:, 1], c='b', label='Q の列1 vs 列2')
        # 原点から各ベクトルへの線を描画
        for i in range(2):
            plt.arrow(0, 0, Q[0, i], Q[1, i],
                     head_width=0.05, head_length=0.1,
                     fc='b', ec='b', alpha=0.5)
        plt.axis('equal')
        plt.grid(True)
        plt.title('Q の列ベクトルの直交性')
        plt.xlabel('第1列ベクトル')
        plt.ylabel('第2列ベクトル')
        plt.legend()
        plt.show()
def qr_applications():
    """QR分解の応用例：最小二乗法による線形回帰"""
    # サンプルデータの生成
    np.random.seed(42)
    n_samples = 100
    # 真の関係: y = 2x + 1 + ノイズ
    X = np.random.uniform(0, 10, n_samples)
    y = 2 * X + 1 + np.random.normal(0, 1, n_samples)
    # 設計行列の作成
    A = np.vstack([X, np.ones(n_samples)]).T
    # QR分解を使用した最小二乗法
    Q, R = np.linalg.qr(A)
    beta = np.linalg.solve(R, Q.T @ y)
    # 結果の表示
    print("\n最小二乗法の結果:")
    print(f"推定された傾き: {beta[0]:.4f}")
    print(f"推定された切片: {beta[1]:.4f}")
    # 結果の可視化
    plt.figure(figsize=(8, 6))
    plt.scatter(X, y, alpha=0.5, label='データ点')
    X_sort = np.sort(X)
    plt.plot(X_sort, beta[0] * X_sort + beta[1], 'r-',
             label='QR分解による回帰直線')
    plt.plot(X_sort, 2 * X_sort + 1, 'g--',
             label='真の関係')
    plt.xlabel('X')
    plt.ylabel('y')
    plt.title('QR分解を用いた線形回帰')
    plt.legend()
    plt.grid(True)
    plt.show()
if __name__ == "__main__":
    qr_analysis()
    qr_applications()

このコードでは，以下の機能を提供している：

NumPyによるQR分解の実行と，教育目的のGram-Schmidt過程の実装
分解結果の視覚化と直交性の検証
QR分解の応用例として線形回帰問題の解法を実装



特に以下の点に注意が必要である：

数値的安定性を確保するため，実装にはNumPyのQR分解を使用することが推奨される
直交性の検証では，数値誤差を考慮して適切な閾値を設定する必要がある
応用例では，QR分解を連立方程式の安定的な解法として活用している


コードを実行するには，以下のライブラリが必要である：
pip install numpy matplotlib




音声信号処理の基礎と実装




音声信号処理の基礎と実装について，Python のscipy，librosa，sounddeviceライブラリを用いた処理方法を示している．WAVファイルの入出力，スペクトログラム生成，基本周波数推定，フォルマント分析，環境音解析まで，説明とプロぐr舞うを示している．環境音解析については，スペクトル特徴量，時間領域特徴量，リズム特徴量などの複数の分析手法を提示している．



1. 音声ファイルの入出力



WAVファイルの読み込み，情報表示，書き出しを行う．scipyのwavfileモジュールやlibrosaを使用する．



サンプリング周波数やビット深度などの音声パラメータも合わせて行う．



import numpy as np
from scipy.io import wavfile
import librosa
import sounddevice as sd
import matplotlib.pyplot as plt

def audio_file_handling():
    """音声ファイルの入出力と情報表示のデモンストレーション"""

    # WAVファイルの読み込み（scipyを使用）
    def read_wav_scipy(filename):
        """scipyによるWAVファイル読み込み"""
        sampling_rate, data = wavfile.read(filename)
        print("=== Scipyでの読み込み結果 ===")
        print(f"サンプリング周波数: {sampling_rate} Hz")
        print(f"チャンネル数: {data.shape[1] if len(data.shape) > 1 else 1}")
        print(f"サンプル数: {len(data)}")
        print(f"時間長: {len(data)/sampling_rate:.2f} 秒")
        return sampling_rate, data

    # WAVファイルの読み込み（librosaを使用）
    def read_wav_librosa(filename):
        """librosaによるWAVファイル読み込み"""
        data, sampling_rate = librosa.load(filename, sr=None)
        print("\n=== Librosaでの読み込み結果 ===")
        print(f"サンプリング周波数: {sampling_rate} Hz")
        print(f"チャンネル数: {data.shape[1] if len(data.shape) > 1 else 1}")
        print(f"サンプル数: {len(data)}")
        print(f"時間長: {len(data)/sampling_rate:.2f} 秒")
        return sampling_rate, data

    # WAVファイルの書き出し
    def write_wav(filename, sampling_rate, data):
        """WAVファイルの書き出し"""
        normalized_data = np.clip(data, -1.0, 1.0)
        int_data = (normalized_data * 32767).astype(np.int16)
        wavfile.write(filename, sampling_rate, int_data)

    # 音声の再生
    def play_audio(data, sampling_rate):
        """音声データの再生"""
        sd.play(data, sampling_rate)
        sd.wait()  # 再生終了まで待機

    # テスト用の音声信号生成（1秒のサイン波）
    def generate_test_signal():
        """テスト用の音声信号生成"""
        duration = 1.0  # 1秒
        sampling_rate = 44100  # 44.1kHz
        t = np.linspace(0, duration, int(sampling_rate * duration))
        signal = 0.5 * np.sin(2 * np.pi * 440 * t)
        return sampling_rate, signal

    # メイン処理
    try:
        # テスト信号の生成と保存
        sampling_rate, test_signal = generate_test_signal()
        write_wav("test_signal.wav", sampling_rate, test_signal)

        # 保存した音声ファイルの読み込みと情報表示
        sr_scipy, data_scipy = read_wav_scipy("test_signal.wav")
        sr_librosa, data_librosa = read_wav_librosa("test_signal.wav")

        # 波形プロット
        plt.figure(figsize=(12, 4))
        time = np.arange(len(test_signal)) / sampling_rate
        plt.plot(time, test_signal)
        plt.xlabel('時間 [秒]')
        plt.ylabel('振幅')
        plt.title('テスト信号の波形')
        plt.grid(True)
        plt.show()

        # 音声の再生
        print("\n音声を再生します...")
        play_audio(test_signal, sampling_rate)

    except Exception as e:
        print(f"エラーが発生しました: {e}")

if __name__ == "__main__":
    audio_file_handling()





2. スペクトログラムの生成と表示



スペクトログラムは，音声信号の時間-周波数表現である．短時間フーリエ変換（STFT）を用いて算出する．



窓関数の選択やフレーム長，オーバーラップ率などのパラメータが結果のスペクトログラムに影響を与える．



import numpy as np
from scipy import signal
import librosa
import librosa.display
import matplotlib.pyplot as plt

def spectrogram_analysis():
    """スペクトログラム分析のデモンストレーション"""

    # テスト信号の生成（周波数が時間とともに変化する信号）
    def generate_chirp_signal():
        """チャープ信号の生成"""
        duration = 3.0  # 3秒
        sampling_rate = 44100  # 44.1kHz
        t = np.linspace(0, duration, int(sampling_rate * duration))
        # 周波数が100Hzから4000Hzまで指数的に変化
        chirp = signal.chirp(t, f0=100, f1=4000, t1=duration, method='exponential')
        return sampling_rate, chirp

    # STFTによるスペクトログラム計算（scipy使用）
    def calculate_spectrogram_scipy(audio_data, sampling_rate):
        """Scipyを使用したスペクトログラム計算"""
        nperseg = 2048  # フレーム長
        noverlap = nperseg // 2  # オーバーラップ

        frequencies, times, Sxx = signal.spectrogram(
            audio_data,
            fs=sampling_rate,
            window='hann',
            nperseg=nperseg,
            noverlap=noverlap,
            scaling='spectrum'
        )
        return frequencies, times, 10 * np.log10(Sxx + 1e-10)

    # STFTによるスペクトログラム計算（librosa使用）
    def calculate_spectrogram_librosa(audio_data, sampling_rate):
        """Librosaを使用したスペクトログラム計算"""
        n_fft = 2048  # FFTサイズ
        hop_length = n_fft // 2  # ホップ長

        D = librosa.stft(audio_data, n_fft=n_fft, hop_length=hop_length,
                        window='hann', center=True)
        S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max)
        return S_db

    # メイン処理
    sampling_rate, chirp = generate_chirp_signal()

    # Scipyによるスペクトログラム
    freqs, times, spec_scipy = calculate_spectrogram_scipy(chirp, sampling_rate)

    # Librosaによるスペクトログラム
    spec_librosa = calculate_spectrogram_librosa(chirp, sampling_rate)

    # 結果の可視化
    plt.figure(figsize=(15, 10))

    # 波形
    plt.subplot(3, 1, 1)
    time = np.arange(len(chirp)) / sampling_rate
    plt.plot(time, chirp)
    plt.xlabel('時間 [秒]')
    plt.ylabel('振幅')
    plt.title('チャープ信号の波形')
    plt.grid(True)

    # Scipyによるスペクトログラム
    plt.subplot(3, 1, 2)
    plt.pcolormesh(times, freqs, spec_scipy, shading='gouraud')
    plt.ylabel('周波数 [Hz]')
    plt.xlabel('時間 [秒]')
    plt.title('Scipyによるスペクトログラム')
    plt.colorbar(label='パワー [dB]')

    # Librosaによるスペクトログラム
    plt.subplot(3, 1, 3)
    librosa.display.specshow(spec_librosa, sr=sampling_rate, x_axis='time',
                           y_axis='hz', hop_length=1024)
    plt.ylabel('周波数 [Hz]')
    plt.xlabel('時間 [秒]')
    plt.title('Librosaによるスペクトログラム')
    plt.colorbar(label='振幅 [dB]')

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    spectrogram_analysis()





3. 基本周波数（F0）推定



基本周波数は音声の高さを表す重要なパラメータである．自己相関法やケプストラム法など，様々な推定手法が存在する．



有声/無声の判定と，オクターブ誤りの処理が重要な課題である．



import numpy as np
import librosa
import matplotlib.pyplot as plt
from scipy.signal import correlate
from scipy.io import wavfile

def f0_estimation_demo():
    """基本周波数推定のデモンストレーション"""

    def generate_test_signal(f0=440, duration=1.0, sr=44100):
        """基本周波数を持つテスト信号の生成"""
        t = np.linspace(0, duration, int(sr * duration))
        # 基本波と倍音を含む信号
        signal = 0.5 * np.sin(2 * np.pi * f0 * t) + \
                0.25 * np.sin(2 * np.pi * 2 * f0 * t) + \
                0.125 * np.sin(2 * np.pi * 3 * f0 * t)
        return signal, sr

    def autocorrelation_f0(signal, sr, frame_length=2048, hop_length=512):
        """自己相関法によるF0推定"""
        f0_values = []
        times = []

        for i in range(0, len(signal) - frame_length, hop_length):
            frame = signal[i:i + frame_length]
            # 自己相関の計算
            corr = correlate(frame, frame, mode='full')
            corr = corr[len(corr)//2:]

            # 最初のピークを探す（直流成分を除く）
            peak_idx = np.argmax(corr[50:]) + 50
            if corr[peak_idx] > 0.1:  # 有声音の閾値
                f0 = sr / peak_idx
                if 50 <= f0 <= 2000:  # 有効なF0範囲
                    f0_values.append(f0)
                else:
                    f0_values.append(0)
            else:
                f0_values.append(0)

            times.append(i / sr)

        return np.array(f0_values), np.array(times)

    def librosa_f0(signal, sr):
        """librosaを使用したF0推定"""
        f0, voiced_flag, voiced_probs = librosa.pyin(
            signal,
            fmin=librosa.note_to_hz('C2'),
            fmax=librosa.note_to_hz('C7'),
            sr=sr
        )
        return f0, voiced_flag

    # テスト信号の生成と分析
    signal, sr = generate_test_signal(f0=440)

    # 自己相関法によるF0推定
    f0_autocorr, times_autocorr = autocorrelation_f0(signal, sr)

    # librosaによるF0推定
    f0_librosa, voiced_flag = librosa_f0(signal, sr)
    times_librosa = librosa.times_like(f0_librosa, sr=sr)

    # 結果の可視化
    plt.figure(figsize=(12, 8))

    # 波形
    plt.subplot(3, 1, 1)
    t = np.linspace(0, len(signal)/sr, len(signal))
    plt.plot(t, signal)
    plt.xlabel('時間 [秒]')
    plt.ylabel('振幅')
    plt.title('テスト信号の波形')
    plt.grid(True)

    # 自己相関法によるF0
    plt.subplot(3, 1, 2)
    plt.plot(times_autocorr, f0_autocorr, 'o-', label='推定F0')
    plt.axhline(y=440, color='r', linestyle='--', label='真のF0')
    plt.xlabel('時間 [秒]')
    plt.ylabel('周波数 [Hz]')
    plt.title('自己相関法によるF0推定結果')
    plt.legend()
    plt.grid(True)

    # librosaによるF0
    plt.subplot(3, 1, 3)
    plt.plot(times_librosa[voiced_flag], f0_librosa[voiced_flag],
             'o-', label='推定F0（有声区間）')
    plt.axhline(y=440, color='r', linestyle='--', label='真のF0')
    plt.xlabel('時間 [秒]')
    plt.ylabel('周波数 [Hz]')
    plt.title('PYIN（librosa）によるF0推定結果')
    plt.legend()
    plt.grid(True)

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    f0_estimation_demo()





4. フォルマント分析



フォルマントは，音声スペクトルの共鳴周波数であり，母音の特徴を表す重要なパラメータである．推定法としては，線形予測分析（LPC）がある．


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

def formant_analysis():
    """フォルマント分析のデモンストレーション"""

    def generate_vowel(formants, duration=1.0, sr=16000):
        """フォルマントを持つ母音の生成"""
        t = np.linspace(0, duration, int(sr * duration))
        f0 = 130  # 基本周波数

        # 声帯波源（のこぎり波）
        source = signal.sawtooth(2 * np.pi * f0 * t)

        # フォルマントフィルタの作成
        vowel = np.zeros_like(t)
        for formant in formants:
            bw = formant * 0.1  # バンド幅
            w0 = 2 * np.pi * formant / sr
            r = np.exp(-np.pi * bw / sr)

            # 共鳴フィルタ
            a = [1, -2*r*np.cos(w0), r**2]
            b = [1]
            vowel += signal.lfilter(b, a, source)

        return vowel / np.max(np.abs(vowel)), sr

    def estimate_formants(signal, sr, order=12):
        """LPCによるフォルマント推定"""
        # LPC係数の計算
        a = librosa.lpc(signal, order=order)

        # 根を求める
        roots = np.roots(a)
        roots = roots[np.imag(roots) >= 0]  # 正の虚部を持つ根のみ

        # 角周波数から周波数へ変換
        angles = np.arctan2(np.imag(roots), np.real(roots))
        freqs = angles * (sr / (2 * np.pi))

        # バンド幅の計算
        bandwidth = -1/2 * (sr/(2*np.pi)) * np.log(np.abs(roots))

        # 周波数でソート
        sorted_idx = np.argsort(freqs)
        freqs = freqs[sorted_idx]
        bandwidth = bandwidth[sorted_idx]

        return freqs, bandwidth

    # 母音/a/のフォルマント
    formants_a = [730, 1090, 2440]
    # 母音/i/のフォルマント
    formants_i = [270, 2290, 3010]

    # 母音の生成
    vowel_a, sr = generate_vowel(formants_a)
    vowel_i, _ = generate_vowel(formants_i)

    # フォルマント推定
    freqs_a, bw_a = estimate_formants(vowel_a, sr)
    freqs_i, bw_i = estimate_formants(vowel_i, sr)

    # 結果の可視化
    plt.figure(figsize=(15, 10))

    # 母音/a/の波形とスペクトル
    plt.subplot(2, 2, 1)
    plt.plot(vowel_a)
    plt.title('母音/a/の波形')
    plt.xlabel('サンプル')
    plt.ylabel('振幅')

    plt.subplot(2, 2, 2)
    freq, spec = signal.welch(vowel_a, sr, nperseg=2048)
    plt.plot(freq, 10 * np.log10(spec))
    plt.vlines(freqs_a[:3], plt.ylim()[0], plt.ylim()[1],
              colors='r', linestyles='--', label='推定フォルマント')
    plt.vlines(formants_a, plt.ylim()[0], plt.ylim()[1],
              colors='g', linestyles=':', label='真のフォルマント')
    plt.title('母音/a/のスペクトル')
    plt.xlabel('周波数 [Hz]')
    plt.ylabel('パワー [dB]')
    plt.legend()

    # 母音/i/の波形とスペクトル
    plt.subplot(2, 2, 3)
    plt.plot(vowel_i)
    plt.title('母音/i/の波形')
    plt.xlabel('サンプル')
    plt.ylabel('振幅')

    plt.subplot(2, 2, 4)
    freq, spec = signal.welch(vowel_i, sr, nperseg=2048)
    plt.plot(freq, 10 * np.log10(spec))
    plt.vlines(freqs_i[:3], plt.ylim()[0], plt.ylim()[1],
              colors='r', linestyles='--', label='推定フォルマント')
    plt.vlines(formants_i, plt.ylim()[0], plt.ylim()[1],
              colors='g', linestyles=':', label='真のフォルマント')
    plt.title('母音/i/のスペクトル')
    plt.xlabel('周波数 [Hz]')
    plt.ylabel('パワー [dB]')
    plt.legend()

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    formant_analysis()





5. 環境音の解析



環境音の解析では，時間-周波数表現に加えて，音響的特徴量の抽出が重要である．特に，以下の特徴量が有用である．



環境音は非定常な性質を持つため，統計的な特徴量や変動パターンの分析が重要である．


import numpy as np
import librosa
import librosa.display
import scipy.stats
from scipy.signal import welch
import matplotlib.pyplot as plt
import warnings

def environmental_sound_analysis():
    """環境音の解析デモンストレーション"""

    def compute_spectral_features(y, sr, n_fft=2048, hop_length=512):
        """スペクトル特徴量の計算（エラー処理追加）"""
        # 無音区間の検出
        if np.all(y == 0):
            return {
                'centroid': np.zeros(1),
                'bandwidth': np.zeros(1),
                'rolloff': np.zeros(1),
                'flatness': np.zeros(1)
            }

        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            # スペクトル重心
            spectral_centroid = librosa.feature.spectral_centroid(
                y=y, sr=sr, n_fft=n_fft, hop_length=hop_length
            )[0]

            # スペクトル帯域幅
            spectral_bandwidth = librosa.feature.spectral_bandwidth(
                y=y, sr=sr, n_fft=n_fft, hop_length=hop_length
            )[0]

            # スペクトルロールオフ
            spectral_rolloff = librosa.feature.spectral_rolloff(
                y=y, sr=sr, n_fft=n_fft, hop_length=hop_length
            )[0]

            # スペクトル平坦度
            spectral_flatness = librosa.feature.spectral_flatness(
                y=y, n_fft=n_fft, hop_length=hop_length
            )[0]

        return {
            'centroid': spectral_centroid,
            'bandwidth': spectral_bandwidth,
            'rolloff': spectral_rolloff,
            'flatness': spectral_flatness
        }

    def compute_temporal_features(y, sr, frame_length=2048, hop_length=512):
        """時間領域特徴量の計算（エラー処理追加）"""
        # 無音区間の検出
        if np.all(y == 0):
            return {
                'rms': np.zeros(1),
                'zcr': np.zeros(1),
                'crest_factor': 0,
                'kurtosis': 0
            }

        # RMSエネルギー
        rms = librosa.feature.rms(
            y=y, frame_length=frame_length, hop_length=hop_length
        )[0]

        # ゼロ交差率
        zero_crossing_rate = librosa.feature.zero_crossing_rate(
            y=y, frame_length=frame_length, hop_length=hop_length
        )[0]

        # 波形の統計量（ゼロ割り防止）
        rms_val = np.sqrt(np.mean(y**2)) + 1e-10
        crest_factor = np.max(np.abs(y)) / rms_val
        kurtosis = scipy.stats.kurtosis(y, fisher=True)

        return {
            'rms': rms,
            'zcr': zero_crossing_rate,
            'crest_factor': crest_factor,
            'kurtosis': kurtosis
        }

    def compute_rhythm_features(y, sr):
        """リズム特徴量の計算"""
        # テンポグラム
        onset_env = librosa.onset.onset_strength(y=y, sr=sr)
        tempogram = librosa.feature.tempogram(
            onset_envelope=onset_env,
            sr=sr
        )

        # オンセット検出
        onset_frames = librosa.onset.onset_detect(
            onset_envelope=onset_env,
            sr=sr
        )

        return {
            'tempogram': tempogram,
            'onset_frames': onset_frames
        }

    def analyze_noise_color(y, sr):
        """ノイズの色分析（パワースペクトル密度の傾き）"""
        freqs, psd = welch(y, sr, nperseg=2048)
        # 対数変換
        log_freqs = np.log10(freqs[1:])  # 0Hzを除外
        log_psd = np.log10(psd[1:])
        # 線形回帰で傾きを計算
        slope, _, _, _, _ = scipy.stats.linregress(log_freqs, log_psd)
        return slope, freqs[1:], psd[1:]

    # テスト信号として雨音をシミュレート
def simulate_rain(duration=5.0, sr=44100):
        """より現実的な雨音のシミュレーション"""
        t = np.linspace(0, duration, int(sr * duration))
        rain = np.zeros_like(t)

        # 雨滴の種類（大小）を考慮
        drop_types = {
            'small': {
                'freq_range': (2000, 4000),
                'duration_range': (0.01, 0.03),
                'amplitude_range': (0.1, 0.3)
            },
            'medium': {
                'freq_range': (1000, 2000),
                'duration_range': (0.03, 0.06),
                'amplitude_range': (0.3, 0.6)
            },
            'large': {
                'freq_range': (500, 1000),
                'duration_range': (0.06, 0.1),
                'amplitude_range': (0.6, 1.0)
            }
        }

        # 各種類の雨滴を生成
        for drop_type, params in drop_types.items():
            # 雨滴の数（ポアソン分布に従う）
            n_drops = np.random.poisson(duration *
                     (100 if drop_type == 'small' else
                      50 if drop_type == 'medium' else 20))

            for _ in range(n_drops):
                # 雨滴の時刻
                t_drop = np.random.random() * duration
                idx = int(t_drop * sr)

                # 雨滴のパラメータ
                freq = np.random.uniform(*params['freq_range'])
                drop_duration = np.random.uniform(*params['duration_range'])
                amplitude = np.random.uniform(*params['amplitude_range'])

                # 雨滴の波形生成
                n_samples = int(drop_duration * sr)
                if idx + n_samples >= len(rain):
                    continue

                # 減衰正弦波として雨滴を生成
                t_drop = np.arange(n_samples) / sr
                envelope = np.exp(-t_drop * 30)  # 指数減衰
                drop = amplitude * envelope * np.sin(2 * np.pi * freq * t_drop)

                # ランダムなフィルタリングで音色を変化
                drop = scipy.signal.lfilter(
                    [1], [1, -0.99], drop
                )

                rain[idx:idx+n_samples] += drop

        # クリッピングを防ぎつつ正規化
        rain = rain / (np.max(np.abs(rain)) + 1e-10)

        return rain, sr

    # 解析の実行
    y, sr = simulate_rain()

    # 各種特徴量の計算
    spectral_features = compute_spectral_features(y, sr)
    temporal_features = compute_temporal_features(y, sr)
    rhythm_features = compute_rhythm_features(y, sr)
    noise_slope, freqs_psd, psd = analyze_noise_color(y, sr)

    # 結果の可視化
    plt.figure(figsize=(15, 10))

    # 波形
    plt.subplot(3, 2, 1)
    plt.plot(np.linspace(0, len(y)/sr, len(y)), y)
    plt.title('雨音波形')
    plt.xlabel('時間 [秒]')
    plt.ylabel('振幅')

    # スペクトログラム
    plt.subplot(3, 2, 2)
    D = librosa.amplitude_to_db(np.abs(librosa.stft(y)), ref=np.max)
    librosa.display.specshow(D, sr=sr, x_axis='time', y_axis='log')
    plt.colorbar(format='%+2.0f dB')
    plt.title('スペクトログラム'
    # スペクトル特徴量の時間変化
    plt.subplot(3, 2, 3)
    times = np.linspace(0, len(y)/sr, len(spectral_features['centroid']))
    plt.plot(times, librosa.amplitude_to_db(spectral_features['flatness']),
             label='平坦度')
    plt.plot(times, spectral_features['centroid']/sr, label='重心')
    plt.legend()
    plt.title('スペクトル特徴量の時間変化')
    plt.xlabel('時間 [秒]')

    # パワースペクトル密度（ノイズの色分析）
    plt.subplot(3, 2, 4)
    plt.loglog(freqs_psd, psd)
    plt.title(f'パワースペクトル密度 (傾き: {noise_slope:.2f})')
    plt.xlabel('周波数 [Hz]')
    plt.ylabel('PSD')

    # テンポグラム
    plt.subplot(3, 2, 5)
    librosa.display.specshow(rhythm_features['tempogram'],
                           sr=sr, x_axis='time')
    plt.title('テンポグラム')
    plt.xlabel('時間 [秒]')
    plt.ylabel('テンポ [BPM]')

    # エネルギーとゼロ交差率
    plt.subplot(3, 2, 6)
    times_rms = np.linspace(0, len(y)/sr, len(temporal_features['rms']))
    plt.plot(times_rms, temporal_features['rms'], label='RMS')
    plt.plot(times_rms, temporal_features['zcr'], label='ZCR')
    plt.legend()
    plt.title('時間領域特徴量')
    plt.xlabel('時間 [秒]')

    plt.tight_layout()
    plt.show()

    # 統計値の表示
    print(f"クレストファクター: {temporal_features['crest_factor']:.2f}")
    print(f"尖度: {temporal_features['kurtosis']:.2f}")
    print(f"ノイズの色（スペクトル傾き）: {noise_slope:.2f}")

if __name__ == "__main__":
    environmental_sound_analysis()



このコードは環境音分析の主要である以下の特徴量を計算している：


スペクトル特徴量（重心，帯域幅，ロールオフ，平坦度）
時間領域特徴量（RMS，ゼロ交差率，クレストファクター，尖度）
リズム特徴量（テンポグラム，オンセット検出）
ノイズの色分析（パワースペクトル密度の傾き）




特に以下の点が重要である：

非定常性の定量化
時間-周波数表現の適切な選択
統計的特徴量の組み合わせ