6. 内積とフーリエ展開

1次元直交射影

code: proj.py

from numpy import random, array, inner, sqrt

import matplotlib.pyplot as plt

a, b = random.normal(0, 1, 2)

for n in range(100):

v = random.uniform(-2, 2, 2)

w = inner(e, v) * e

plt.axis('scaled'), plt.xlim(-2, 2), plt.ylim(-2, 2), plt.show()

https://gyazo.com/f27e30b65bfab0825c7ca83caf13f64a

グラム・シュミットの直交化法

code: gram_schmidt.py

from numpy import array, vdot, sqrt

def proj(x, E, inner=vdot):

def gram_schmidt(A, inner=vdot):

E = []

while A != []:

a = array(A.pop(0))

b = a - proj(a, E, inner)

c = sqrt(inner(b, b))

if c >= 1.0e-15:

E.append(b / c)

return E

if __name__ == '__main__':

E = gram_schmidt(A)

for n, e in enumerate(E):

print(f'e{n+1} = {e}')

2次元直交射影

code: proj2d.py

from vpython import proj, curve, vec, hat, arrow, label, box

def proj2d(x, e): return x - proj(x, e)

def draw_fig(c):

o, e, u = vec(0, 0, 0), hat(vec(1, 2, 3)), vec(5, 5, 5)

x, y, z = vec(1, 0, 0), vec(0, 1, 0), vec(0, 0, 1)

box(axis=e, up=proj2d(y, e),

width=20, height=20, length=0.1, opacity=0.5)

arrow(axis=3 * e)

label(pos=10 * ax, text=lbl)

label(pos=proj2d(10 * ax, e), text=f"{lbl}'")

draw_fig(c1), draw_fig(c2)

https://gyazo.com/24978fd2af3ebfdb570e6939140f2190

code: screen.py

from numpy import array

import matplotlib.pyplot as plt

from gram_schmidt import gram_schmidt

def curve(pos, color=(0, 0, 0)):

A = array(pos)

def draw_fig(c):

plt.text(*E.dot(10 * v), t, fontsize=24)

draw_fig(c2), plt.axis('equal'), plt.show()

https://gyazo.com/a72b8e5edc1c6fbca70d24f60b7c2f3e

問6.4.

code: lena5.py

from vpython import *

from gram_schmidt2 import gram_schmidt

canvas(background=color.white)

x, y, z, v = vec(*x), vec(*y), vec(*z), vec(*v)

arrow(axis=v, color=color.yellow)

with open('lena.txt', 'r') as fd:

Data = eval(fd.read())

for x, y in Data:

color=color.black, radius=2)

https://gyazo.com/660bac51da718b8039373a36dca30adc

関数空間

code: integral.py

from numpy import array, sqrt

def integral(f, D):

N = len(D) - 1

return sum(f(x)) * w

def inner(f, g, D):

return integral(lambda x: f(x).conj() * g(x), D)

norm = {

'L1': lambda f, D: integral(lambda x: abs(f(x)), D),

'L2': lambda f, D: sqrt(inner(f, f, D)),

'Loo': lambda f, D: max(abs(f(D))),

}

if __name__ == '__main__':

from numpy import linspace, pi, sin, cos

D = linspace(0, pi, 1001)

print(f'<sin|cos> = {inner(sin, cos, D)}')

実行結果

code: python

<sin|cos> = 4.5760562204146845e-17

||f_1||_2 = 3.214875266047432

最小2乗法

code: lstsqr.py

from numpy import linspace, cumsum, vdot, sort

from numpy.random import uniform, normal

from gram_schmidt import gram_schmidt, proj

import matplotlib.pyplot as plt

n = 20

x = sort(uniform(-1, 1, n))

z = 4 * x**3 - 3 * x

sigma = 0.2

y = z + normal(0, sigma, n)

y0 = proj(z, E)

plt.figure(figsize=(15,5))

plt.errorbar(x, z, yerr=sigma, fmt='ro')

plt.plot(x, y0, color='g'), plt.plot(x, y, color='b'), plt.show()

https://gyazo.com/3a14e792d27381cc00c6b344909a3148

三角級数

code: trigonometric.py

from numpy import inner, pi, sin, cos, sqrt, ones

from gram_schmidt import proj

def f(k, t):

if k < 0:

return sin(2 * k * pi * t) * sqrt(2)

elif k == 0:

return ones(len(t))

elif k > 0:

return cos(2 * k * pi * t) * sqrt(2)

def lowpass(K, t, g):

n = len(t)

return proj(g, FK, inner=lambda x, y: inner(x, y) / n)

if __name__ == '__main__':

from numpy import arange

import matplotlib.pyplot as plt

t = arange(0, 1, 1 / 1000)

fig, ax = plt.subplots(1, 2, figsize=(15, 5))

for k in range(-3, 0):

for k in range(4):

plt.show()

https://gyazo.com/335f25cfffe212d421da9078211ff0c6

ローパスフィルタ

code: brown.py

from numpy import arange, cumsum, sqrt

from numpy.random import normal

from trigonometric import lowpass

import matplotlib.pyplot as plt

n = 1000

dt = 1 / n

t = arange(0, 1, dt)

g = cumsum(normal(0, sqrt(dt), n))

fig, ax = plt.subplots(3, 3, figsize=(16, 8), dpi=100)

i, j = divmod(k, 3)

gK = lowpass(K, t, g)

plt.show()

https://gyazo.com/f39a5bd9396d4052599f66763f2765a5

フーリエ級数

code: fourier.py

from numpy import arange, exp, pi, vdot

from gram_schmidt import proj

def f(k, t): return exp(2j * pi * k * t)

def lowpass(K, t, z):

dt = 1 / len(t)

return proj(z, FK, inner=lambda x, y: vdot(x, y) * dt)

if __name__ == '__main__':

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

k = 2

t = arange(0, 1, 1 / 1000)

fig = plt.figure()

ax = Axes3D(fig)

ax.scatter(t, x, y, s=1)

ax.plot(t, x, y)

plt.show()

https://gyazo.com/b8d402b41170eaf18cfbeece6b802615

直交多項式

code: poly_np1.py

from numpy import array, linspace, sqrt, ones, pi

import matplotlib.pyplot as plt

from gram_schmidt import gram_schmidt

m = 10000

D = linspace(-1, 1, m + 1)

inner = {

'Ledendre': lambda f, g: f.dot(g) * 2 / m,

'Chebyshev1': lambda f, g: f.dot(g / sqrt(1 - x**2)) * 2 / m,

'Chebyshev2': lambda f, g: f.dot(g * sqrt(1 - x**2)) * 2 / m,

}

for e in E:

plt.plot(x, e)

plt.show()

https://gyazo.com/a669ce2b984610dddf028c6a09049c35

code: poly_np2.py

from numpy import linspace, exp, sqrt

from numpy.polynomial.chebyshev import Chebyshev

from numpy.polynomial.legendre import Legendre

from numpy.polynomial.laguerre import Laguerre

from numpy.polynomial.hermite import Hermite

import matplotlib.pyplot as plt

x1 = linspace(-1, 1, 1001)

x3 = linspace(0, 10, 1001)

x4 = linspace(-3, 3, 1001)

f, x, w = Legendre, x1, 1

# f, x, w = Chebyshev, x2, 1

# f, x, w = Chebyshev, x2, 1/sqrt(1 - x2**2)

# f, x, w = Laguerre, x3, 1

# f, x, w = Laguerre, x3, exp(-x3)

# f, x, w = Hermite, x4, 1

# f, x, w = Hermite, x4, exp(-x4**2)

for n in range(6):

e = f.basis(n)(x)

plt.plot(x, e * w)

plt.show()

https://gyazo.com/616ab4f27a0d337c6f2bcd7e449c8143

code: poly_sp1.py

from sympy import integrate, sqrt, exp, oo

from sympy.abc import x

D = {

'Ledendre': ((x, -1, 1), 1),

'Chebyshev1': ((x, -1, 1), 1 / sqrt(1 - x**2)),

'Chebyshev2': ((x, -1, 1), sqrt(1 - x**2)),

'Laguerre': ((x, 0, oo), exp(-x)),

'Hermite': ((x, -oo, oo), exp(-x**2)),

}

def inner(expr1, expr2):

return integrate(f'({expr1}) * ({expr2}) * ({weight})', dom)

def gram_schmidt(A):

E = []

while A != []:

a = A.pop(0)

c = sqrt(inner(b, b))

E.append(b / c)

return E

for n, e in enumerate(E):

print(f'e{n}(x) = {e}')

実行結果

code: python

e0(x) = sqrt(2)/2

e1(x) = sqrt(6)*x/2

e2(x) = 3*sqrt(10)*(x**2 - 1/3)/4

e3(x) = 5*sqrt(14)*(x**3 - 3*x/5)/4

code: poly_sp2.py

from sympy.polys.orthopolys import (

legendre_poly,

chebyshevt_poly,

chebyshevu_poly,

laguerre_poly,

hermite_poly,

)

from sympy.abc import x

e = legendre_poly

for n in range(4):

print(f'e{n}(x) = {e(n, x)}')

実行結果

code: python

e0(x) = 1

e1(x) = x

e2(x) = 3*x**2/2 - 1/2

e3(x) = 5*x**3/2 - 3*x/2

ベクトル列の収束

code: limit.py

from numpy import array, sin, cos, pi, inf

from numpy.linalg import norm

import matplotlib.pyplot as plt

def A(t):

P, L1, L2, Loo = [], [], [], []

M = range(1, 100)

for m in M:

xm = A(pi / 2 / m).dot(x)

P.append(xm)

L1.append(norm(x - xm, ord=1))

L2.append(norm(x - xm))

Loo.append(norm(x - xm, ord=inf))

fig, ax = plt.subplots(1, 2, figsize=(10, 5))

for p in P:

plt.show()

https://gyazo.com/8ebe0a7862289d68842379dea6c3adeb

フーリエ解析

code: sound.py

from numpy import arange, fft

import scipy.io.wavfile as wav

class Sound:

def __init__(self, wavfile):

self.file = wavfile

self.rate, Data = wav.read(f'{wavfile}.wav')

dt = 1 / self.rate

self.len = len(Data)

self.tmax = self.len / self.rate

self.time = arange(0, self.tmax, dt)

self.data = Data.astype('float') / 32768

self.fft = fft.fft(self.data)

def power_spectrum(self, rng=None):

spectrum = abs(self.fft) ** 2

if rng is None:

r1, r2 = -self.len / 2, self.len / 2

else:

R = arange(int(r1), int(r2))

def lowpass(self, K):

k = int(K * self.tmax)

U = self.fft.copy()

V = fft.ifft(U).real

Data = (V * 32768).astype('int16')

wav.write(f'{self.file}{K}.wav', self.rate, Data)

return V

code: spectrum.py

import matplotlib.pyplot as plt

from sound import Sound

sound1, sound2 = Sound('CEG'), Sound('mono')

fig, ax = plt.subplots(1, 2, figsize=(15, 5))

plt.show()

https://gyazo.com/da0b6fa20fd1e6c36d29739c006c54ba

code: lowpass.py

import matplotlib.pyplot as plt

from sound import Sound

sound = Sound('mono')

X, Y = sound.time, sound.data

Y3000 = sound.lowpass(3000)

fig, ax = plt.subplots(1, 2, figsize=(10, 5), dpi=100)

plt.show()

https://gyazo.com/f4944dd03f85d7cb2f5001098fabdbc9