サンプリング定理とエイリアシング：ナイキスト周波数とアンチエイリアシングフィルタ

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, sosfilt, resample, freqz

# --- パラメータ設定 ---
f_signal = 130       # 信号周波数 [Hz]
fs_high   = 1000     # 十分高いサンプリング周波数（エイリアシングなし）
fs_low    = 200      # 不十分なサンプリング周波数（エイリアシング発生）
duration  = 0.1      # 表示時間 [s]

# --- 元の信号（高サンプリングレートで生成）---
t_high = np.arange(0, duration, 1 / fs_high)
x_orig = np.sin(2 * np.pi * f_signal * t_high)

# --- 低サンプリングレートでサンプリング ---
t_low = np.arange(0, duration, 1 / fs_low)
x_alias = np.sin(2 * np.pi * f_signal * t_low)

# エイリアス周波数を計算式で確認
f_alias = abs(f_signal - round(f_signal / fs_low) * fs_low)
print(f"信号周波数: {f_signal} Hz")
print(f"サンプリング周波数: {fs_low} Hz (ナイキスト: {fs_low//2} Hz)")
print(f"エイリアス周波数: {f_alias} Hz")
# => エイリアス周波数: 70 Hz

# --- プロット（時間領域） ---
fig, axes = plt.subplots(2, 1, figsize=(10, 6))

axes[0].plot(t_high * 1000, x_orig, 'b-', label=f'元信号 ({f_signal} Hz)', lw=1.5)
axes[0].plot(t_low * 1000, x_alias, 'ro', ms=8, label=f'低fs={fs_low}Hzでのサンプル')
# エイリアス波形を重ねて表示
t_fine = np.linspace(0, duration, 2000)
axes[0].plot(t_fine * 1000,
             np.sin(2 * np.pi * f_alias * t_fine),
             'r--', lw=1.5, label=f'エイリアス波形 ({f_alias} Hz)')
