YU2TA7KA's BLOG ~take one step at a time~

派生開発、組み込み開発周りのこと。

離散フーリエ変換による波形データの周波数解析【No.3 離散フーリエ変換連載】

はじめに

前回は、DFTの実装とその結果解釈を行いました。今回は、そのDFT実装を用いて、種々の波形データを用意/入力し、より実践的な周波数解析を実施します。

波形データの準備

DFTに入力する波形データを用意します。

1個だけの周波数成分を含むcos波

変数名 記号 設定値 概要
サンプリング時間 To 1.0sec 入力データの取得開始から終了までの時間
サンプリング数 N 100 入力データの取得数
周波数 f 10Hz DFTで得たい数値


2個の周波数成分を含むcos波の合成波

変数名 記号 設定値 概要
サンプリング時間 To 1.0sec 入力データの取得開始から終了までの時間
サンプリング数 N 100 入力データの取得数
周波数 f 10Hz, 19Hz DFTで得たい数値

1個だけの周波数成分を含むcos波にノイズを含んだ波

変数名 記号 設定値 概要
サンプリング時間 To 1.0sec 入力データの取得開始から終了までの時間
サンプリング数 N 100 入力データの取得数
周波数 f 10Hz DFTで得たい数値

ノイズ成分として、各サンプリング値にランダムな値を加算しました。

DFTによる周波数解析

用意した波形をDFTした結果が以下になります。それぞれ周波数成分が抽出できています。

1個だけの周波数成分を含むcos波


設定した周波数(10Hz)の点にピークが出ています。

2個の周波数成分を含むcos波の合成波


設定した周波数(10, 19Hz)の点にピークが出ています。

1個だけの周波数成分を含むcos波にノイズを含んだ波


ノイズ成分があるため、いろいろな周波数成分が抽出されています。ノイズ成分は控えめにしていたので、抽出対象の波形の周波数成分(信号成分)が一番強く出ています。実際の入力波形データはノイズ成分が含まれることを考慮して、対策をしつつ周波数解析を行っていくはずで、ノイズ対策に関することも別途勉強できればと思います。今はあまり解っていません。窓関数とかそういうのだと思っています。

おわりに

DFTによって周波数成分を抽出することが確認できました。つまり、例えば波形データとして表現できる音を入力して、雑音除去などができるようになります。しかし、DFTは計算量がO(n^2)のため、波形データの数が増えると現実的な時間以内に計算することができなくなります。そこで計算量がO(nlog(n))である高速フーリエ変換(FFT)の出番です!ということで、いよいよ本丸のFFTを学んでいきたいと思います。



おまけ

1個だけの周波数成分を含むcos波にノイズを含んだ波を生成して、DFTするコードは以下です。

import numpy as np
import matplotlib.pyplot as plt

def generate_dft_matrix(size):
    n = size
    w = np.exp(-2j * np.pi * np.arange(n) / n * np.arange(n)[:, None])
    w.real[abs(w.real) <= 1.0e-10] = 0
    w.imag[abs(w.imag) <= 1.0e-10] = 0
    print(f"DFT matrix:\n{w}")
    return w

def dft(x):
    n = len(x)
    w = generate_dft_matrix(n)
    x_dft = np.dot(w, x) / n
    x_dft.real[abs(x_dft.real) <= 1.0e-10] = 0
    x_dft.imag[abs(x_dft.imag) <= 1.0e-10] = 0
    print(f"DFT result: {x_dft}")
    return x_dft

def plot_dft_result(dft_result, period, sampling_number):
    frequency_resolution = 1 / period
    k_array = np.arange(0, sampling_number, 1)
    frequency_array = np.array([k * frequency_resolution if k <= sampling_number / 2 else (k - sampling_number) * frequency_resolution for k in range(sampling_number)])
    plt.xticks(frequency_array[::10])
    plt.stem(frequency_array, dft_result.real, linefmt="blue")
    plt.stem(frequency_array, dft_result.imag, linefmt="green")
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Amplitude")
    plt.grid()
    plt.show()

def generate_cos_wave(frequency, amplitude, sampling_rate, duration):
    """Generate a cosine wave signal.

    Args:
        frequency (float): The frequency of the wave in Hz.
        amplitude (float): The amplitude of the wave.
        sampling_rate (int): The number of samples per second.
        duration (float): The duration of the signal in seconds.

    Returns:
        numpy.ndarray: A numpy array containing the signal.
    """
    # calculate the number of samples
    num_samples = int(sampling_rate * duration)

    # generate the time array
    time_array = np.arange(num_samples) / sampling_rate

    # generate the signal
    signal = amplitude * np.cos(2 * np.pi * frequency * time_array)

    return time_array, signal

def plot_signal(time, signal):
    plt.plot(time, signal)
    plt.xlabel('Time (seconds)')
    plt.ylabel('Amplitude')
    plt.title('Waveform')
    plt.show()

# Unit data
f = 10
amplitude = 1
period = 1 # period[sec]
sampling_rate = 100

# Generate the signal
time, signal = generate_cos_wave(f, amplitude, sampling_rate, period)

# Add noise component
noise = np.random.randn(len(time)) * 0.2

plot_signal(time, signal + noise)

# DFT of the signal
dft_result = dft(signal + noise)

# Plot DFT result
plot_dft_result(dft_result, period, sampling_rate)