takecの気まぐれブログ

プログラミング等の話題を気まぐれに

窓関数の種類について

DFTを行う場合、入力信号がDFT領域に整数周期で入っていないとリーケージ誤差(漏れ誤差)が発生する。 この誤差の影響を抑えるために窓関数を掛けるが、種類が色々ありどのような特性を持っているのか把握しきれていない。 とりあえず、scipy.signalに含まれる窓関数とよく使っている7次ブラックマンハリスについて特性の一覧を出してみた。

f:id:takecccc:20170703053301p:plain

import numpy as np
from scipy import signal
from scipy.fftpack import fft, fftfreq, fftshift
import matplotlib as mpl
import matplotlib.pyplot as plt
#import subprocess
#print("subprocess imported")

N = 64
FFTN = 4096
filename = "window_function.png"
n = np.linspace(0,N,N,False)

w = 2*np.pi*n/N

# generate window
w_boxcar = signal.boxcar(N, sym=False)
w_barthann = signal.barthann(N, sym=False)
w_bartlett = signal.bartlett(N, sym=False)
w_blackman = signal.blackman(N, sym=False)
w_4_blackmanharris = signal.blackmanharris(N, sym=False)
w_7_blackmanharris = (
        0.27105140069342
        -0.43329793923448*np.cos(w)
        +0.21812299954311*np.cos(2*w)
        -0.06592544638803*np.cos(3*w)
        +0.01081174209837*np.cos(4*w)
        -0.00077658482522*np.cos(5*w)
        +0.00001388721735*np.cos(6*w)
        )
w_bohman = signal.bohman(N, sym=False)
w_chebwin = signal.chebwin(N, at=100, sym=False)
w_cosine = signal.cosine(N, sym=False)
w_exponential = signal.exponential(N, tau=30, sym=False)
w_flattop = signal.flattop(N, sym=False)
w_gaussian = signal.gaussian(N, std=N/8, sym=False)
w_general_gaussian = signal.general_gaussian(N, p=1.5, sig=N/6, sym=False)
w_hamming = signal.hamming(N, sym=False)
w_hann = signal.hann(N, sym=False)
w_kaiser = signal.kaiser(N, beta=np.pi*4, sym=False)
w_nuttall = signal.nuttall(N, sym=False)
w_parzen = signal.parzen(N, sym=False)
w_slepian = signal.slepian(N, width=0.1, sym=False)
w_triang = signal.triang(N, sym=False)
w_tukey = signal.tukey(N, sym=False)


#描画設定
#mpl.rcParams['font.family'] = 'sans-serif'
mpl.rcParams['font.family']         = 'Times New Roman'
mpl.rcParams['mathtext.default']    = 'regular'
mpl.rcParams['xtick.top']           = 'True'
mpl.rcParams['ytick.right']         = 'True'
mpl.rcParams['xtick.direction']     = 'in'
mpl.rcParams['ytick.direction']     = 'in'
mpl.rcParams['axes.grid']           = 'True'
mpl.rcParams['axes.xmargin']        = '0'
mpl.rcParams['axes.ymargin']        = '.05'
mpl.rcParams['savefig.facecolor']   = '#ffffff'
mpl.rcParams['savefig.edgecolor']   = 'None'
mpl.rcParams['savefig.bbox']        = 'tight'
mpl.rcParams['axes.titlesize']      = 12
mpl.rcParams['axes.labelsize']      = 12
mpl.rcParams['lines.linewidth']     = 1.0
mpl.rcParams['lines.markersize']    = 10
mpl.rcParams['xtick.labelsize']     = 10
mpl.rcParams['ytick.labelsize']     = 10

fig, axs = plt.subplots(8,3,figsize=(14,28))
for ax in axs[:,0]:
    ax.set_xlim((0,2*np.pi))
    ax.set_xticks(np.linspace(0, 2*np.pi, 11))
    ax.set_xticklabels(["{0:.1f}$\pi$".format(x/np.pi) for x in ax.get_xticks()])
    ax.set_xlabel("Phase [rad]")
    ax.set_ylim((-0.2,1.2))
    ax.set_ylabel("Amplitude")
    ax.minorticks_on()
    ax.tick_params(which='both', direction='in', zorder=1.8)
    ax.grid(True,which='major')

for ax in axs[:,1]:
    ax.set_xlim((-10,10))
    ax.set_xlabel("bins")
    ax.set_ylim((-200,10))
    ax.set_ylabel("Magnitude [dB]")
    ax.grid(True,which='both')

for ax in axs[:,2]:
    ax.set_xlim((-N/2,N/2))
    ax.set_xlabel("bins")
    ax.set_ylim((-200,10))
    ax.set_ylabel("Magnitude [dB]")
    ax.grid(True,which='major')

#freq = np.linspace(-N/2, N/2, FFTN,False)
freq = fftshift(fftfreq(FFTN,1/N))
def mag(w):
    w_normalize = w/(np.sum(w)/len(w))
    A = fft(w_normalize, FFTN) / (len(w))
    responce = 20 * np.log10(np.abs(fftshift(A)+1E-32))
    #responce[responce<-200] = None
    return responce