axes[0].set_xlabel('Time [ms]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('エイリアシング：時間領域')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# --- FFTによる周波数領域の比較 ---
N_high = len(x_orig)
N_low  = len(x_alias)

X_high = np.fft.rfft(x_orig)
X_low  = np.fft.rfft(x_alias)

freqs_high = np.fft.rfftfreq(N_high, 1 / fs_high)
freqs_low  = np.fft.rfftfreq(N_low,  1 / fs_low)

axes[1].plot(freqs_high, 2 / N_high * np.abs(X_high),
             'b-', label=f'fs={fs_high} Hz（エイリアシングなし）')
axes[1].plot(freqs_low, 2 / N_low * np.abs(X_low),
             'r-', label=f'fs={fs_low} Hz（エイリアシングあり）')
axes[1].axvline(fs_low / 2, color='k', ls='--', label=f'ナイキスト周波数 ({fs_low//2} Hz)')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('エイリアシング：周波数領域')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

fs = 1000   # サンプリング周波数 [Hz]
fn = fs / 2  # ナイキスト周波数

# 様々な入力周波数に対するエイリアス周波数を計算
f_input = np.linspace(0, 3 * fs, 3000)
f_alias = np.abs(f_input - np.round(f_input / fs) * fs)

plt.figure(figsize=(10, 5))
plt.plot(f_input, f_alias, 'b-', lw=2)
plt.axhline(fn, color='r', ls='--', label=f'ナイキスト周波数 = {fn:.0f} Hz')
plt.fill_between([0, fn], [0, 0], [fn, fn],
                 alpha=0.1, color='green', label='エイリアスなし領域')
plt.xlabel('入力周波数 [Hz]')
plt.ylabel('観測周波数（エイリアス周波数）[Hz]')
plt.title(f'周波数折り返し特性 (fs = {fs} Hz)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(0, fn * 1.1)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, sosfilt, freqz

np.random.seed(42)

# --- 帯域外成分を含む信号の生成 ---
fs_orig = 10000   # 元の高サンプリングレート [Hz]
fs_target = 1000  # ダウンサンプリング後のレート [Hz]
M = fs_orig // fs_target  # ダウンサンプリング比

duration = 0.5
t = np.arange(0, duration, 1 / fs_orig)

# 帯域内信号（100Hz）+ 帯域外信号（600Hz: ダウン後にエイリアシング）
x = np.sin(2 * np.pi * 100 * t) + 0.5 * np.sin(2 * np.pi * 600 * t)

# --- アンチエイリアシングフィルタの設計（バターワース LPF）---
# カットオフ: ナイキスト周波数の0.9倍
fc = (fs_target / 2) * 0.9  # = 450 Hz
sos = butter(8, fc, btype='low', fs=fs_orig, output='sos')
x_filtered = sosfilt(sos, x)

# --- ダウンサンプリング（AAFなし vs あり）---
x_down_noaaf = x[::M]          # AAFなし
x_down_aaf   = x_filtered[::M] # AAFあり

# --- 周波数応答の比較 ---
N_down = len(x_down_noaaf)
freqs_down = np.fft.rfftfreq(N_down, 1 / fs_target)

X_noaaf = np.fft.rfft(x_down_noaaf)
X_aaf   = np.fft.rfft(x_down_aaf)

fig, axes = plt.subplots(3, 1, figsize=(10, 10))

# フィルタの周波数応答
w, h = freqz(sos, worN=2048, fs=fs_orig, whole=False)
axes[0].plot(w, 20 * np.log10(np.abs(h) + 1e-12), 'b-', lw=2)
axes[0].axvline(fc, color='r', ls='--', label=f'カットオフ fc = {fc:.0f} Hz')
axes[0].axvline(fs_target / 2, color='g', ls='--',
                label=f'ナイキスト周波数 = {fs_target//2} Hz')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('アンチエイリアシングフィルタの周波数応答（8次バターワース）')
axes[0].set_xlim(0, 2000)
axes[0].set_ylim(-80, 5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# AAFなしのスペクトル
axes[1].plot(freqs_down, 20 * np.log10(2 / N_down * np.abs(X_noaaf) + 1e-12),
             'r-', lw=1.5)
axes[1].axvline(fs_target / 2, color='k', ls='--', alpha=0.5)
axes[1].set_title('ダウンサンプリング後のスペクトル（AAFなし）— 600Hzが400Hzにエイリアス')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Magnitude [dB]')
axes[1].set_xlim(0, fs_target / 2)
axes[1].grid(True, alpha=0.3)

# AAFありのスペクトル
axes[2].plot(freqs_down, 20 * np.log10(2 / N_down * np.abs(X_aaf) + 1e-12),
             'b-', lw=1.5)
axes[2].axvline(fs_target / 2, color='k', ls='--', alpha=0.5)
axes[2].set_title('ダウンサンプリング後のスペクトル（AAFあり）— エイリアスが除去')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Magnitude [dB]')
axes[2].set_xlim(0, fs_target / 2)
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import resample, butter, sosfilt

# --- 元信号の生成 ---
fs_orig = 1000   # 元のサンプリングレート
duration = 1.0
t_orig = np.arange(0, duration, 1 / fs_orig)
x_orig = (np.sin(2 * np.pi * 50 * t_orig)
          + 0.3 * np.sin(2 * np.pi * 120 * t_orig))

# --- ダウンサンプリング（1000 Hz → 400 Hz、デシメーション比 M=2.5）---
fs_down = 400
n_down = int(len(x_orig) * fs_down / fs_orig)
x_down = resample(x_orig, n_down)  # FFTベース：内部でAAFを適用
t_down = np.linspace(0, duration, n_down)

# --- アップサンプリング（元の1000 Hz → 4000 Hz）---
fs_up = 4000
n_up = int(len(x_orig) * fs_up / fs_orig)
x_up = resample(x_orig, n_up)      # 補間によるアップサンプリング
t_up = np.linspace(0, duration, n_up)

# --- スペクトル比較 ---
def spectrum(x, fs):
    N = len(x)
    X = np.fft.rfft(x)
    f = np.fft.rfftfreq(N, 1 / fs)
    return f, 2 / N * np.abs(X)

fig, axes = plt.subplots(3, 1, figsize=(10, 9))

f_orig, A_orig = spectrum(x_orig, fs_orig)
f_down, A_down = spectrum(x_down, fs_down)
f_up,   A_up   = spectrum(x_up,   fs_up)

axes[0].plot(f_orig, A_orig, 'b-')
axes[0].set_title(f'元信号のスペクトル (fs = {fs_orig} Hz)')
axes[0].set_ylabel('Amplitude')
axes[0].set_xlim(0, fs_orig / 2)
axes[0].grid(True, alpha=0.3)

axes[1].plot(f_down, A_down, 'r-')
axes[1].set_title(f'ダウンサンプリング後 (fs = {fs_down} Hz)')
axes[1].set_ylabel('Amplitude')
axes[1].set_xlim(0, fs_down / 2)
axes[1].grid(True, alpha=0.3)

axes[2].plot(f_up, A_up, 'g-')
axes[2].set_title(f'アップサンプリング後 (fs = {fs_up} Hz)')
axes[2].set_ylabel('Amplitude')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_xlim(0, fs_up / 2)
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"元信号:             {len(x_orig)} サンプル ({fs_orig} Hz)")
print(f"ダウンサンプリング: {len(x_down)} サンプル ({fs_down} Hz)")
print(f"アップサンプリング: {len(x_up)} サンプル ({fs_up} Hz)")

from scipy.signal import resample_poly

# 1000 Hz → 441 Hz（L=441, M=1000 のポリフェーズ）
x_poly = resample_poly(x_orig, 441, 1000)
print(f"ポリフェーズリサンプリング: {len(x_orig)} → {len(x_poly)} サンプル")