はじめに
前回は、DFTの実装とその結果解釈を行いました。今回は、そのDFT実装を用いて、種々の波形データを用意/入力し、より実践的な周波数解析を実施します。
波形データの準備
DFTに入力する波形データを用意します。
１個だけの周波数成分を含むcos波
| 変数名 | 記号 | 設定値 | 概要 |
|---|---|---|---|
| サンプリング時間 | To | 1.0sec | 入力データの取得開始から終了までの時間 |
| サンプリング数 | N | 100 | 入力データの取得数 |
| 周波数 | f | 10Hz | DFTで得たい数値 |
2個の周波数成分を含むcos波の合成波
| 変数名 | 記号 | 設定値 | 概要 |
|---|---|---|---|
| サンプリング時間 | To | 1.0sec | 入力データの取得開始から終了までの時間 |
| サンプリング数 | N | 100 | 入力データの取得数 |
| 周波数 | f | 10Hz, 19Hz | DFTで得たい数値 |
１個だけの周波数成分を含むcos波にノイズを含んだ波
| 変数名 | 記号 | 設定値 | 概要 |
|---|---|---|---|
| サンプリング時間 | To | 1.0sec | 入力データの取得開始から終了までの時間 |
| サンプリング数 | N | 100 | 入力データの取得数 |
| 周波数 | f | 10Hz | DFTで得たい数値 |
ノイズ成分として、各サンプリング値にランダムな値を加算しました。
DFTによる周波数解析
用意した波形をDFTした結果が以下になります。それぞれ周波数成分が抽出できています。
１個だけの周波数成分を含むcos波
設定した周波数(10Hz)の点にピークが出ています。
2個の周波数成分を含むcos波の合成波
設定した周波数(10, 19Hz)の点にピークが出ています。
１個だけの周波数成分を含むcos波にノイズを含んだ波
ノイズ成分があるため、いろいろな周波数成分が抽出されています。ノイズ成分は控えめにしていたので、抽出対象の波形の周波数成分(信号成分)が一番強く出ています。実際の入力波形データはノイズ成分が含まれることを考慮して、対策をしつつ周波数解析を行っていくはずで、ノイズ対策に関することも別途勉強できればと思います。今はあまり解っていません。窓関数とかそういうのだと思っています。
おわりに
DFTによって周波数成分を抽出することが確認できました。つまり、例えば波形データとして表現できる音を入力して、雑音除去などができるようになります。しかし、DFTは計算量がO(n^2)のため、波形データの数が増えると現実的な時間以内に計算することができなくなります。そこで計算量がO(nlog(n))である高速フーリエ変換(FFT)の出番です！ということで、いよいよ本丸のFFTを学んでいきたいと思います。
おまけ
１個だけの周波数成分を含む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)
