フーリエ変換

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
N = 1024  # サンプル数
T = 1.0 / 800.0  # サンプリング間隔
x = np.linspace(0.0, N*T, N, endpoint=False)

# 元の信号を生成（例：2つの異なる周波数の正弦波を合成）
freq1 = 50.0  # 1つ目の正弦波の周波数
freq2 = 120.0  # 2つ目の正弦波の周波数
signal = np.sin(freq1 * 2.0 * np.pi * x) + 0.5 * np.sin(freq2 * 2.0 * np.pi * x)

# フーリエ変換
fft_signal = np.fft.fft(signal)

# 逆フーリエ変換で信号を復元
reconstructed_signal = np.fft.ifft(fft_signal)

# プロット
plt.figure(figsize=(12, 6))

# 元の信号をプロット
plt.subplot(2, 1, 1)
plt.plot(x, signal, label='Original Signal')
plt.title('Original Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()

# 復元された信号をプロット
plt.subplot(2, 1, 2)
plt.plot(x, reconstructed_signal.real, label='Reconstructed Signal', linestyle='--')
plt.title('Reconstructed Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# サンプル数
N = 1024
# サンプリング周波数
Fs = 1000.0
# 時間軸
t = np.arange(N) / Fs

# 周波数成分を持つ信号を生成
freq1 = 50.0  # 50 Hz
freq2 = 120.0  # 120 Hz
signal = 0.7 * np.sin(2 * np.pi * freq1 * t) + np.sin(2 * np.pi * freq2 * t)

# 離散フーリエ変換 (DFT) の計算
dft = np.fft.fft(signal)

# 高速フーリエ変換 (FFT) の計算
fft = np.fft.fft(signal)

# 周波数軸の計算
freqs = np.fft.fftfreq(N, 1/Fs)

# プロット
plt.figure(figsize=(12, 6))

# 元の信号をプロット
plt.subplot(3, 1, 1)
plt.plot(t, signal)
plt.title('Original Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')

# DFTの結果をプロット
plt.subplot(3, 1, 2)
plt.plot(freqs[:N//2], np.abs(dft)[:N//2] * 1/N)  # 1/Nで正規化
plt.title('Discrete Fourier Transform (DFT)')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')

# FFTの結果をプロット
plt.subplot(3, 1, 3)
plt.plot(freqs[:N//2], np.abs(fft)[:N//2] * 1/N)  # 1/Nで正規化
plt.title('Fast Fourier Transform (FFT)')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def initial_temperature_distribution(x):
    """ 初期温度分布を定義する関数 """
    return np.sin(np.pi * x)  # 例として、sin(πx)の形の初期分布を使用

def heat_conduction_solution(x, t, alpha, n_terms):
    """ フーリエ級数を用いて熱伝導方程式の解を計算する関数 """
    solution = np.zeros_like(x)
    for n in range(1, n_terms + 1):
        coefficient = (2 / n) * (1 - (-1) ** n)
        solution += coefficient * np.exp(-alpha * (n * np.pi) ** 2 * t) * np.sin(n * np.pi * x)
    return solution

# パラメータ設定
L = 1.0  # 棒の長さ
alpha = 0.01  # 熱拡散率
n_terms = 50  # フーリエ級数の項数

# 空間と時間の離散化
x = np.linspace(0, L, 100)
times = [0, 0.01, 0.1, 0.5, 1.0]  # 時刻のリスト

# 初期温度分布をプロット
plt.figure(figsize=(12, 8))
for t in times:
    u_xt = heat_conduction_solution(x, t, alpha, n_terms)
    plt.plot(x, u_xt, label=f't={t:.2f}')

plt.title('Heat Conduction Solution using Fourier Series')
plt.xlabel('Position x')
plt.ylabel('Temperature u(x,t)')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def rectangle_wave(x, period=2*np.pi):
    """ 矩形波を定義する関数 """
    return np.where(np.mod(x, period) < period / 2, 1, -1)

def fourier_series(x, n_terms, period=2*np.pi):
    """ 矩形波のフーリエ級数展開を計算する関数 """
    result = np.zeros_like(x)
    for n in range(1, n_terms + 1):
        result += (4 / (np.pi * (2*n - 1))) * np.sin((2*n - 1) * x)
    return result

# パラメータ設定
x = np.linspace(-2*np.pi, 2*np.pi, 1000)  # xの範囲
n_terms_list = [1, 3, 5, 10, 20]  # 使用するフーリエ級数の項数

# 矩形波をプロット
plt.figure(figsize=(12, 10))
for n_terms in n_terms_list:
    y_approx = fourier_series(x, n_terms)
    plt.plot(x, y_approx, label=f'{n_terms} terms')

# 元の矩形波をプロット
plt.plot(x, rectangle_wave(x), linestyle='--', color='black', label='Rectangle Wave', linewidth=2)

plt.title('Fourier Series Convergence')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid(True)
plt.ylim(-1.5, 1.5)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータ設定
L = 10.0  # 領域のサイズ
N = 64  # グリッドの分割数

# グリッドの作成
x = np.linspace(0, L, N, endpoint=False)
y = np.linspace(0, L, N, endpoint=False)
z = np.linspace(0, L, N, endpoint=False)
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')

# 3D関数の生成（例: サイン波関数）
F = np.sin(2 * np.pi * X) * np.sin(2 * np.pi * Y) * np.sin(2 * np.pi * Z)

# フーリエ変換
F_hat = np.fft.fftn(F)

# 周波数軸の作成
kx = np.fft.fftfreq(N, d=(L/N))
ky = np.fft.fftfreq(N, d=(L/N))
kz = np.fft.fftfreq(N, d=(L/N))
KX, KY, KZ = np.meshgrid(kx, ky, kz, indexing='ij')

# フーリエ変換結果の振幅
F_hat_mag = np.fft.fftshift(np.abs(F_hat))  # 中心化

# プロット設定
fig = plt.figure(figsize=(14, 6))

# 元の関数のプロット
ax1 = fig.add_subplot(121, projection='3d')
slice_index = N // 2
ax1.plot_surface(X[:, :, slice_index], Y[:, :, slice_index], F[:, :, slice_index], cmap='viridis')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('F(x, y, z)')
ax1.set_title('Original Function (Slice in z)')

# フーリエ変換結果のプロット
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_surface(KX[:, :, slice_index], KY[:, :, slice_index], F_hat_mag[:, :, slice_index], cmap='viridis')
ax2.set_xlabel('kx')
ax2.set_ylabel('ky')
ax2.set_zlabel('Magnitude')
ax2.set_title('Magnitude of Fourier Transform (Slice in kz)')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.stats import norm

# 定数の設定
fs = 48000  # サンプリング周波数
f = 1000  # 正弦波の周波数
A = 10**(-90/20)  # 正弦波の振幅（-90dB）
N = fs  # サンプル数

# 正弦波の生成
t = np.arange(N) / fs
sine_wave = A * np.sin(2 * np.pi * f * t)

# 端数処理関数
def round_zero(x):
    return np.floor(x)

def round_infinity(x):
    return np.ceil(x)

def round_nearest(x):
    return np.round(x)

def round_nearest_even(x):
    return np.round(x/2)*2

# 端数処理の適用
sine_wave_rz = round_zero(sine_wave)
sine_wave_ri = round_infinity(sine_wave)
sine_wave_rn = round_nearest(sine_wave)
sine_wave_rne = round_nearest_even(sine_wave)

# ディザリング関数
def dither_uniform(x):
    noise = np.random.uniform(0, 1, len(x))
    return np.round(x + noise) - noise

def dither_gaussian(x):
    noise = norm.rvs(scale=1, size=len(x))
    return np.round(x + noise) - noise

def dither_pink(x):
    noise = np.random.randn(len(x))
    noise = np.convolve(noise, [0.05, 0.25, 0.5, 0.25, 0.05], mode='same')
    return np.round(x + noise) - noise

# ディザリングの適用
sine_wave_dither_uniform = dither_uniform(sine_wave)
sine_wave_dither_gaussian = dither_gaussian(sine_wave)
sine_wave_dither_pink = dither_pink(sine_wave)

# FFTの計算
def compute_fft(x):
    return np.abs(fft(x))[:len(x)//2]

# 結果のプロット
def plot_results(signals, titles):
    plt.figure(figsize=(12, 10))
    for i, (signal, title) in enumerate(zip(signals, titles)):
        plt.subplot(4, 2, i+1)
        fft_result = compute_fft(signal)
        fft_result[fft_result == 0] = 1e-12  # 非常に小さい値に置き換え
        plt.plot(np.linspace(0, fs/2, len(fft_result)), 20 * np.log10(fft_result))
        plt.title(title)
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Magnitude (dB)')
        plt.grid(True)
    plt.tight_layout()
    plt.show()

signals = [sine_wave_rz, sine_wave_ri, sine_wave_rn, sine_wave_rne, sine_wave_dither_uniform, sine_wave_dither_gaussian, sine_wave_dither_pink]
titles = ['RZ (Zero Direction)', 'RI (Infinity Direction)', 'RN (Round Nearest)', 'RNE (Round Nearest Even)', 'Dither (Uniform)', 'Dither (Gaussian)', 'Dither (Pink Noise)']

plot_results(signals, titles)

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import stft

# サンプル信号の作成
fs = 1000  # サンプリング周波数
t = np.linspace(0, 2, 2 * fs, endpoint=False)
x = np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 120 * t) + np.sin(2 * np.pi * 300 * t)

# STFTの計算
f, t, Zxx = stft(x, fs=fs, nperseg=256)

# STFTのプロット
plt.figure(figsize=(10, 6))
plt.pcolormesh(t, f, np.abs(Zxx), shading='gouraud')
plt.title('STFT Magnitude')
plt.xlabel('Time [sec]')
plt.ylabel('Frequency [Hz]')
plt.colorbar(label='Magnitude')
plt.show()