-
-
Save takatakamanbou/476727bd3943060a723d0cc4ba81cfba to your computer and use it in GitHub Desktop.
PIP2020-09-note2.ipynb
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "nbformat": 4, | |
| "nbformat_minor": 0, | |
| "metadata": { | |
| "colab": { | |
| "name": "PIP2020-09-note2.ipynb", | |
| "provenance": [], | |
| "collapsed_sections": [], | |
| "include_colab_link": true | |
| }, | |
| "kernelspec": { | |
| "name": "python3", | |
| "display_name": "Python 3" | |
| } | |
| }, | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "view-in-github", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "<a href=\"https://colab.research.google.com/gist/takatakamanbou/476727bd3943060a723d0cc4ba81cfba/pip2020-09-note2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "NI3mY1FGZWgV", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "# 2020年度パターン情報処理第9回講義資料その2\n", | |
| "\n", | |
| "この科目のウェブサイトへ https://www-tlab.math.ryukoku.ac.jp/wiki/?PIP/2020\n", | |
| "\n", | |
| "" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "ffaGVK2o5VWe", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "## はじめに" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "0IpujRN_Wuib", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### これは何？\n", | |
| "\n", | |
| "これは，Google Colaboratory（以下 Google Colab) という，Google が提供しているサービスを利用して作成した資料です．Notebook と呼びます．クラウド上に仮想的に構築された Linux マシン上で Python のプログラムを実行することができます．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "QRNjWashW4pr", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### Notebook の動かし方\n", | |
| "\n", | |
| "この Notebook では，上の方のセルで作った変数や関数を後のセルで使うことがあります．そのため，上の方のセルを実行せずに下の方を実行すると，エラーになることがあります．上から順に実行していきましょう．\n", | |
| "\n", | |
| "また，メニューの「ランタイム」から，「すべてのセルを実行」したりすることも可能です．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "1jfSXIY6VNxP", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "## 準備\n", | |
| "\n", | |
| "以下のセルは，プログラムを実行するための準備を行うためのものです．このセルを実行してから，先へ進みましょう．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "IFYQk7ONZENz", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# 科学技術計算のライブラリ NumPy のモジュールを np という名前で使えるようにする\n", | |
| "import numpy as np\n", | |
| "# 科学技術計算のライブラリ SciPy の中の WAVE ファイルを扱うモジュール を wavfile という名前で使えるようにする\n", | |
| "import scipy.io.wavfile as wavfile\n", | |
| "# グラフを描くためのライブラリ matplotlib の pyplot を plt という名前でインポート\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "\n", | |
| "# matplotlib が描くグラフの文字サイズをデフォルト 12pt から 18pt へ\n", | |
| "plt.rcParams[\"font.size\"] = 18" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "807R4hQI5oH9", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# WAVE ファイルを読み込む関数\n", | |
| "#\n", | |
| "def readWAVE(filename):\n", | |
| " \n", | |
| " framerate, data = wavfile.read(filename)\n", | |
| "\n", | |
| " # チャンネル数とフレーム数（データ点の数）を求める\n", | |
| " if data.ndim == 1:\n", | |
| " nchannels = 1\n", | |
| " else:\n", | |
| " nchannels = data.shape[1]\n", | |
| " nframes = data.shape[0]\n", | |
| " \n", | |
| " print('filename = ', filename)\n", | |
| " print('nchannels = ', nchannels) # チャンネル数\n", | |
| " print('framerate = ', framerate) # 標本化周波数\n", | |
| " print('nframes = ', nframes) # フレーム数\n", | |
| " print('duration = {0:.2f}[sec]'.format(nframes / framerate)) # 長さ（秒単位）\n", | |
| " print('dtype = ', data.dtype) # データ型（量子化ビット数に対応）\n", | |
| "\n", | |
| " assert data.dtype == 'uint8' or data.dtype == 'int16' or data.dtype == 'int32' or data.dtype == 'float32'\n", | |
| " \n", | |
| " # 最初の10秒分だけ取り出す（元がそれより短ければそのまま）\n", | |
| " nframesNew = min(framerate * 10, nframes) \n", | |
| " if nchannels == 1:\n", | |
| " dataNew = data[:nframesNew]\n", | |
| " else:\n", | |
| " # 多チャンネル信号なら0番目と1番目の平均値を取り出す\n", | |
| " if data.dtype == 'float32': # 浮動小数点数のときは [0, 1] の値なので普通に足して2で割る\n", | |
| " dataNew = (data[:nframesNew, 0] + data[:nframesNew, 1])/2\n", | |
| " else: # 整数型のときはオーバーフローする可能性があるので，いったん64bit整数にしてから\n", | |
| " data64 = (data[:nframesNew, 0].astype(np.int64) + data[:nframesNew, 1].astype(np.int64))//2\n", | |
| " dataNew = data64.astype(data.dtype)\n", | |
| " \n", | |
| " return framerate, dataNew\n", | |
| "\n", | |
| "\n", | |
| "# WAVE ファイルを書き込む関数\n", | |
| "#\n", | |
| "def writeWAVE(filename, framerate, data):\n", | |
| "\n", | |
| " wavfile.write(filename, framerate, data)" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "82S39H2VB8oX", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "## ★ 7.2 標本化とエイリアシング" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "39pG5akMPrP8", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### 講義資料 Q1\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "w9YIPxtAqXsN", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "### 下ごしらえ\n", | |
| "#\n", | |
| "tr0 = np.linspace(0.0, 6.0*np.pi, num=1000)\n", | |
| "tr1 = np.linspace(0.0, 6.0*np.pi, num=5)\n", | |
| "tr2 = np.linspace(0.0, 6.0*np.pi, num=7)\n", | |
| "tr3 = np.linspace(0.0, 6.0*np.pi, num=8)\n", | |
| "s1 = np.sin(tr1)\n", | |
| "s2 = np.sin(tr2)\n", | |
| "s3 = np.sin(tr3)" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "vqMvTsOvFWYZ", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "fig, ax = plt.subplots(3, facecolor=\"white\", figsize=(8, 12))\n", | |
| "ax[0].plot(tr1, s1, \"o\", color=\"red\")\n", | |
| "ax[0].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].plot(tr2, s2, \"o\", color=\"red\")\n", | |
| "ax[1].set_ylim(-1.1, 1.1)\n", | |
| "ax[2].plot(tr3, s3, \"o\", color=\"red\")\n", | |
| "ax[2].set_ylim(-1.1, 1.1)\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "e0KK4HPb6X4T", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "上図は，ある正弦波をそれぞれ異なる標本化間隔（左から順に$\\frac{3}{2}\\pi, \\pi, \\frac{5}{6}\\pi$）で標本化したものである．これらの点をつなぐ曲線を描いて，標本化される前の正弦波がどのようなものだったか予想してみなさい．\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "OH_aJw4RD4_W", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "fig, ax = plt.subplots(2, facecolor=\"white\", figsize=(8, 8))\n", | |
| "ax[0].plot(tr1, s1, \"o\", color=\"red\")\n", | |
| "ax[0].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].plot(tr1, s1, \"o\", color=\"red\")\n", | |
| "ax[1].plot(tr0, 0*tr0, \"-\", color=\"gray\", label=\"0\")\n", | |
| "ax[1].plot(tr0, -np.sin(tr0/3), \"-\", color=\"blue\", label=\"$-\\\\sin\\\\frac{1}{3}t$\")\n", | |
| "ax[1].plot(tr0, np.sin(tr0), \"-\", color=\"green\", label=\"$\\\\sin t$\")\n", | |
| "ax[1].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].legend()\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "colab_type": "code", | |
| "id": "rzPBczjCIqi6", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "fig, ax = plt.subplots(2, facecolor=\"white\", figsize=(8, 8))\n", | |
| "ax[0].plot(tr2, s2, \"o\", color=\"red\")\n", | |
| "ax[0].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].plot(tr2, s2, \"o\", color=\"red\")\n", | |
| "ax[1].plot(tr0, 0*tr0, \"-\", color=\"gray\", label=\"0\")\n", | |
| "ax[1].plot(tr0, -np.sin(tr0/3), \"-\", color=\"blue\", label=\"$-\\\\sin\\\\frac{1}{3}t$\")\n", | |
| "ax[1].plot(tr0, np.sin(tr0), \"-\", color=\"green\", label=\"$\\\\sin t$\")\n", | |
| "ax[1].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].legend()\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "O93uFQ30GQ8W", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "fig, ax = plt.subplots(2, facecolor=\"white\", figsize=(8, 8))\n", | |
| "ax[0].plot(tr3, s3, \"o\", color=\"red\")\n", | |
| "ax[0].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].plot(tr3, s3, \"o\", color=\"red\")\n", | |
| "ax[1].plot(tr0, 0*tr0, \"-\", color=\"gray\", label=\"0\")\n", | |
| "ax[1].plot(tr0, -np.sin(tr0/3), \"-\", color=\"blue\", label=\"$-\\\\sin\\\\frac{1}{3}t$\")\n", | |
| "ax[1].plot(tr0, np.sin(tr0), \"-\", color=\"green\", label=\"$\\\\sin t$\")\n", | |
| "ax[1].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].legend()\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "6s4TK05p6dMy", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "実は，元の信号は周期 $2\\pi$ の正弦波 $\\sin{t}$ だった．標本化間隔を $\\frac{3}{2}\\pi$ や $\\pi$ にすると，そこから元の信号の形を推測するのはむずかしい\b．標本化間隔 $\\frac{5}{6}\\pi$ だといけそう．\n", | |
| "\n", | |
| "⇒ 正弦波の場合，元の信号の周期の半分より短い周期で標本化すれば，そこから元の信号の形がわかる．↑の例では，元の信号の周期が $2\\pi$ で，$\\frac{5}{6}\\pi < \\pi$．\n", | |
| "\n", | |
| "以下は，標本化周期を $0.95\\pi$ とした場合．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "O3uxOgdlNbeS", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "tr4 = np.arange(0.0, 6.0*np.pi+0.01, 0.95*np.pi)\n", | |
| "s4 = np.sin(tr4)\n", | |
| "\n", | |
| "fig, ax = plt.subplots(2, facecolor=\"white\", figsize=(8, 8))\n", | |
| "ax[0].plot(tr4, s4, \"o\", color=\"red\")\n", | |
| "ax[0].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].plot(tr4, s4, \"o\", color=\"red\")\n", | |
| "ax[1].plot(tr0, np.sin(tr0), \"-\", color=\"green\", label=\"$\\\\sin t$\")\n", | |
| "ax[1].set_ylim(-1.1, 1.1)\n", | |
| "ax[1].legend()\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "35atLHsK6iOI", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### エイリアシング" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "Q-PeW0cfOnfS", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "### 設定\n", | |
| "#\n", | |
| "sampfreq = 8000 # 標本化周波数 8000 Hz\n", | |
| "valmax = 2**15 - 1 # 16ビット量子化したデータを記録する．そのときのとり得る値の最大値\n", | |
| "t = np.arange(0, 2, 1/sampfreq)\n", | |
| "\n", | |
| "fList = [2000, 3000, 3800, 4000, 4200, 5000] # 作成する正弦波の周波数" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "1R4gsGcOFahy", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "### 正弦波の WAVE データを作成\n", | |
| "#\n", | |
| "signal = np.empty((len(fList), len(t)), dtype=np.int16)\n", | |
| "for i, f in enumerate(fList):\n", | |
| " raw = 0.5*valmax*np.sin(2*np.pi*f*t) # f [Hz] の正弦波を生成\n", | |
| " signal[i, :] = raw\n", | |
| " writeWAVE(f'sin{f:04d}.wav', sampfreq, signal[i, :])" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "uPwko5FB95TG", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "以下のセルのコメントの指示にしたがって修正してから実行すると，上記で作った音のファイル（WAVE形式，拡張子 `.wav`）をダウンロードできます．自分のPCで鳴らしてみよう．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "7XdEkvRH935z", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# ファイル一覧\n", | |
| "! ls -l\n", | |
| "\n", | |
| "# ファイルをダウンロード\n", | |
| "from google.colab import files\n", | |
| "if 1 == 0: # ←の 0 を 1 に修正\n", | |
| " for f in fList:\n", | |
| " files.download(f'sin{f:04d}.wav')" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "28wmQMPmSkax", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# グラフに描いてみる\n", | |
| "\n", | |
| "nt = sampfreq//1000\n", | |
| "t2 = np.linspace(0, nt/sampfreq, num=1000)\n", | |
| "\n", | |
| "fig, ax = plt.subplots(len(fList), facecolor=\"white\", figsize=(8, 3*len(fList)))\n", | |
| "for i, f in enumerate(fList):\n", | |
| " ax[i].plot(t[:nt], signal[i, :nt], \"o\", color=\"red\")\n", | |
| " ax[i].plot(t2, 0.5*valmax*np.sin(2*np.pi*f*t2), \"-\", color=\"blue\", label=f\"{f}[Hz]\")\n", | |
| " if f == sampfreq//2:\n", | |
| " ax[i].plot(t2, 0*t2, \"-\", color=\"gray\", label=\"0\")\n", | |
| " elif f > sampfreq//2:\n", | |
| " fd = sampfreq - f\n", | |
| " ax[i].plot(t2, -0.5*valmax*np.sin(2*np.pi*fd*t2), \"-\", color=\"gray\", label=f\"{fd}[Hz]\") \n", | |
| " ax[i].legend()\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "d9MsZDVBqF_Z", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### もふもふの音の再標本化" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "VXmXseUJdvUn", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "以前使ったもふもふの音を再び使ってみます．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "colab_type": "code", | |
| "id": "aHxP1ppLqR9N", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# WAVE ファイルをダウンロード\n", | |
| "#\n", | |
| "### こちらのサイトの素材を利用させてもらってます http://www.ne.jp/asahi/music/myuu/wave/wave.htm\n", | |
| "\n", | |
| "! wget -nc http://www.ne.jp/asahi/music/myuu/wave/cat1.wav\n", | |
| "! ls -l\n", | |
| "fnCat = 'cat1.wav'\n", | |
| "\n", | |
| "import os\n", | |
| "\n", | |
| "if not os.path.exists(fnCat):\n", | |
| " print(f'{fnCat}のダウンロードがうまくできていないようです')" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "colab_type": "code", | |
| "id": "yCzndUjYqR9S", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# WAVEファイルの内容を読み込む\n", | |
| "sampfreqCat, signalCat = readWAVE(fnCat)\n", | |
| "N = signalCat.shape[0]" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "colab_type": "code", | |
| "id": "hTPNzVets3k7", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "### 高速フーリエ変換して振幅スペクトルを求め，グラフに描く\n", | |
| "#\n", | |
| "freq = np.arange(N) * sampfreqCat / N # 横軸の値．単位は Hz\n", | |
| "fig, ax = plt.subplots(2, facecolor=\"white\", figsize=(8, 6))\n", | |
| "C = np.fft.fft(signalCat)/N\n", | |
| "spec_amp = np.abs(C)\n", | |
| "ax[0].plot(freq, spec_amp, \"-\", color=\"red\", label=\"cat1\")\n", | |
| "ax[0].set_xlim(0, sampfreqCat//2)\n", | |
| "ax[0].legend()\n", | |
| "ax[1].plot(freq[:N//6], spec_amp[:N//6], \"-\", color=\"red\", label=\"cat1\")\n", | |
| "ax[1].set_xlim(0, freq[N//6])\n", | |
| "ax[1].legend()\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "-zlBV4o0d2i3", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "この音は 44100 [Hz] で標本化されており，振幅スペクトルは上記のようになっています（下段は，描く周波数範囲を狭めたもの）．4500 [Hz] くらいまでの成分が含まれているようです．\n", | |
| "\n", | |
| "これを無理やり 6300 [Hz] で標本化し直してみます．6300/2 = 3150 [Hz] 以上の高い周波数の成分がちゃんと表せなくなりそうですが...．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "FZ5trDMHYPMs", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "### 7個ごとにサンプリング => 標本化周波数 44100/7 = 6300 [Hz] で再標本化\n", | |
| "#\n", | |
| "signalCat2 = signalCat[::7]\n", | |
| "sampfreqCat2 = sampfreqCat // 7\n", | |
| "N2 = len(signalCat2)" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "CnZSCm7YYtXK", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "### 高速フーリエ変換して振幅スペクトルを求め，グラフに描く\n", | |
| "#\n", | |
| "freq2 = np.arange(N2) * sampfreqCat2 / N2 # 横軸の値．単位は Hz\n", | |
| "fig, ax = plt.subplots(2, facecolor=\"white\", figsize=(8, 6))\n", | |
| "C = np.fft.fft(signalCat2)/N2\n", | |
| "spec_amp2 = np.abs(C)\n", | |
| "ax[0].plot(freq, spec_amp, \"-\", color=\"red\", label=\"cat1\")\n", | |
| "ax[0].set_xlim(0, sampfreqCat2//2)\n", | |
| "ax[0].legend()\n", | |
| "ax[1].plot(freq2, spec_amp2, \"-\", color=\"red\", label=\"cat1(resampled)\")\n", | |
| "ax[1].set_xlim(0, sampfreqCat2//2)\n", | |
| "ax[1].legend()\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "0_6acvvMcTTF", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "上の図（cat1）は元の信号の振幅スペクトル，下の図（cat1_resampled）はそれを低い標本化周波数で再標本化したものです．エイリアシングが起こっており，ところどころに元の信号のスペクトルにはなかった成分が出現しているのがわかります．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "colab_type": "text", | |
| "id": "9_he9-nkqR9V" | |
| }, | |
| "source": [ | |
| "以下のセルのコメントの指示にしたがって修正してから実行すると，2つの音のファイル（WAVE形式，拡張子 `.wav`）をダウンロードできます．自分のPCで鳴らしてみよう．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "G7CcC1_dZj6Z", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "fnCat2 = 'cat1_resampled.wav'\n", | |
| "writeWAVE(fnCat2, sampfreqCat2, signalCat2)\n", | |
| "! ls -l\n", | |
| "\n", | |
| "# ファイルをダウンロード\n", | |
| "from google.colab import files\n", | |
| "if 1 == 0: # ←の 0 を 1 に修正\n", | |
| " files.download(fnCat) # 元のもふもふの声\n", | |
| " files.download(fnCat2) # 再標本化したもの" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "WPcLHE0ZKgsB", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "## ★ 7.3 離散フーリエ変換" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "L8v386B-K1q6", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### 講義資料に載ってる $N=4$ の例" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "ZRw0-PD_MIeL", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "講義資料 Q3" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "bDZFxmL6wmP4", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# 講義資料 3/4 Q3\n", | |
| "f = np.array([4.0, 4.0, 0.0, 0.0])\n", | |
| "N = f.shape[0]\n", | |
| "# 離散フーリエ変換（資料の式の N 倍の値が出てくるので，1/N してる）\n", | |
| "C = np.fft.fft(f) / N\n", | |
| "# フーリエ係数 C[k] を表示． ここでは虚数単位が i ではなく j と表されていることに注意．\n", | |
| "print(C)" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "JPLkKpaKK58R", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# 逆変換\n", | |
| "ft = np.fft.ifft(C) * N\n", | |
| "print(ft)" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "LrlLTdIxMMbe", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "講義資料 Q5" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "eIu34JCxLDOI", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# 講義資料 3/4 Q5 の例\n", | |
| "f = np.array([0.0, 1.0, 0.0, 2.0])\n", | |
| "# 離散フーリエ変換\n", | |
| "C = np.fft.fft(f) / f.shape[0]\n", | |
| "# フーリエ係数 C[k] を表示． ここでは虚数単位が i ではなく j と表されていることに注意．\n", | |
| "print(\"フーリエ係数:\", C)\n", | |
| "# 逆変換\n", | |
| "ft = np.fft.ifft(C) * f.shape[0]\n", | |
| "print(\"逆変換の結果:\", ft)" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "fT9XkLRwMP_H", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### 気温データのフーリエ解析 (1)\n", | |
| "\n", | |
| "第8回講義資料の気温データ．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "ZWDqZMTbLV7k", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# 気温のデータを入手する\n", | |
| "! wget -nc https://www-tlab.math.ryukoku.ac.jp/~takataka/course/PIP/temp2004.txt\n", | |
| "! ls -l\n", | |
| "# 入手したファイルを読み込んで temp2004 という変数に入れておく\n", | |
| "temp2004 = np.loadtxt('temp2004.txt')\n", | |
| "N = temp2004.shape[0]\n", | |
| "print(temp2004)\n", | |
| "print(temp2004.shape) # 366次元ベクトル" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "b4ZElzy3_4-u", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "#### 気温変化のグラフを描いてみる" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "8LWqJNR1Ax1A", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# t と temp2004 でグラフを描く\n", | |
| "t = np.arange(0, 366)\n", | |
| "fig, ax = plt.subplots(facecolor=\"white\", figsize=(8, 6))\n", | |
| "ax.plot(t, temp2004, \"-\", label=\"temp2004\")\n", | |
| "ax.legend()\n", | |
| "ax.set_xlim([0, 365]) # 横軸の範囲\n", | |
| "ax.set_ylim([-5, 35]) # 縦軸の範囲\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "j9sEb1OpyjwC", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "#### このデータを離散フーリエ変換してみる" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "Gz4Kg7qkvr9n", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# temp2004 のフーリエ変換\n", | |
| "C = np.fft.fft(temp2004) / N\n", | |
| "print(C.shape) # C は366次元の複素ベクトル" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "n-_1d_JEA6dI", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# Cの0番目は定数成分の係数\n", | |
| "print(C[0])" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "XCLLOcSDxac9", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# Cの1番目は周期366日の正弦波成分の係数\n", | |
| "print(C[1])" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "wLGMHxd_0_tH", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# C[1] の絶対値と偏角から周期366日の正弦波成分の振幅と位相を求める\n", | |
| "amp = np.abs(C[1])*2\n", | |
| "theta = np.angle(C[1])\n", | |
| "print('振幅 = ', amp, '位相 = ', theta)" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "6DKZ4TAmNZRR", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "$k=1$の波は，周期 $N=266$，振幅 $2|C_1|$，位相 $\\arg C_1$ の正弦波，つまり\n", | |
| "$$\n", | |
| "2|C_1|\\cos \\left(\\frac{2\\pi}{366}t + \\arg C_1 \\right)\\\\\n", | |
| "$$\n", | |
| "上記の振幅と位相の値を代入して，さらに $C_0$ を加えた式でグラフを描くと...． " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "YwZvoBBJc5aZ", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# 上記の振幅，位相をもつ周期366の正弦波を生成\n", | |
| "sig = amp * np.cos(2*np.pi/366*t + theta)\n", | |
| "# それに定数成分を加える\n", | |
| "sig += C[0].real" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "FBzcZiKi1sNo", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# グラフを描く\n", | |
| "fig, ax = plt.subplots(facecolor=\"white\", figsize=(8, 6))\n", | |
| "ax.plot(t, temp2004, \"-\", label = \"temp2004\")\n", | |
| "ax.plot(t, sig, \"-\", label = \"sig\")\n", | |
| "ax.legend()\n", | |
| "ax.set_xlim([0, 365]) # 横軸の範囲\n", | |
| "ax.set_ylim([-5, 35]) # 縦軸の範囲\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "7ePQKt2IP3mu", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "1年（366日）周期の気温変化の成分のみを眺めることができる．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "69IeLxp5dGMv", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### 気温データのフーリエ解析 (2)\n", | |
| "\n", | |
| "#### 一部の係数を 0 に置き換えて逆変換してみる\n", | |
| "\n", | |
| "ここでやってることの詳しいところは，次回説明の予定です．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "bXEkMdkNdFbh", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "### 0, 1, 2, ..., N-1 個の C の要素のうち，0から H番目までと(N-H+1)からN番目までを 0 にしたものを返す関数\n", | |
| "#\n", | |
| "def lowpass(C, H):\n", | |
| " N = C.shape[0] # C の要素数を N とする\n", | |
| " if H >= N//2:\n", | |
| " return C\n", | |
| " else:\n", | |
| " Ct = np.copy(C)\n", | |
| " Ct[H+1:N-H] = 0\n", | |
| " return Ct" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "qqmtLBnQegdK", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# 5番目までの成分のみを残す\n", | |
| "H = 5\n", | |
| "Ct = lowpass(C, H)\n", | |
| "# 逆変換\n", | |
| "rec = np.fft.ifft(Ct).real * Ct.shape[0]\n", | |
| "#print(rec)" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "ws1gMO3cfO6s", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "fig, ax = plt.subplots(facecolor=\"white\", figsize=(8, 6))\n", | |
| "ax.plot(t, temp2004, \"-\", label = \"temp2004\")\n", | |
| "ax.plot(t, rec, \"-\", label = f\"rec(H={H})\")\n", | |
| "ax.legend()\n", | |
| "ax.set_xlim([0, 365]) # 横軸の範囲\n", | |
| "ax.set_ylim([-5, 35]) # 縦軸の範囲\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "Tm9TwrfqSc1Q", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "上記は，$k = 0, 1, \\dots , 5$ までの成分のみを使って元の信号を表そうとした（他の成分の係数は $0$ とみなした）もの．$k$ の小さい成分は周期の大きい波なので，ゆっくりした気温変化は表せているけれど，短い周期での変化は再現できていない．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "_Plktu3sgc-b", | |
| "colab_type": "code", | |
| "cellView": "both", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "#@title スライダを動かすと H （何番目の成分までつかうか）を変えて近似させられます．H を変えると近似結果はどう変わるか観察しよう { run: \"auto\" }\n", | |
| "H = 0 #@param {type:\"slider\", min:0, max:20, step:1}\n", | |
| "# H をもっと大きくしてみたいひとは，↑のmax: の数字を 366 にしてみたらよい\n", | |
| "\n", | |
| "Ct = lowpass(C, H)\n", | |
| "rec = np.fft.ifft(Ct).real * Ct.shape[0]\n", | |
| "t = np.arange(0, 366)\n", | |
| "fig, ax = plt.subplots(facecolor=\"white\", figsize=(8, 6))\n", | |
| "ax.plot(t, temp2004, \"-\", label = \"temp2004\")\n", | |
| "ax.plot(t, rec, \"-\", label = f\"rec(H={H})\")\n", | |
| "ax.legend()\n", | |
| "ax.set_xlim([0, 365]) # 横軸の範囲\n", | |
| "ax.set_ylim([-5, 35]) # 縦軸の範囲\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "l3qsveG9ijtS", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "### フーリエ変換されたデータを入手して，逆変換してみよう" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "SJtdyiypa6MX", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "#### データの入手\n", | |
| "\n", | |
| "signal1, signal2 という変数に，それぞれとある音の信号を離散フーリエ変換して得たフーリエ係数値を格納する．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "h65z8eG--ndY", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# 音響データのファイルを入手する\n", | |
| "! wget -nc https://www-tlab.math.ryukoku.ac.jp/~takataka/course/PIP/hw08.npz\n", | |
| "! ls -l\n", | |
| "# ファイルの中身を取り出す\n", | |
| "hw08 = np.load('hw08.npz')\n", | |
| "#signal1, signal2, signal3 = hw08['signal1'], hw08['signal2'], hw08['signal3'] # 3つの音響データをフーリエ変換した値\n", | |
| "signal1, signal2 = hw08['signal1'], hw08['signal2'] # 2つの音響データをフーリエ変換した値\n", | |
| "signal1 /= 5 # 1/5 倍 （少し音を小さくする）\n", | |
| "signal2 /= 5\n", | |
| "# これらの音響データの標本化周波数\n", | |
| "framerate = 44100" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "-mcS3qr2a-2U", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "#### 入手したデータから振幅スペクトルを求めてみよう\n", | |
| "\n", | |
| "signal1 と signal2 の振幅スペクトルをグラフに描いてみる．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "6TamJCsf_eIw", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# signal1 と signal2 の振幅スペクトルを求める\n", | |
| "amp1 = np.abs(signal1)\n", | |
| "amp2 = np.abs(signal2)\n", | |
| "# それらの振幅スペクトルをグラフに描く\n", | |
| "N = amp1.shape[0]\n", | |
| "freq = np.arange(N) * framerate / N # 横軸の値．単位は Hz\n", | |
| "fig, ax = plt.subplots(facecolor=\"white\", figsize=(8, 6))\n", | |
| "ax.plot(freq, amp1, \"-\", label = \"amp1\")\n", | |
| "ax.plot(freq, amp2, \"-\", label = \"amp2\")\n", | |
| "ax.legend()\n", | |
| "ax.set_xlim([0, 8000]) # 横軸の範囲\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "riwfCVBSXMOE", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "#### 逆変換して元の信号に戻してみる\n", | |
| "\n", | |
| "signal1, signal2 をIDFTしたものを f1, f2 へ" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "yAgnVskqXRm5", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# signal1 と signal2 を逆変換して元の信号を求める\n", | |
| "f1 = np.fft.ifft(signal1).real * signal1.shape[0]\n", | |
| "f2 = np.fft.ifft(signal2).real * signal2.shape[0]" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "LVPooOrcbJj4", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "#### 波形を描いてみる\n", | |
| "\n", | |
| "f1, f2 のグラフを描こう" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "zWx6pII1VmJl", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# それらの最初の 0.01 秒分をグラフに描く\n", | |
| "t = np.arange(0, 0.01, 1.0/framerate)\n", | |
| "fig, ax = plt.subplots(facecolor=\"white\", figsize=(8, 6))\n", | |
| "ax.plot(t, f1[:t.shape[0]], \"-\", label = \"signal1\")\n", | |
| "ax.plot(t, f2[:t.shape[0]], \"-\", label = \"signal2\")\n", | |
| "ax.legend()\n", | |
| "ax.set_xlim([0, 0.01]) # 横軸の範囲\n", | |
| "ax.set_ylim([-32768, 32768]) # 縦軸の範囲（16bit量子化で表せる範囲全体）\n", | |
| "plt.show()" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "oYDbylYhXFpu", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "上記の振幅スペクトルおよび波形のグラフからすると，signal1 と signal2 はどちらの方が高い音に聞こえるはずか考えよう．" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "Fe4-EbFUYwFv", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "# WAVE 形式のファイルとして保存\n", | |
| "wavfile.write('signal1.wav', framerate, f1.astype(np.int16))\n", | |
| "wavfile.write('signal2.wav', framerate, f2.astype(np.int16))\n", | |
| "\n", | |
| "# 実際に音を聞いてみる場合は，以下の 0 を 1 に\n", | |
| "if 1 == 0:\n", | |
| " import os\n", | |
| " if os.path.exists('signal1.wav'):\n", | |
| " files.download('signal1.wav')\n", | |
| " if os.path.exists('signal2.wav'):\n", | |
| " files.download('signal2.wav')" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "metadata": { | |
| "id": "9zBtlpwBcJX-", | |
| "colab_type": "code", | |
| "colab": {} | |
| }, | |
| "source": [ | |
| "" | |
| ], | |
| "execution_count": 0, | |
| "outputs": [] | |
| } | |
| ] | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment