TSP信号とインパルス応答の取得

このドキュメントでは、TSP（Time Stretched Pulse）信号の生成方法と、インパルス応答を取得する方法について説明します。Pythonを使用して、TSP信号を生成し、その信号を逆TSP信号と畳み込むことでインパルス応答を取得します。また、インパルス応答を窓関数で短縮する方法についても説明します。

TSP信号の生成

以下のコードは、TSP信号を生成する関数と逆TSP信号を生成する関数を定義しています。

import numpy as np
import scipy.signal as signal

def generate_tsp(N):
    tsp = signal.chirp(np.linspace(0, 1, N), f0=0.1, f1=0.9, t1=N, method='linear')
    return tsp

def generate_inverse_tsp(tsp):
    N = len(tsp)
    inverse_tsp = np.flip(tsp)
    inverse_tsp = inverse_tsp / np.sum(tsp ** 2)
    return inverse_tsp

# TSP信号の生成
N = 1024
tsp = generate_tsp(N)
inverse_tsp = generate_inverse_tsp(tsp)

インパルス応答の取得

以下のコードは、TSP信号を使用して録音した信号と逆TSP信号を畳み込むことでインパルス応答を取得します。また、取得したインパルス応答を窓関数で短縮し、FIRフィルタのパラメータとして利用する方法も示しています。

import numpy as np

class NLMS:
    def __init__(self, filter_order, step_size, epsilon=1e-8):
        """
        filter_order: フィルタの次数
        step_size: ステップサイズ（学習率）
        epsilon: ゼロ除算を防ぐための小さな定数
        """
        self.filter_order = filter_order
        self.step_size = step_size
        self.epsilon = epsilon
        self.weights = np.zeros(filter_order)

    def adapt(self, desired, input_signal):
        """
        desired: 望ましい出力信号（ターゲット信号）
        input_signal: フィルタの入力信号
        """
        input_signal = np.asarray(input_signal)
        if len(input_signal) < self.filter_order:
            input_signal = np.concatenate([np.zeros(self.filter_order - len(input_signal)), input_signal])
        
        x = input_signal[-self.filter_order:]  # 入力信号の一部（フィルタの次数に対応）
        y = np.dot(self.weights, x)  # 現在のフィルタ係数で出力を計算
        error = desired - y  # エラーを計算

        # フィルタ係数の更新
        norm_factor = np.dot(x, x) + self.epsilon
        self.weights += (self.step_size / norm_factor) * error * x

        return y, error

    def filter(self, input_signal):
        """
        input_signal: フィルタの入力信号
        """
        input_signal = np.asarray(input_signal)
        if len(input_signal) < self.filter_order:
            input_signal = np.concatenate([np.zeros(self.filter_order - len(input_signal)), input_signal])
        
        x = input_signal[-self.filter_order:]  # 入力信号の一部（フィルタの次数に対応）
        y = np.dot(self.weights, x)  # 現在のフィルタ係数で出力を計算
        
        return y

# 使用例
nlms = NLMS(filter_order=4, step_size=0.1)

input_signals = [-0.9, -0.7, 0.3, 0.5, -0.1, 0.2, -0.8, 0.7, -0.4, 0.6]
desired_signals = [-0.8, -0.6, 0.2, 0.4, -0.05, 0.15, -0.7, 0.6, -0.3, 0.5]

for desired, input_signal in zip(desired_signals, input_signals):
    y, error = nlms.adapt(desired, input_signal)
    print(f"入力: {input_signal}, 望ましい出力: {desired}, フィルタ出力: {y}, エラー: {error}")

# フィルタを使用して新しい入力信号を処理
new_input_signal = [0.2, -0.5, 0.3, -0.1]
output = nlms.filter(new_input_signal)
print(f"新しい入力信号のフィルタ出力: {output}")

import numpy as np

class FilteredXLMS:
    def __init__(self, filter_order, step_size, secondary_path, epsilon=1e-8):
        """
        filter_order: フィルタの次数
        step_size: ステップサイズ（学習率）
        secondary_path: 二次経路のインパルス応答
        epsilon: ゼロ除算を防ぐための小さな定数
        """
        self.filter_order = filter_order
        self.step_size = step_size
        self.epsilon = epsilon
        self.weights = np.zeros(filter_order)
        self.secondary_path = np.asarray(secondary_path)
        self.secondary_path_order = len(self.secondary_path)
        self.x_buffer = np.zeros(self.filter_order + self.secondary_path_order - 1)
    
    def adapt(self, desired, input_signal):
        """
        desired: 望ましい出力信号（ターゲット信号）
        input_signal: フィルタの入力信号
        """
        self.x_buffer = np.roll(self.x_buffer, -1)
        self.x_buffer[-1] = input_signal

        x = self.x_buffer[-self.filter_order:]  # 入力信号の一部（フィルタの次数に対応）
        filtered_x = np.convolve(self.x_buffer, self.secondary_path, mode='valid')[:self.filter_order]
        
        y = np.dot(self.weights, x)  # フィルタの出力は入力信号（x）に基づく
        error = desired - y  # エラーを計算

        norm_factor = np.dot(filtered_x, filtered_x) + self.epsilon
        self.weights += (self.step_size / norm_factor) * error * filtered_x

        return y, error

    def filter(self, input_signal):
        """
        input_signal: フィルタの入力信号
        """
        self.x_buffer = np.roll(self.x_buffer, -1)
        self.x_buffer[-1] = input_signal

        x = self.x_buffer[-self.filter_order:]  # 入力信号の一部（フィルタの次数に対応）
        y = np.dot(self.weights, x)  # フィルタの出力は入力信号（x）に基づく
        
        return y

