From 09779a182f033ba767a916caa8be565b268c3265 Mon Sep 17 00:00:00 2001 From: rkatotmu Date: Mon, 10 Jul 2023 16:12:59 +0900 Subject: [PATCH] =?UTF-8?q?4=E7=AB=A0=E5=89=8D=E5=8D=8A=E3=82=92=E8=BF=BD?= =?UTF-8?q?=E5=8A=A0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- rkato/chapter04/q01.py | 17 +++++++ rkato/chapter04/q02.py | 33 ++++++++++++ rkato/chapter04/q03.py | 33 ++++++++++++ rkato/chapter04/q04.py | 111 +++++++++++++++++++++++++++++++++++++++++ rkato/chapter04/q05.py | 60 ++++++++++++++++++++++ 5 files changed, 254 insertions(+) create mode 100644 rkato/chapter04/q01.py create mode 100644 rkato/chapter04/q02.py create mode 100644 rkato/chapter04/q03.py create mode 100644 rkato/chapter04/q04.py create mode 100644 rkato/chapter04/q05.py diff --git a/rkato/chapter04/q01.py b/rkato/chapter04/q01.py new file mode 100644 index 0000000..0d3f7d4 --- /dev/null +++ b/rkato/chapter04/q01.py @@ -0,0 +1,17 @@ +import numpy as np + + +# L-Sこのゼロ詰め　Sの倍数になるようにゼロ詰め +def padding(x, L, S): + x_pad = np.pad(x, [L - S, L - S]) + re = np.mod(len(x_pad), S) + if re != 0: + x_pad = np.pad(x_pad, [0, S - re]) + return x_pad + + +##### 確認コード ##### +L = 4 +S = 3 +x = np.ones(8) +print(padding(x, L, S)) diff --git a/rkato/chapter04/q02.py b/rkato/chapter04/q02.py new file mode 100644 index 0000000..75acfa7 --- /dev/null +++ b/rkato/chapter04/q02.py @@ -0,0 +1,33 @@ +import numpy as np + + +# L-Sこのゼロ詰め　Sの倍数になるようにゼロ詰めする関数 +def padding(x, L, S): + x_pad = np.pad(x, [L - S, L - S]) # 先頭と末尾をゼロ詰め + re = np.mod(len(x), S) # Sの倍数か調べる + if re != 0: + x_pad = np.pad(x_pad, [0, S - re]) # Sの余り分付け足す + return x_pad + + +def frame_div(x, L, S): + x_pad = padding(x, L, S) + T = int(np.floor((x_pad.size - L) / S)) + 1 + x_t = np.array([x_pad[t * S : t * S + L] for t in range(T)]) + return x_t + + +# # フレーム分割 +# def frame_div(x, L, S): +# x_pad = padding(x, L, S) +# T = int(np.floor((len(x_pad) - L) / S)) + 1 # T行分抽出 +# x_t = np.array([x_pad[t * S : t * S + L]] for t in range(T)) +# return x_t + + +##### 確認コード ##### +L = 4 +S = 3 +x = np.ones(8) +x_t = frame_div(x, L, S) +print(x_t) diff --git a/rkato/chapter04/q03.py b/rkato/chapter04/q03.py new file mode 100644 index 0000000..4045181 --- /dev/null +++ b/rkato/chapter04/q03.py @@ -0,0 +1,33 @@ +import numpy as np + + +# L-Sこのゼロ詰め　Sの倍数になるようにゼロ詰めする関数 +def padding(x, L, S): + x_pad = np.pad(x, [L - S, L - S]) # 先頭と末尾をゼロ詰め + re = np.mod(len(x), S) # Sの倍数か調べる + if re != 0: + x_pad = np.pad(x_pad, [0, S - re]) # Sの余り分付け足す + return x_pad + + +# フレーム分割 +def frame_div(x, L, S): + x_pad = padding(x, L, S) + T = int(np.floor((x_pad.size - L) / S)) + 1 + x_t = np.array([x_pad[t * S : t * S + L] for t in range(T)]) + return x_t + + +def stft(x, L, S, w): + x_t = frame_div(x, L, S) + T = len(x_t) + X = np.array([np.fft.rfft(x_t[t] * w) for t in range(T)], dtype="complex") + return X.T + + +#####確認コード##### +L = 4 +S = 3 +x = np.ones(8) +wnd = np.hamming(L) +x_stft = stft(x, L, S, wnd) diff --git a/rkato/chapter04/q04.py b/rkato/chapter04/q04.py new file mode 100644 index 0000000..02bae2f --- /dev/null +++ b/rkato/chapter04/q04.py @@ -0,0 +1,111 @@ +from pathlib import Path +from typing import Union +import numpy as np +from matplotlib import pyplot as plt + +config = { + "figure.subplot.bottom": 0.15, + "figure.subplot.left": 0.15, + "figure.figsize": [10, 8], + "font.size": 24, + "mathtext.fontset": "cm", + "pdf.fonttype": 42, + "legend.borderaxespad": 1, + "lines.linewidth": 2, + "savefig.transparent": True, +} + + +LABELS = { + 0: ("viridis", "Power spectrogram"), + 1: ("viridis", "Phase spectrogram"), + 2: (plt.cm.Reds, "Cosine of phase spectrogram"), + 3: (plt.cm.Blues, "Sine of hase spectrogram"), +} + + +# L-Sこのゼロ詰め　Sの倍数になるようにゼロ詰めする関数 +def pad(x, L, S): + x_pad = np.pad(x, [L - S, L - S]) + re = np.mod(x_pad.size, S) + if re != 0: + x_pad = np.pad(x_pad, [0, S - re]) + return x_pad + + +# フレーム分割 +def frame_div(x, L, S): + x_pad = pad(x, L, S) + T = int(np.floor((x_pad.size - L) / S)) + 1 + x_t = np.array([x_pad[t * S : t * S + L] for t in range(T)]) + return x_t + + +# 短時間フーリエ変換 +def stft(x, L, S, w): + x_div = frame_div(x, L, S) + T = x_div.shape[0] + X = np.empty([L // 2 + 1, T], dtype="complex") + for t in range(0, T): + X[:, t] = np.fft.rfft(x_div[t, :] * w) + return X + + +def main(output_dir: Union[Path, str]): + out_p = Path(output_dir) + + # 正弦波の生成　振幅:１　音高:440Hz 信号長:0.1秒　サンプリングレート:16kHz + A = 1 # amplitude + f = 440 # frequency (Hz) + sec = 0.1 # time (sec) + fs = 16000 # sampling rate (Hz) + t = np.arange(0, sec, 1 / fs) + y = A * np.sin(2 * np.pi * f * t) + # 窓幅L,シフト幅S + L = 1000 + S = L // 2 + + w = np.hamming(L) + y_stft = stft(y, L, S, w) + print(f"y_stft.shape: {y_stft.shape}") + spec = [ + np.abs(y_stft), + np.angle(y_stft), + np.cos(np.angle(y_stft)), + np.sin(np.angle(y_stft)), + ] + + # # ここからグラフ描画 + plt.rcParams.update(config) + fig, ax = plt.subplots(4, 1) + for i in range(len(spec)): + # ax[i].pcolormesh(spec[i]) + ax[i].imshow( + spec[i], + cmap=LABELS[i][0], + interpolation="nearest", + aspect="auto", + origin="lower", + extent=[0, y_stft.shape[1], 0, y_stft.shape[0]], + ) + ax[i].set_title(LABELS[i][1]) + ax[i].set_rasterized(True) + fig.supxlabel("Frame") + fig.supylabel("Frequency (Hz)") + plt.tight_layout() + fig.subplots_adjust( + left=0.13, right=1 - 0.13, bottom=0.11, top=1 - 0.11, hspace=0.8 + ) + plt.savefig(out_p / "04.pdf", dpi=300) + plt.clf() + plt.close() + + +if __name__ == "__main__": + out_p = Path.cwd() / "outputs" + if not out_p.exists(): + out_p.mkdir(parents=True) + + main(out_p) + + print("Finished") diff --git a/rkato/chapter04/q05.py b/rkato/chapter04/q05.py new file mode 100644 index 0000000..6781172 --- /dev/null +++ b/rkato/chapter04/q05.py @@ -0,0 +1,60 @@ +from pathlib import Path +from typing import Union +import numpy as np +import matplotlib.pyplot as plt + + +config = { + "figure.subplot.bottom": 0.15, + "figure.subplot.left": 0.15, + "figure.figsize": [10, 8], + "font.size": 24, + "mathtext.fontset": "cm", + "pdf.fonttype": 42, + "legend.borderaxespad": 1, + "lines.linewidth": 2, + "savefig.transparent": True, +} +# 合成窓 +def comp_wnd(wnd, S): + L = wnd.size + Q = L // S + for l in range(0, L): + sigma = 0 + for m in range(0, Q): + sigma += wnd[l - m * S] ** 2 + wnd_s = wnd / sigma + + return wnd_s + + +def main(output_dir: Union[Path, str]): + out_p = Path(output_dir) + + L = 1000 + S = 500 + wnd = np.hamming(L) + cwnd = comp_wnd(wnd, S) + + plt.rcParams.update(config) + fig, ax = plt.subplots(2, 1) + ax[0].plot(wnd) + ax[0].grid() + ax[0].set_title("Hamming window") + ax[1].plot(cwnd) + ax[1].set_title("Composite window") + ax[1].grid() + plt.tight_layout() + plt.savefig(out_p / "05.pdf") + plt.clf() + plt.close() + + +if __name__ == "__main__": + out_p = Path.cwd() / "outputs" + if not out_p.exists(): + out_p.mkdir(parents=True) + + main(out_p) + + print("Finished") \ No newline at end of file