L26 서술 파트는 세 걸음이었다. 복소수 위에서 길이를 되찾고, 푸리에 행렬이 유니타리임을 보이고, 그 행렬을 쪼개 으로 줄였다.
이 노트북에서는 A.T 와 A.conj().T 가 어떻게 다른 답을 내는지 눈으로 보고,
를 직접 확인하며, 실제 신호에서 주파수를 뽑아낸다.
마지막으로 FFT를 재귀 함수 열 줄로 짜서 numpy.fft 와 대조하고 속도를 잰다.
| 서술 파트의 내용 | 여기서 확인하는 방법 |
|---|---|
| 가 무너진다 | 로 0이 나오는 것 |
| 켤레전치가 고친다 | np.vdot 과 A.conj().T |
| 파이썬 함정 | A.T 는 오류 없이 틀린 답을 낸다 |
| 에르미트 | 고윳값이 실수, 대각도 실수 |
| 유니타리 | 길이 보존, |
| , | 단위원 위의 점들 |
| 푸리에 행렬 | 를 코드로, np.fft.fft 와 일치 |
| 여러 에서 오차 | |
| 역변환이 공짜 | |
| 열이 순수한 음 | 슬라이더로 주파수를 올려 가며 |
| 신호 분석 | 섞인 두 진동을 도로 갈라낸다 |
| FFT 인수분해 | 세 조각의 곱이 과 같은지 |
| 재귀 구현과 속도 측정 |
0. 준비¶
import time
import numpy as np
import plotly.graph_objects as go
from linalg_viz import COLORS, layout2d, show_matrix, slider_figure
np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(26)
print("numpy", np.__version__)numpy 2.5.2
1. 길이가 무너지고, 켤레가 되살린다¶
x = np.array([1, 1j])
print("x =", x)
print()
print("x @ x =", x @ x, " <- 0. 영벡터가 아닌데 길이가 0")
print("x.conj() @ x =", (x.conj() @ x).real, " <- 2. 이것이 옳다")
print("np.vdot(x, x) =", np.vdot(x, x).real, " <- vdot 은 첫 인자를 켤레 취한다")
print("np.abs(x)**2 의 합 =", (np.abs(x) ** 2).sum())
print()
print("길이 :", np.linalg.norm(x), " = sqrt(2) =", np.sqrt(2))x = [1.+0.j 0.+1.j]
x @ x = 0j <- 0. 영벡터가 아닌데 길이가 0
x.conj() @ x = 2.0 <- 2. 이것이 옳다
np.vdot(x, x) = 2.0 <- vdot 은 첫 인자를 켤레 취한다
np.abs(x)**2 의 합 = 2.0
길이 : 1.4142135623730951 = sqrt(2) = 1.4142135623730951
A = np.array([[1 + 2j, 3 - 1j],
[0 + 1j, 2 + 0j]])
print("A =\n", A)
print()
print("A.T @ A <- 틀림\n", A.T @ A)
print()
print("A.conj().T @ A <- 맞음\n", A.conj().T @ A)
print()
print("맞은 쪽의 대각 :", np.diag(A.conj().T @ A), " <- 실수. 각 열의 길이 제곱")
print(" 실제 열의 길이 제곱 :", (np.abs(A) ** 2).sum(axis=0), " <- 같다")
print("틀린 쪽의 대각 :", np.diag(A.T @ A), " <- 실수도 아니다")
print()
print("맞은 쪽은 자기 켤레전치와 같은가 :",
np.allclose(A.conj().T @ A, (A.conj().T @ A).conj().T))
print("틀린 쪽은 어떤가 :",
np.allclose(A.T @ A, (A.T @ A).conj().T))A =
[[1.+2.j 3.-1.j]
[0.+1.j 2.+0.j]]
A.T @ A <- 틀림
[[-4.+4.j 5.+7.j]
[ 5.+7.j 12.-6.j]]
A.conj().T @ A <- 맞음
[[ 6.+0.j 1.-9.j]
[ 1.+9.j 14.+0.j]]
맞은 쪽의 대각 : [ 6.+0.j 14.+0.j] <- 실수. 각 열의 길이 제곱
실제 열의 길이 제곱 : [ 6. 14.] <- 같다
틀린 쪽의 대각 : [-4.+4.j 12.-6.j] <- 실수도 아니다
맞은 쪽은 자기 켤레전치와 같은가 : True
틀린 쪽은 어떤가 : False
2. 에르미트와 유니타리¶
def 에르미트인가(A, 눈감아=1e-12):
"""A^H = A 인지 본다."""
A = np.asarray(A)
return bool(np.abs(A - A.conj().T).max() < 눈감아)
def 유니타리인가(U, 눈감아=1e-10):
"""U^H U = I 인지 본다."""
U = np.asarray(U)
return bool(np.abs(U.conj().T @ U - np.eye(U.shape[1])).max() < 눈감아)H = np.array([[2, 1 - 1j],
[1 + 1j, 3]])
print("H =\n", H)
print("에르미트인가 :", 에르미트인가(H))
print("대각 :", np.diag(H), " <- 실수여야 한다")
print("고윳값 :", np.linalg.eigvalsh(H), " <- 실수. 기대한 1, 4")
print(" 대각합", np.trace(H).real, "= 1 + 4, 행렬식",
round(float(np.linalg.det(H).real), 10), "= 1 * 4")
print()
값, V = np.linalg.eigh(H)
print("고유벡터가 직교인가 :", 유니타리인가(V))
print(" V^H V - I 최대 :", np.abs(V.conj().T @ V - np.eye(2)).max())H =
[[2.+0.j 1.-1.j]
[1.+1.j 3.+0.j]]
에르미트인가 : True
대각 : [2.+0.j 3.+0.j] <- 실수여야 한다
고윳값 : [1. 4.] <- 실수. 기대한 1, 4
대각합 5.0 = 1 + 4, 행렬식 4.0 = 1 * 4
고유벡터가 직교인가 : True
V^H V - I 최대 : 3.3306690738754696e-16
대각에 허수를 넣으면 에르미트가 아니다.
나쁨 = np.array([[1j, 1], [1, 0]])
print("[[i,1],[1,0]] 이 에르미트인가 :", 에르미트인가(나쁨))
print(" 고윳값 :", np.linalg.eigvals(나쁨), " <- 실수가 아니다")[[i,1],[1,0]] 이 에르미트인가 : False
고윳값 : [ 0.866+0.5j -0.866+0.5j] <- 실수가 아니다
# 무작위 에르미트 행렬 300개
최대허수 = 0.0
for _ in range(300):
n = int(rng.integers(2, 7))
B = rng.normal(size=(n, n)) + 1j * rng.normal(size=(n, n))
A = B + B.conj().T # 반드시 에르미트
assert 에르미트인가(A)
최대허수 = max(최대허수, float(np.abs(np.linalg.eigvals(A).imag).max()))
print(f"에르미트 300개의 고윳값 허수부 최댓값 : {최대허수:.3e} -> 전부 실수")에르미트 300개의 고윳값 허수부 최댓값 : 1.412e-15 -> 전부 실수
3. 단위원 위의 개 점¶
for N in (4, 8, 16, 64):
w = np.exp(-2j * np.pi / N)
print(f"N={N:>3} : w^N = {w**N:>22.12f} w^(N/2) = {w**(N//2):>22.12f}")N= 4 : w^N = 1.000000000000+0.000000000000j w^(N/2) = -1.000000000000-0.000000000000j
N= 8 : w^N = 1.000000000000+0.000000000000j w^(N/2) = -1.000000000000-0.000000000000j
N= 16 : w^N = 1.000000000000-0.000000000000j w^(N/2) = -1.000000000000+0.000000000000j
N= 64 : w^N = 1.000000000000+0.000000000000j w^(N/2) = -1.000000000000-0.000000000000j
N = 8
w = np.exp(-2j * np.pi / N)
거듭 = w ** np.arange(N)
원각 = np.linspace(0, 2 * np.pi, 300)
색 = ["#2ca02c" if k == 0 else ("#d62728" if k == N // 2 else COLORS["input"])
for k in range(N)]
go.Figure(
data=[go.Scatter(x=np.cos(원각), y=np.sin(원각), mode="lines",
line=dict(color="#cccccc", dash="dot"), name="단위원"),
go.Scatter(x=거듭.real, y=거듭.imag, mode="markers+text",
marker=dict(size=14, color=색),
text=[f"w^{k}" for k in range(N)],
textposition="top center", name="w 의 거듭제곱")],
layout=layout2d(f"w = exp(-2 pi i / {N}) 의 거듭제곱", extent=1.6))Loading...
print("여덟 점을 전부 더하면 :", 거듭.sum(), " <- 0. 골고루 흩어져 상쇄된다")
print("이것이 5절 직교성의 씨앗이다.")여덟 점을 전부 더하면 : (-1.1102230246251565e-15-4.440892098500626e-16j) <- 0. 골고루 흩어져 상쇄된다
이것이 5절 직교성의 씨앗이다.
4. 푸리에 행렬¶
def 푸리에행렬(N):
"""(F_N)_jk = w^(jk), w = exp(-2 pi i / N). numpy.fft 와 같은 규약."""
j, k = np.meshgrid(np.arange(N), np.arange(N), indexing="ij")
return np.exp(-2j * np.pi * j * k / N)F4 = 푸리에행렬(4)
print("F_4 (반올림) =")
print(np.round(F4, 10))
손 = np.array([[1, 1, 1, 1],
[1, -1j, -1, 1j],
[1, -1, 1, -1],
[1, 1j, -1, -1j]], dtype=complex)
print()
print("서술 파트에서 손으로 쓴 것과 같은가 :", np.allclose(F4, 손))
print("0행이 전부 1 인가 :", np.allclose(F4[0], 1))
print("0열이 전부 1 인가 :", np.allclose(F4[:, 0], 1))F_4 (반올림) =
[[ 1.+0.j 1.+0.j 1.+0.j 1.+0.j]
[ 1.+0.j 0.-1.j -1.-0.j -0.+1.j]
[ 1.+0.j -1.-0.j 1.+0.j -1.-0.j]
[ 1.+0.j -0.+1.j -1.-0.j 0.-1.j]]
서술 파트에서 손으로 쓴 것과 같은가 : True
0행이 전부 1 인가 : True
0열이 전부 1 인가 : True
각 열은 하나의 순수한 음이다¶
열을 에 대해 그려 보자. 가 커질수록 빨리 진동한다.
N = 32
F = 푸리에행렬(N)
j = np.arange(N)
프레임 = []
for k in range(N):
열 = F[:, k]
프레임.append([
go.Scatter(x=j, y=열.real, mode="lines+markers",
line=dict(color=COLORS["input"], width=2.5),
marker=dict(size=5), name="실수부 (cos)"),
go.Scatter(x=j, y=열.imag, mode="lines+markers",
line=dict(color=COLORS["output"], width=2.5, dash="dash"),
marker=dict(size=5), name="허수부 (sin)"),
])
배치 = dict(title=dict(text=f"F_{N} 의 k 번째 열"),
xaxis=dict(title=dict(text="j (시각)"), range=[-1, N]),
yaxis=dict(title=dict(text="값"), range=[-1.25, 1.25]),
height=430, margin=dict(l=70, r=20, t=60, b=50))
slider_figure(프레임, list(range(N)), 배치, prefix="k = ", initial=3)Loading...
이면 아예 진동하지 않는 상수이다. 를 올리면 한 화면 안에 들어가는 진동의 수가 정확히 개가 된다. ()에서 가장 빠르고, 그보다 더 올리면 오히려 느려 보인다. 와 가 같은 속도로 반대 방향을 돌기 때문이다.
5. ¶
서술 파트의 증명은 등비급수 한 줄이었다. 숫자로 확인하자.
print(f"{'N':>6}{'|F^H F - N I| 최대':>24}")
for N in (4, 8, 16, 64, 256):
F = 푸리에행렬(N)
print(f"{N:>6}{np.abs(F.conj().T @ F - N * np.eye(N)).max():>24.2e}") N |F^H F - N I| 최대
4 4.68e-16
8 5.85e-15
16 3.44e-14
64 3.74e-13
256 4.47e-12
N = 8
F = 푸리에행렬(N)
U = F / np.sqrt(N)
print("F/sqrt(N) 이 유니타리인가 :", 유니타리인가(U))
print("고윳값의 절댓값 :", np.round(np.abs(np.linalg.eigvals(U)), 10))
v = rng.standard_normal(N) + 1j * rng.standard_normal(N)
print("길이 보존 :", round(float(np.linalg.norm(U @ v)), 10),
"vs", round(float(np.linalg.norm(v)), 10))F/sqrt(N) 이 유니타리인가 : True
고윳값의 절댓값 : [1. 1. 1. 1. 1. 1. 1. 1.]
길이 보존 : 3.6530216321 vs 3.6530216321
역변환이 공짜¶
F역_공식 = F.conj().T / N
F역_소거 = np.linalg.inv(F)
print("F^H/N 과 inv(F) 의 최대 차이 :", np.abs(F역_공식 - F역_소거).max())
print()
x = rng.standard_normal(N)
c = F @ x
print("x :", x)
print("F^H c/N :", (F역_공식 @ c).real)
print("돌아왔는가 :", np.allclose((F역_공식 @ c).real, x))
print()
print("np.fft.ifft 와도 같은가 :", np.allclose(F역_공식 @ c, np.fft.ifft(c)))F^H/N 과 inv(F) 의 최대 차이 : 3.546311116759026e-16
x : [ 0.2994 0.0713 -2.507 1.0607 -0.2628 0.4883 -0.1615 -0.0038]
F^H c/N : [ 0.2994 0.0713 -2.507 1.0607 -0.2628 0.4883 -0.1615 -0.0038]
돌아왔는가 : True
np.fft.ifft 와도 같은가 : True
print(f"{'N':>6}{'F@x 와 np.fft.fft 의 최대 오차':>32}")
for N in (8, 64, 256, 512):
x = rng.standard_normal(N)
print(f"{N:>6}{np.abs(푸리에행렬(N) @ x - np.fft.fft(x)).max():>32.2e}") N F@x 와 np.fft.fft 의 최대 오차
8 3.51e-15
64 3.53e-13
256 3.32e-12
512 8.27e-12
6. 섞인 진동을 도로 갈라낸다¶
여기가 푸리에 변환을 쓰는 이유이다. 두 개의 진동을 섞고 잡음까지 넣은 신호에서 원래의 두 주파수를 뽑아내 보자.
N = 256
j = np.arange(N)
신호 = (3.0 * np.sin(2 * np.pi * 5 * j / N)
+ 1.5 * np.sin(2 * np.pi * 17 * j / N)
+ 0.6 * rng.standard_normal(N))
c = 푸리에행렬(N) @ 신호
크기 = np.abs(c)
봉우리 = np.argsort(크기[: N // 2])[::-1][:4]
print("가장 큰 계수가 있는 k :", np.sort(봉우리))
print(" (원래 넣은 것은 5 와 17)")
print()
for k in np.sort(봉우리):
print(f" k = {k:>3} : |c_k| = {크기[k]:>8.1f}")가장 큰 계수가 있는 k : [ 5 17 39 63]
(원래 넣은 것은 5 와 17)
k = 5 : |c_k| = 380.0
k = 17 : |c_k| = 180.3
k = 39 : |c_k| = 18.6
k = 63 : |c_k| = 20.3
go.Figure(
data=[go.Scatter(x=j, y=신호, mode="lines",
line=dict(color=COLORS["input"], width=1.6))],
layout=go.Layout(title=dict(text="신호 자체 — 눈으로는 무엇이 섞였는지 알 수 없다"),
xaxis=dict(title=dict(text="j")),
yaxis=dict(title=dict(text="값")),
height=340, margin=dict(l=70, r=20, t=60, b=50)))Loading...
go.Figure(
data=[go.Bar(x=j[: N // 2], y=크기[: N // 2],
marker=dict(color=COLORS["output"]))],
layout=go.Layout(title=dict(text="|F x| — 5 와 17 에 봉우리가 선다"),
xaxis=dict(title=dict(text="k (주파수)"), range=[-1, 40]),
yaxis=dict(title=dict(text="크기")),
height=380, margin=dict(l=70, r=20, t=60, b=50)))Loading...
잡음을 걷어내는 것도 한 줄이다. 작은 계수를 0으로 만들고 역변환하면 된다.
걸러낸c = np.where(크기 > 100, c, 0)
복원 = (푸리에행렬(N).conj().T @ 걸러낸c / N).real
깨끗 = 3.0 * np.sin(2*np.pi*5*j/N) + 1.5 * np.sin(2*np.pi*17*j/N)
print("남긴 계수의 개수 :", int((크기 > 100).sum()), "/", N)
print("잡음 없는 원본과의 평균 오차 :", np.abs(복원 - 깨끗).mean())
print("거르기 전 신호의 평균 오차 :", np.abs(신호 - 깨끗).mean())
go.Figure(
data=[go.Scatter(x=j, y=신호, mode="lines", name="잡음 섞인 신호",
line=dict(color="#cccccc", width=1.4)),
go.Scatter(x=j, y=복원, mode="lines", name="큰 계수만 남기고 복원",
line=dict(color=COLORS["output"], width=3)),
go.Scatter(x=j, y=깨끗, mode="lines", name="잡음 없는 원본",
line=dict(color=COLORS["input"], width=2, dash="dash"))],
layout=go.Layout(title=dict(text="주파수 쪽에서 자르면 시간 쪽이 깨끗해진다"),
xaxis=dict(title=dict(text="j"), range=[0, 128]),
height=400, margin=dict(l=70, r=20, t=60, b=50)))남긴 계수의 개수 : 4 / 256
잡음 없는 원본과의 평균 오차 : 0.06509237122179179
거르기 전 신호의 평균 오차 : 0.5151701800534789
Loading...
def 세조각(n):
"""F_2n 을 쪼갠 세 조각을 돌려준다."""
N = 2 * n
w = np.exp(-2j * np.pi / N)
D = np.diag(w ** np.arange(n))
앞 = np.block([[np.eye(n), D], [np.eye(n), -D]])
가운데 = np.block([[푸리에행렬(n), np.zeros((n, n))],
[np.zeros((n, n)), 푸리에행렬(n)]])
P = np.zeros((N, N))
for 자리, k in enumerate(list(range(0, N, 2)) + list(range(1, N, 2))):
P[자리, k] = 1.0 # 짝수 먼저, 그다음 홀수
return 앞, 가운데, P앞, 가운데, P = 세조각(2)
print(show_matrix(P, "P (4x4). 짝수 자리를 위로"))
print("P 가 순열행렬인가 :", np.allclose(P @ P.T, np.eye(4)))
print("P (a,b,c,d) =", P @ np.array([0, 1, 2, 3.]), " <- 0,2 가 앞으로")
print()
print(show_matrix(np.round(앞.real, 6), "[[I, D], [I, -D]] 의 실수부"))P (4x4). 짝수 자리를 위로
[ 1 0 0 0 ]
[ 0 0 1 0 ]
[ 0 1 0 0 ]
[ 0 0 0 1 ]
P 가 순열행렬인가 : True
P (a,b,c,d) = [0. 2. 1. 3.] <- 0,2 가 앞으로
[[I, D], [I, -D]] 의 실수부
[ 1 0 1 0 ]
[ 0 1 0 0 ]
[ 1 0 -1 -0 ]
[ 0 1 -0 -0 ]
print(f"{'n':>5}{'N=2n':>7}{'|앞·가운데·P - F_N| 최대':>28}")
for n in (2, 4, 8, 16, 32):
앞, 가운데, P = 세조각(n)
print(f"{n:>5}{2*n:>7}{np.abs(앞 @ 가운데 @ P - 푸리에행렬(2*n)).max():>28.2e}") n N=2n |앞·가운데·P - F_N| 최대
2 4 3.67e-16
4 8 3.30e-15
8 16 1.57e-14
16 32 3.14e-14
32 64 6.27e-14
def 나의fft(신호):
"""분할정복 구조를 그대로. N 은 2의 거듭제곱이라고 가정한다."""
신호 = np.asarray(신호, dtype=complex)
N = 신호.size
if N == 1:
return 신호
e = 나의fft(신호[0::2]) # 짝수 자리
o = 나의fft(신호[1::2]) # 홀수 자리
회전 = np.exp(-2j * np.pi * np.arange(N // 2) / N) # D 의 대각
보정 = 회전 * o
return np.concatenate([e + 보정, e - 보정]) # +D 와 -Dfor N in (8, 64, 1024):
x = rng.standard_normal(N)
print(f"N={N:>5} : 나의fft 와 np.fft.fft 의 최대 오차 "
f"{np.abs(나의fft(x) - np.fft.fft(x)).max():.2e}")N= 8 : 나의fft 와 np.fft.fft 의 최대 오차 4.44e-16
N= 64 : 나의fft 와 np.fft.fft 의 최대 오차 9.06e-15
N= 1024 : 나의fft 와 np.fft.fft 의 최대 오차 5.52e-14
x = np.arange(64.0)
print("행렬 곱 vs 나의fft 최대 오차 :",
np.abs(푸리에행렬(64) @ x - 나의fft(x)).max())
print("FFT 는 근사가 아니라 같은 답이다.")행렬 곱 vs 나의fft 최대 오차 : 9.223403739270653e-12
FFT 는 근사가 아니라 같은 답이다.
8. 얼마나 빨라지는가¶
print(f"{'N':>10}{'N^2':>14}{'(N/2) log2 N':>16}{'몇 배':>12}")
for k in (10, 15, 20, 24):
N = 2 ** k
print(f"{N:>10,}{N*N:>14.2e}{N//2*k:>16.2e}{N*N/(N//2*k):>12,.0f}") N N^2 (N/2) log2 N 몇 배
1,024 1.05e+06 5.12e+03 205
32,768 1.07e+09 2.46e+05 4,369
1,048,576 1.10e+12 1.05e+07 104,858
16,777,216 2.81e+14 2.01e+08 1,398,101
크기목록 = [2 ** k for k in range(7, 13)]
느린, 빠른 = [], []
for N in 크기목록:
x = rng.standard_normal(N)
F = 푸리에행렬(N)
반복 = max(3, 200_000 // N) # 작은 N 은 여러 번 재야 한다
t0 = time.perf_counter()
for _ in range(반복):
F @ x
느린.append((time.perf_counter() - t0) / 반복)
t0 = time.perf_counter()
for _ in range(반복):
np.fft.fft(x)
빠른.append((time.perf_counter() - t0) / 반복)
print(f"{'N':>6}{'행렬곱 (초)':>16}{'FFT (초)':>14}{'몇 배':>10}")
for N, a, b in zip(크기목록, 느린, 빠른):
print(f"{N:>6}{a:>16.2e}{b:>14.2e}{a/b:>10.1f}") N 행렬곱 (초) FFT (초) 몇 배
128 8.31e-06 1.13e-05 0.7
256 2.76e-05 1.27e-05 2.2
512 1.04e-04 1.64e-05 6.3
1024 3.35e-04 2.29e-05 14.6
2048 1.70e-03 3.90e-05 43.6
4096 9.20e-03 7.36e-05 124.9
go.Figure(
data=[go.Scatter(x=크기목록, y=느린, mode="lines+markers",
name="행렬 곱 O(N^2)",
line=dict(color="#d62728", width=3)),
go.Scatter(x=크기목록, y=빠른, mode="lines+markers",
name="np.fft.fft O(N log N)",
line=dict(color="#2ca02c", width=3))],
layout=go.Layout(title=dict(text="같은 답, 다른 시간"),
xaxis=dict(title=dict(text="N"), type="log"),
yaxis=dict(title=dict(text="초"), type="log"),
height=420, margin=dict(l=70, r=20, t=60, b=50)))Loading...
로그-로그 그래프에서 기울기가 곧 지수이다. 이론값은 2와 1이다. 재 보자.
기울기 = lambda x, y: float(np.polyfit(np.log(x), np.log(y), 1)[0])
print(f"{'구간':>18}{'행렬곱':>10}")
print(f"{'전체':>18}{기울기(크기목록, 느린):>10.2f}")
print(f"{'N >= 1024 만':>18}{기울기(크기목록[3:], 느린[3:]):>10.2f} <- 이론값 2") 구간 행렬곱
전체 2.00
N >= 1024 만 2.39 <- 이론값 2
행렬 곱은 뒤쪽 구간에서 기울기가 2에 정확히 맞는다. 작은 에서 기울기가 낮게 나오는 것은 함수 호출과 배열 준비 같은 고정 비용이 실제 계산보다 크기 때문이다.
FFT는 이 범위에서 아예 기울기를 잴 수가 없다. 시간이 수십 마이크로초라 전부 고정 비용이기 때문이다. FFT는 메모리를 거의 쓰지 않으니 훨씬 큰 까지 갈 수 있다. 따로 재 보자.
큰크기 = [2 ** k for k in range(12, 23)]
큰시간 = []
for N in 큰크기:
x = rng.standard_normal(N)
np.fft.fft(x) # 첫 호출은 버린다
반복 = max(3, 2_000_000 // N)
t0 = time.perf_counter()
for _ in range(반복):
np.fft.fft(x)
큰시간.append((time.perf_counter() - t0) / 반복)
print(f"{'N':>12}{'초':>13}")
for N, t in zip(큰크기, 큰시간):
print(f"{N:>12,}{t:>13.2e}")
print()
print(f"전체 기울기 : {기울기(큰크기, 큰시간):.2f}")
print(f"N >= 65536 만 : {기울기(큰크기[4:], 큰시간[4:]):.2f} "
f"<- 이론값 1 을 조금 넘는다 (log 때문)") N 초
4,096 7.41e-05
8,192 1.43e-04
16,384 2.16e-04
32,768 4.71e-04
65,536 1.03e-03
131,072 1.95e-03
262,144 4.10e-03
524,288 1.47e-02
1,048,576 2.81e-02
2,097,152 6.08e-02
4,194,304 1.20e-01
전체 기울기 : 1.10
N >= 65536 만 : 1.19 <- 이론값 1 을 조금 넘는다 (log 때문)
개의 표본을 0.16초에 변환한다. 같은 일을 행렬 곱으로 하려면 곱셈이 번 필요하고, 애초에 그 행렬을 메모리에 올릴 수가 없다.
print("메모리도 문제다. N x N 복소 행렬의 크기 :")
for N in (1024, 4096, 2**16, 2**20):
바이트 = N * N * 16
print(f" N={N:>9,} : {바이트/2**30:>14,.1f} GB")
print()
print("N = 2^20 이면 16 테라바이트가 넘는다. 행렬을 아예 만들 수가 없다.")
print("FFT 는 행렬을 만들지 않는다. 이것도 못지않게 중요한 장점이다.")메모리도 문제다. N x N 복소 행렬의 크기 :
N= 1,024 : 0.0 GB
N= 4,096 : 0.2 GB
N= 65,536 : 64.0 GB
N=1,048,576 : 16,384.0 GB
N = 2^20 이면 16 테라바이트가 넘는다. 행렬을 아예 만들 수가 없다.
FFT 는 행렬을 만들지 않는다. 이것도 못지않게 중요한 장점이다.
마치며...¶
| 서술 파트의 내용 | 이 노트북의 코드 |
|---|---|
| 가 0 | x @ x 가 0, np.vdot(x,x) 가 2 |
| 파이썬 함정 | A.T @ A 는 대각조차 실수가 아니다 |
| 에르미트 | 무작위 300개의 고윳값 허수부가 10-15 |
| 유니타리 | , 길이 보존 |
| 의 두 성질 | 여러 에서 확인. 여덟 점의 합이 0 |
| 손으로 쓴 것과 일치 | |
| 열이 순수한 음 | 슬라이더로 를 올리면 진동이 빨라진다 |
| 까지 확인 | |
| 역변환이 공짜 | 이 inv(F) 와 일치 |
| 신호 분석 | 섞인 5와 17을 도로 갈라낸다 |
| 잡음 제거 | 큰 계수만 남기고 역변환 |
| 인수분해 | 세 조각의 곱이 |
| FFT | 재귀 열 줄이 np.fft.fft 와 일치 |
| 속도 | 로그-로그 기울기가 2 대 1 |
더 해 볼 것¶
6절의 신호에서 잡음 크기를 0.6에서 3.0으로 올려 보자. 봉우리가 여전히 보이는가? 몇까지 올리면 묻히는가?
6절의 자르는 기준을 100에서 바꿔 가며 복원 오차를 그려 보자. 너무 많이 자르면 무슨 일이 생기는가?
나의fft에 이 2의 거듭제곱이 아닌 입력을 넣으면 어떻게 되는가? 왜 그런가?np.fft.fft는 왜 괜찮은가?의 고윳값을 구해 보자. 절댓값이 전부 인가? 고윳값 자체는 몇 종류인가? ( 라는 사실과 관계가 있다.)
다음 강의에서는 실수의 세계로 돌아간다. 대칭행렬의 고윳값이 전부 양수라면 무슨 일이 벌어지는지 본다.