ナイキスト周波数を超えたらどうなるかアニメーションで見る

はじめに

　サンプリング定理(標本化定理)というものがあります。1波長に2回サンプリングすれば、元の信号を完全に復元できるというものです。このサンプリングする間隔をサンプリングレートと呼び、その周波数の半分の値をナイキスト周波数と呼びます。計測対象の信号がナイキスト周波数を超えるとエイリアシングが起こります。
　このようなことを知識として知っていても、ナイキスト周波数を超えるとどのように見えるのかアニメーションで見たことがないなと思ったので作りました。

2020/06/06 追記
実空間のアニメーションも追加しました。

参考

サンプリング定理を思い返すのに読みました。
小野測器:FFTアナライザについて(page4)

プログラム

周波数空間

import numpy as np
import matplotlib.pyplot as plt

N = 512 #FFTのサンプル点数
fs = 2000.0 #サンプリングレート
t = np.linspace(0.0,N/fs,N)
fft_freq = np.linspace(0.0,fs,N) #FFTの周波数軸
f_array = np.linspace(10.0,4000.0,200)

for i,f in enumerate(f_array):
    sig = np.sin(2.0*np.pi*f*t) * np.hamming(N)
    fft_amp = np.abs(np.fft.fft(sig)) / N * 2.0 / 0.54
    
    fig,axes = plt.subplots()
    axes.plot(fft_freq,fft_amp,label="{} Hz".format(np.round(f)))
    axes.set_xlabel("Frequency Hz")
    axes.set_ylabel("Amplitude")
    axes.set_ylim(0.0,1.1)
    axes.legend(loc="upper right")
    fig.savefig("{}.png".format(i))
    plt.close()

実空間

def aliasing_freq(f,fs):
    """
    折り返された周波数
    """
    sign = (-1) ** int(f/(fs/2.0)) #周波数の符号
    if int(f/(fs/2.0))%2 == 0:
        return sign * (f - fs/2.0 * int(f/(fs/2.0+1.0e-6)))
    else:
        return sign * (fs/2.0-(f - fs/2.0 * int(f/(fs/2.0+1.0e-6))))

for i,f in enumerate(f_array):
    sig = np.sin(2.0*np.pi*f*t)
    sig_true = np.sin(2.0*np.pi*f*t_true)
    
    fig,axes = plt.subplots()
    axes.plot(t_true,sig_true,label="{} Hz".format(np.round(f)),alpha=0.5)
    axes.scatter(t,sig)
    axes.plot(t_true,np.sin(2.0*np.pi*aliasing_freq(f,fs)*t_true))
    axes.set_xlim(0.0,0.0+1.0/f*10)
    axes.set_ylim(-1.2,1.2)
    axes.set_xlabel("Time s")
    axes.set_ylabel("Amplitude")
    axes.legend(loc="upper right")
    fig.savefig("graph/{}.png".format(i))
    plt.close()

連番写真はGIMPでアニメーションに変換しました。

結果

周波数空間ではナイキスト周波数を境目に鏡のように動きます。

実空間でも同様です。

こうしてみるとなぜナイキスト周波数を超えてはいけないかがよくわかります。今回は周波数がわかりますが、結果だけ見れば折り返されたものなのか、正しく測れたものなのか見分けが付きません。

まとめ

ナイキスト周波数を超えてはいけないと覚えてはいても、超えたらどうなるのかを見る機会は少ないと思います。誰かの役に立てば幸いです。