# 使用例
filter_order = 4
step_size = 0.1
secondary_path = [0.2, 0.3, 0.5]

filtered_x_lms = FilteredXLMS(filter_order, step_size, secondary_path)

input_signals = [0.5, -0.4, 0.3, 0.7, -0.2, 0.4, -0.6, 0.8, -0.1, 0.5]
desired_signals = [0.45, -0.35, 0.25, 0.65, -0.15, 0.35, -0.55, 0.75, -0.05, 0.45]

for desired, input_signal in zip(desired_signals, input_signals):
    y, error = filtered_x_lms.adapt(desired, input_signal)
    print(f"入力: {input_signal}, 望ましい出力: {desired}, フィルタ出力: {y}, エラー: {error}")

# フィルタを使用して新しい入力信号を処理
new_input_signal = [0.2, -0.5, 0.3, -0.1]
for signal in new_input_signal:
    output = filtered_x_lms.filter(signal)
    print(f"新しい入力信号のフィルタ出力: {output}")

結果と注意点

Filtered-X LMSアルゴリズムにおいて、フィルタのタップ数（フィルタの長さ）を大きくしすぎると、フィルタの出力が-1.0～1.0の範囲を超える可能性があります。これは、フィルタのタップ数が増えると、フィルタの出力が入力信号の影響を強く受けやすくなるためです。 具体的には、以下の要因がフィルタの出力範囲に影響を与えます。 フィルタ係数の合計値: フィルタ係数の合計値が大きくなると、入力信号が増幅されることがあります。 入力信号の振幅: 入力信号が-1.0～1.0の範囲内である場合でも、フィルタ係数によってスケーリングされるため、出力信号がその範囲を超えることがあります。 タップ数の影響: フィルタのタップ数が増えると、フィルタ係数の合計値が増加する可能性があり、それによって出力信号が増幅されることがあります。 これを防ぐためには、フィルタ係数の正規化やクリッピングを行うことが有効です。フィルタ係数の正規化は、フィルタ係数の合計値が一定の範囲内に収まるように調整する方法です。また、出力信号が-1.0～1.0の範囲を超えた場合には、クリッピングによって信号をその範囲内に制限することができます。 以下に、フィルタ係数の正規化と出力信号のクリッピングを追加したFiltered-X LMSアルゴリズムの実装例を示します。

import numpy as np

class FilteredXLMS:
    def __init__(self, filter_order, step_size, secondary_path, epsilon=1e-8):
        """
        filter_order: フィルタの次数
        step_size: ステップサイズ（学習率）
        secondary_path: 二次経路のインパルス応答
        epsilon: ゼロ除算を防ぐための小さな定数
        """
        self.filter_order = filter_order
        self.step_size = step_size
        self.epsilon = epsilon
        self.weights = np.zeros(filter_order)
        self.secondary_path = np.asarray(secondary_path)
        self.secondary_path_order = len(self.secondary_path)
        self.x_buffer = np.zeros(self.filter_order + self.secondary_path_order - 1)
    
    def adapt(self, desired, input_signal):
        """
        desired: 望ましい出力信号（ターゲット信号）
        input_signal: フィルタの入力信号
        """
        self.x_buffer = np.roll(self.x_buffer, -1)
        self.x_buffer[-1] = input_signal

        x = self.x_buffer[-self.filter_order:]  # 入力信号の一部（フィルタの次数に対応）
        filtered_x = np.convolve(self.x_buffer, self.secondary_path, mode='valid')[:self.filter_order]
        
        y = np.dot(self.weights, x)  # フィルタの出力は入力信号（x）に基づく
        error = desired - y  # エラーを計算

        norm_factor = np.dot(filtered_x, filtered_x) + self.epsilon
        self.weights += (self.step_size / norm_factor) * error * filtered_x

        # フィルタ係数の正規化
        self.weights = self.weights / np.linalg.norm(self.weights)

        return y, error

    def filter(self, input_signal):
        """
        input_signal: フィルタの入力信号
        """
        self.x_buffer = np.roll(self.x_buffer, -1)
        self.x_buffer[-1] = input_signal

        x = self.x_buffer[-self.filter_order:]  # 入力信号の一部（フィルタの次数に対応）
        y = np.dot(self.weights, x)  # フィルタの出力は入力信号（x）に基づく
        
        # 出力信号のクリッピング
        y = np.clip(y, -1.0, 1.0)

        return y

# 使用例
filter_order = 4
step_size = 0.1
secondary_path = [0.2, 0.3, 0.5]

filtered_x_lms = FilteredXLMS(filter_order, step_size, secondary_path)

input_signals = [0.5, -0.4, 0.3, 0.7, -0.2, 0.4, -0.6, 0.8, -0.1, 0.5]
desired_signals = [0.45, -0.35, 0.25, 0.65, -0.15, 0.35, -0.55, 0.75, -0.05, 0.45]

for desired, input_signal in zip(desired_signals, input_signals):
    y, error = filtered_x_lms.adapt(desired, input_signal)
    print(f"入力: {input_signal}, 望ましい出力: {desired}, フィルタ出力: {y}, エラー: {error}")

# フィルタを使用して新しい入力信号を処理
new_input_signal = [0.2, -0.5, 0.3, -0.1]
for signal in new_input_signal:
    output = filtered_x_lms.filter(signal)
    print(f"新しい入力信号のフィルタ出力: {output}")