def mags(ws):
    return np.array([mag(w) for w in ws])

g0 = np.array([
    w_boxcar,
    w_hann,
    w_hamming,
    ])
g0_legend = [
        "boxcar",
        "hann",
        "hamming",
        ]
g1 = np.array([
    w_bohman,
    w_chebwin,
    ])
g1_legend = [
        "bohman",
        "chebwin",
        ]
g2 = np.array([
    w_blackman,
    w_nuttall,
    w_4_blackmanharris,
    w_kaiser,
    w_7_blackmanharris,
    ])
g2_legend = [
        "blackman",
        "nuttall",
        "4 term blackman-harris",
        "kaiser",
        "7 term blackman-harris",
        ]
g3 = np.array([
    w_flattop,
    ])
g3_legend = [
        "flattop",
        ]
g4 = np.array([
    w_tukey,
    w_parzen,
    ])
g4_legend = [
        "tukey",
        "parzen",
        ]
g5 = np.array([
    w_triang,
    w_bartlett,
    w_barthann,
    w_hann,
    ])
g5_legend = [
        "triang",
        "bartlett",
        "barthann",
        "hann",
        ]
g6 = np.array([
    w_gaussian,
    w_general_gaussian,
    ])
g6_legend = [
        "gaussian",
        "general_gaussian",
        ]
g7 = np.array([
    w_cosine,
    w_slepian,
    ])
g7_legend = [
        "cosine",
        "slepian",
        ]
axs[0,1].set_ylim((-150,10))
axs[0,2].set_ylim((-150,10))
axs[0,0].plot(w, g0.T)
axs[0,1].plot(freq, mags(g0).T)
axs[0,2].plot(freq, mags(g0).T)
axs[0,2].legend(
        g0_legend,
        bbox_to_anchor = (1.05, 1),
        loc = 'upper left',
        borderaxespad = 0
        )
axs[1,1].set_ylim((-150,10))
axs[1,2].set_ylim((-150,10))
axs[1,0].plot(w, g1.T)
axs[1,1].plot(freq, mags(g1).T)
axs[1,2].plot(freq, mags(g1).T)
axs[1,2].legend(
        g1_legend,
        bbox_to_anchor = (1.05, 1),
        loc = 'upper left', borderaxespad = 0
        )
axs[2,1].set_ylim((-200,10))
axs[2,2].set_ylim((-200,10))
axs[2,0].plot(w, g2.T)
axs[2,1].plot(freq, mags(g2).T)
axs[2,2].plot(freq, mags(g2).T)
axs[2,2].legend(
        g2_legend,
        bbox_to_anchor = (1.05, 1),
        loc = 'upper left', borderaxespad = 0
        )
axs[3,1].set_ylim((-0.5,0.5))
axs[3,2].set_ylim((-100,10))
axs[3,1].set_xlim((-2,2))
axs[3,0].plot(w, g3.T)
axs[3,1].plot(freq, mags(g3).T)
axs[3,2].plot(freq, mags(g3).T)
axs[3,2].legend(
        g3_legend,
        bbox_to_anchor = (1.05, 1),
        loc = 'upper left', borderaxespad = 0
        )
axs[4,1].set_ylim((-150,10))
axs[4,2].set_ylim((-150,10))
axs[4,0].plot(w, g4.T)
axs[4,1].plot(freq, mags(g4).T)
axs[4,2].plot(freq, mags(g4).T)
axs[4,2].legend(
        g4_legend,
        bbox_to_anchor = (1.05, 1),
        loc = 'upper left', borderaxespad = 0
        )
axs[5,1].set_ylim((-120,10))
axs[5,2].set_ylim((-120,10))
axs[5,0].plot(w, g5.T)
axs[5,1].plot(freq, mags(g5).T)
axs[5,2].plot(freq, mags(g5).T)
axs[5,2].legend(
        g5_legend,
        bbox_to_anchor = (1.05, 1),
        loc = 'upper left', borderaxespad = 0
        )
axs[6,1].set_ylim((-100,10))
axs[6,2].set_ylim((-100,10))
axs[6,0].plot(w, g6.T)
axs[6,1].plot(freq, mags(g6).T)
axs[6,2].plot(freq, mags(g6).T)
axs[6,2].legend(
        g6_legend,
        bbox_to_anchor = (1.05, 1),
        loc = 'upper left', borderaxespad = 0
        )
axs[7,1].set_ylim((-80,10))
axs[7,2].set_ylim((-80,10))
axs[7,0].plot(w, g7.T)
axs[7,1].plot(freq, mags(g7).T)
axs[7,2].plot(freq, mags(g7).T)
axs[7,2].legend(
        g7_legend,
        bbox_to_anchor = (1.05, 1),
        loc = 'upper left', borderaxespad = 0
        )

fig.savefig(filename)