from numpy import random, array, inner, sqrt import matplotlib.pyplot as plt a, b = random.normal(0, 1, 2) e = array([b, -a]) / sqrt(a ** 2 + b ** 2) for n in range(100): v = random.uniform(-2, 2, 2) w = inner(e, v) * e plt.plot([v[0], w[0]], [v[1], w[1]]) plt.plot([v[0], w[0]], [v[1], w[1]], marker='.') plt.axis('scaled'), plt.xlim(-2, 2), plt.ylim(-2, 2), plt.show()

from numpy import array, vdot, sqrt def proj(x, E, inner=vdot): return sum([inner(e, x) * e for e in E]) 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__': A = [[1, 2, 3], [2, 3, 4], [3, 4, 5]] E = gram_schmidt(A) for n, e in enumerate(E): print(f'e{n+1} = {e}') print(array([[vdot(e1, e2) for e2 in E] for e1 in E]))

from vpython import proj, curve, vec, hat, arrow, label, box def proj2d(x, e): return x - proj(x, e) def draw_fig(c): curve(pos=c), curve(pos=[c[1], c[6]]) curve(pos=[c[2], c[7]]), curve(pos=[c[3], c[8]]) 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) for ax, lbl in [(x, 'x'), (y, 'y'), (z, 'z')]: curve(pos=[-5 * ax, 10 * ax], color=vec(1, 1, 1) - ax) label(pos=10 * ax, text=lbl) curve(pos=[proj2d(-5 * ax, e), proj2d(10 * ax, e)], color=ax) label(pos=proj2d(10 * ax, e), text=f"{lbl}'") c0 = [o, x, x + y, y, o, z, x + z, x + y + z, y + z, z] c1 = [u + 2 * v for v in c0] c2 = [proj2d(v, e) for v in c1] draw_fig(c1), draw_fig(c2)

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) plt.plot(A[:, 0], A[:, 1], color=color) def draw_fig(c): curve(pos=c), curve(pos=[c[1], c[6]]) curve(pos=[c[2], c[7]]), curve(pos=[c[3], c[8]]) o, e, u = array([0, 0, 0]), array([1, 2, 3]), array([5, 5, 5]) x, y, z = array([1, 0, 0]), array([0, 1, 0]), array([0, 0, 1]) E = array(gram_schmidt([e, x, y, z])[1:]) for v, t in [(x, "x'"), (y, "y'"), (z, "z'")]: curve(pos=[E.dot(-5 * v), E.dot(10 * v)], color=v) plt.text(*E.dot(10 * v), t, fontsize=24) c0 = [o, x, x + y, y, o, z, x + z, x + y + z, y + z, z] c1 = [u + 2 * v for v in c0] c2 = [E.dot(v) for v in c1] draw_fig(c2), plt.axis('equal'), plt.show()

from vpython import * from gram_schmidt2 import gram_schmidt canvas(background=color.white) x, y, z, v = [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 2, 3] E = gram_schmidt([v, x, y, z]) x, y, z, v = vec(*x), vec(*y), vec(*z), vec(*v) e0, e1, e2 = vec(*E[0]), vec(*E[1]), vec(*E[2]) for e in [x, y, z]: curve(pos=[e * (-5), e * 10], color=e) arrow(axis=v, color=color.yellow) with open('lena.txt', 'r') as fd: Data = eval(fd.read()) for x, y in Data: points(pos=[10 * x * e1 + 10 * y * e2], color=color.black, radius=2)

from numpy import array, sqrt def integral(f, D): N = len(D) - 1 w = (D[-1] - D[0]) / N x = array([(D[n] + D[n + 1]) / 2 for n in range(N)]) 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)}') print(f'||f_1||_2 = {norm["L2"](lambda x: x, D)}')

<sin|cos> = 4.5760562204146845e-17 ||f_1||_2 = 3.214875266047432

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) E = gram_schmidt([x**0, x**1, x**2, x**3]) 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()

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) FK = [f(k, t) for k in range(-K, K + 1)] 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): ax[0].plot(t, f(k, t)) for k in range(4): ax[1].plot(t, f(k, t)) plt.show()

from numpy import arange, cumsum, sqrt from numpy.random import normal from trigonometric import lowpass #from fourier 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) for k, K in enumerate([0, 1, 2, 4, 8, 16, 32, 64, 128]): i, j = divmod(k, 3) gK = lowpass(K, t, g) ax[i][j].plot(t, g), ax[i][j].plot(t, gK) ax[i][j].text(0.1, 0.5, f'K = {K}', fontsize = 20) ax[i][j].set_xlim(0, 1) plt.show()

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) FK = [f(k, t) for k in range(-K, K + 1)] 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) x = [z.real for z in f(k, t)] y = [z.imag for z in f(k, t)] fig = plt.figure() ax = Axes3D(fig) ax.scatter(t, x, y, s=1) ax.plot(t, x, y) plt.show()

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) x = array([(D[n] + D[n + 1]) / 2 for n in range(m)]) 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, } A = [x ** n for n in range(6)] E = gram_schmidt(A, inner=inner['Ledendre']) for e in E: plt.plot(x, e) plt.show()

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) x2 = x1[1:-1] 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()

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)), } dom, weight = D['Ledendre'] def inner(expr1, expr2): return integrate(f'({expr1}) * ({expr2}) * ({weight})', dom) def gram_schmidt(A): E = [] while A != []: a = A.pop(0) b = a - sum([inner(e, a) * e for e in E]) c = sqrt(inner(b, b)) E.append(b / c) return E E = gram_schmidt([1, x, x**2, x**3]) for n, e in enumerate(E): print(f'e{n}(x) = {e}')

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

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)}')

e0(x) = 1 e1(x) = x e2(x) = 3*x**2/2 - 1/2 e3(x) = 5*x**3/2 - 3*x/2

from numpy import array, sin, cos, pi, inf from numpy.linalg import norm import matplotlib.pyplot as plt def A(t): return array([[cos(t), -sin(t)], [sin(t), cos(t)]]) x = array([1, 0]) 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: ax[0].plot(p[0], p[1], marker='.', color='black') ax[1].plot(M, L1), plt.text(1, L1[0], 'L1', fontsize=16) ax[1].plot(M, L2), plt.text(1, L2[0], 'L2', fontsize=16) ax[1].plot(M, Loo), plt.text(1, Loo[0], 'Loo', fontsize=16) ax[1].set_xlim(0, 20) plt.show()

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: r1, r2 = rng[0] * self.tmax, rng[1] * self.tmax R = arange(int(r1), int(r2)) return R / self.tmax, spectrum[R] def lowpass(self, K): k = int(K * self.tmax) U = self.fft.copy() U[range(k + 1, self.len - k)] = 0 V = fft.ifft(U).real Data = (V * 32768).astype('int16') wav.write(f'{self.file}{K}.wav', self.rate, Data) return V

import matplotlib.pyplot as plt from sound import Sound sound1, sound2 = Sound('CEG'), Sound('mono') fig, ax = plt.subplots(1, 2, figsize=(15, 5)) ax[0].plot(*sound1.power_spectrum((-1000, 1000))) ax[1].plot(*sound2.power_spectrum()) plt.show()

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) ax[0].plot(X, Y), ax[0].plot(X, Y3000) ax[0].set_ylim(-1, 1) ax[1].plot(X, Y), ax[1].plot(X, Y3000) ax[1].set_xlim(0.2, 0.21), ax[1].set_ylim(-1, 1) plt.show()