久しぶりにmacのvs codeでpythonを実行したらエラーが発生。

ごくごく簡単なスクリプトで以前にも実行していたのに、なぜかsegmentation fault 11が発生している。

デバッグしたところ、plot.show()で発生している。

何かしたかな?と考えたところ、そういえば、OSをBIg Surにあげたと。

検索してみると、同様の問題は発生していて、matplotlibをアンインストールして再度インストールしろと。

で、ターミナルから
conda uninstall matplotlib
でアンインストールして
conda install matplotlib
でインストールした。

これで動くかと思ったら今度は
no module named scipy
とでてしまう。

conda install scipy
でインストールしたら、とりあえずスクリプトは実行できるようになった。

アップサンプリングする際に、サンプリング周波数を上げるためにゼロ点を挿入する。

ゼロ点を挿入した信号のスペクトラムをみてみると、もとの信号のスペクトラムが繰り返えされている。

ゼロ点を2点追加
resample_interpolate
スペクトラムが繰り返される
resample_fft_after_interpolate

このことを考えてみる。

サンプリング周波数をTとし、周波数をk/Tとする信号S \[ \cos(2 \pi \frac{k}{T} n ) \] を考える。

信号Sのフーリエ変換は \[ F(k') = \sum_{n=-T/2} ^{T/2}S[n] \cos(2 \pi \frac{k'}{T} n) \] と表現できる

信号Sのサンプル点間にa個のゼロ点を挿入してみる。ゼロ点を挿入した信号をSmとする。信号Smは \[ \cos(2 \pi \frac{k}{aT} n ) \ \ n=0,1,...,aT-1\] とあらわされる。

信号Smはn=am (m=0,1,...,T-1)のときのみ値を持つので \[ \cos(2 \pi \frac{k}{aT} am ) \ \ m=0,1,...,T-1 \] と表現できる。整理すると、 \[ \cos(2 \pi \frac{k}{T} m ) \ \ m=0,1,...,T-1 \] となる。

さて、信号Smのフーリエ変換を考えると \[ F(k') = \sum_{n=-aT/2}^{aT/2} Sm[n] \cos(2 \pi \frac{k'}{aT} n) \] となる。ただ、信号Sm[n]はn=am(m=0,1,...,T-1)のときのみ値を持つので \[ F(k') = \sum_{m=-T/2}^{T/2} Sm[m] \cos(2 \pi \frac{k'}{T} m) \] と表現することもできる。

信号Smは \[ \cos(2 \pi \frac{k}{T} m )\] であったので、フーリエ変換の結果として値を持つのは \[ \cos(2 \pi \frac{k'}{T} m) \] において \[ 2 \pi \frac{k'}{T} \pm 2\pi l = \pm 2 \pi \frac{k}{T} \] が成り立つ$k'$の値のときとなる。

したがって \[k' = \pm k \pm T l \] となり、この$k'$がナイキスト周波数以下の範囲$0 \le k' \le \frac{aT}{2}$を満たす範囲内で値を持つことになる。

例えば、a=3のとき$k'=k,\frac{aT}{2}\pm k$で値を持つことになる。

つまり、3回スペクトラムが繰り返されることになる。

リサンプリングの方法としてゼロ点を挿入・サンプル点を間引く方法と、フーリエ変換を使う方法の二種類が考えられる。

ゼロ点挿入・サンプル点間引き

N倍にアップサンプリングすることを考える。

まず、N-1のゼロを元信号のサンプル点の間に挿入して(インターポレーション)、サンプリング周波数をあげる。

以下の図は、3倍にアップサンプリングする場合の例。元信号のサンプル点の間に2個ゼロ点を挿入している。

元信号
resample_ori
サンプル点の間にゼロ点を2個入れた状態
resample_interpolate

上記をフーリエ変換して周波数領域で見ると、ゼロ点を挿入したほうは元信号の周波数表示がN回繰り返しされた状態になっている。

元信号
resample_fft_ori
サンプル点の間にゼロ点を2個入れた状態
resample_fft_after_interpolate

そこで、LPFで不要な繰り返しの部分をカットしてあげることで、元信号のサンプリング周波数をN倍した信号が得られる。

LPFで繰り返しの部分をカット
resample_LPF

元信号とアップサンプリングした信号を重ねて表示してみる。ほぼ重なっている。アップサンプリングすることで曲線がなめらかになっている分、元信号が少しカクカクにみえる。

resample_via_interploration

一方、N分の1にダウンサンプリングする場合は、N個ずつサンプルを抽出すればよい。ただし、サンプリング周波数が小さくなることで、ナイキスト周波数も小さくなるので、折り返しが発生しないように、サンプルを抽出する前に、LPFで帯域をナイキスト周波数におさまるようにしておく必要がある。ダウンサンプリング限定であればscipyにdecimateという関数がある。

ただ、上記の方法だと、整数倍の周波数変換しかできない。例えば、8kサンプリングを10kサンプリングにすることは単にゼロ点を挿入したり間引くことではできない。 ただ、ゼロ点挿入+間引きを合わせることで可能になる。例えば、例えば、8kサンプリングを10kサンプリングにする場合、まず、8kを40kをアップサンプリングして、その後40kを10kにダウンサンプリングすることで所望の周波数変換が実現できる。

フーリエ変換を用いた方法

フーリエ変換を用いた方法では、対象の信号をフーリエ変換して、周波数領域でサンプル点を追加または削除して、逆フーリエ変換をするこで、サンプル数の変換を実現している。

例えば、サンプル数Nの信号をMにアップサンプリングする場合、信号をフーリエ変換してデータ点Fをえる。Fの中央にM-Nのゼロ点を挿入しFのサンプル数をNにし、Fを逆フーリエ変換する。

サンプル点をNからMにダウンサンプリングするときは、フーリエ変換後のデータ点の中央からN-M個の点を削除して、逆フーリエ変換する。

フーリエ変換を用いる方法の利点は、任意のサンプル数を追加・削除できるので、任意の周波数変換が行える。

ただ、フーリエ変換を用いているので、元信号は周期的であることが前提。したがって、フーリエ変換を用いて周波数変換して得られた信号の終点は、始点につながるように変換される。

例えば、サンプリング周波数64の信号をサンプリング周波数150にアップサンプルすることを考える。

元信号をフーリエ変換した点列の中央に150-64=86個のゼロ点を挿入する。

元信号をフーリエ変換した状態
resample_fft_ori_150
中央にゼロ点を挿入しサンプル数を64から150に変更した状態
resample_fft_interpolate_150

逆フーリエ変換する。元信号とアップサンプルした信号を重ねて表示してみる。

resample_fft_result_150

信号が周期的であることを前提としているため、アップサンプルした信号の最初と最後は元信号とだいぶずれてしまっている。



"""
resample 
"""
# -*- coding: utf-8 -*-
import numpy as np
import scipy.signal as signal
import scipy
from scipy import io
from scipy.io import wavfile
import matplotlib.pyplot as plt
import sys, os

#read wav file
wav_path = "/Users/XXX/Downloads/samplewav.wav"
samplerate, wav_data = wavfile.read(wav_path)
fs = samplerate
data = wav_data[2*fs:3*fs]/(2**15)
#bc.showSpectrum(data, title='original')
#bc.showStem(data[:dip_sa], title='original')
"""
upsampleing via interpolation + LPF
"""
#
#interpolation
#
uprate = 3
data_up = np.zeros(len(data)*3, 'float')
for i, d in enumerate(data):
    data_up[i*uprate] = d
#bc.showSpectrum(data_up, title='after interpolation')
#bc.showStem(data_up[:dip_sa*uprate], title='after interpolation')
#
#LPF
#
tap_num = 151
cutoff=40 
b=scipy.signal.firwin(numtaps=tap_num, cutoff=fs/2, fs = 3*fs)
a = 1
data_lpf = filtered = signal.lfilter(b, a, data_up)
delay_sa = int((tap_num-1)/2)
data_lpf = data_lpf[delay_sa:]
#adjust Power
data_lpf *= uprate
#
#display result
#
start_sa = 500
disp_sa = 50
bc.showSpectrum(data_lpf, title='after LPF')
#bc.showPlot(data[:dip_sa])
#bc.showPlot(data_lpf[:3*dip_sa]*uprate)
data_ori = data[start_sa:start_sa+disp_sa]
t_list = range(0, disp_sa*uprate, uprate)
plt.plot(t_list, data_ori)
t_u_list = range(0, disp_sa*uprate)
plt.plot(t_u_list, data_lpf[start_sa*uprate:(start_sa+disp_sa)*uprate])
plt.legend(['original', 'upsample'])
plt.title('upsample via interpolation + LPF')
plt.show()


"""
upsampling via fft
"""
fftpt = 64
fftdt = np.fft.fft(data[:fftpt])
plt.plot(np.abs(fftdt))
plt.title('fft')
plt.show()
#upsampling 
sig_len = fftpt
upsamplerate = 150
insert_point_num = upsamplerate - sig_len
#insert zero around center
zero_point = np.zeros(insert_point_num, 'complex')
cnt_pt = int(fftpt/2)
fftdt_inserted = np.append(fftdt[:cnt_pt], zero_point)
fftdt_inserted = np.append(fftdt_inserted, fftdt[cnt_pt:])
plt.plot(np.abs(fftdt_inserted))
plt.title('fft insertesd')
plt.show()
#inverse fft
data_inv = np.fft.ifft(fftdt_inserted)
data_inv *= upsamplerate / sig_len
t_list = np.linspace(0, upsamplerate, sig_len, endpoint=False)
plt.plot(t_list, data[:fftpt])
t_u_list = range(0, upsamplerate)
plt.plot(t_u_list, data_inv.real)
plt.legend(['original', 'upsample'])
plt.title('upsample via FFT')
plt.show()

矩形をフーリエ変換した関数はsinc関数と名前がついている。

矩形区間の長さを2aとするとsinc関数は

\[ \frac{1}{\pi f}\sin2\pi a f \]

となる。

sinc_rect
sinc_sinc_tap_20

この式から、$f=\frac{n}{2a}$で零点となることがわかる。

scipyのsinc関数をみてみると

\[ \frac{1 }{\pi t}\sin\pi t \]

となってる。

したがって、sinc関数は$n=1,2,...$の整数値のところで零点となる。

sinc_discrete

離散化した場合をみてみる。

サンプリング周波数をNとして、矩形の長さを2aとする。

sinc_rect 2

この矩形のフーリエ変換は

\[ \frac{1}{\pi \frac{n}{N}}\sin2\pi a \frac{n}{N} \]

となる。

したがって、零点の位置$n$は

\[ 2\pi a \frac{n}{N} = \pi \]

から

\[ n = \frac{N}{2a} \times 定数\]

となる。

例えば、$N=1024$、$a=20$のとき、$n=\frac{1024}{20}\times\frac{1}{2} = 25.6 \times 定数$

となる。

sinc_dicrete

周波数領域での積は、時間領域では畳み込みに対応する。この関係性を利用すると、周波数を好みにカットする関数を作って、それを逆フーリエ変換するという方法で、フィルタを作ることができる。

逆フーリエ変換_フィルタ

このようなフィルタの設計法を窓関数設計という。

理想的なフィルタは無限に続きまた因果性も満たさない

ローパスフィルタを考えてみる。

理想的なローパスフィルタは、ある周波数以上の信号をバッサリ遮断するものだから、矩形をかけることにあたる。

矩形を逆フーリエ変換することでフィルタが得られるわけだが、矩形の逆フーリエ変換はsinc関数になる。

cut_sinc_sinc

sinc関数は無限に続いているし、時刻が負の時も値をもっているので因果性をみたさない。したがって理想的に周波数を遮断するフィルタは物理的に不可能といえる。

sinc関数を途中でカットする

妥協案として考えられるのが、無限に続くsinc関数をあるところでカットして有限にすること。あと、時刻が負の時に値も持つ件を解決するために、時間が正の方向にシフトする。そうすることで、実現可能なフィルタとなる。ただ、当然ながら周波数を遮断する性能は落ちる。では、上記のようなカットしたりシフトしたりする影響はどこにでてくるのか。

その影響は2つある。ひとつは遮断特性、もうひとつは振幅が一定でなくなる(リプル)。

cut_sinc_rec

理想的なフィルタでは通過する周波数と遮断される周波数が境界で完全にわかれていたが、途中でカットしたフィルタではなだらかに遮断されていくことになる。

また、通過する周波数において振幅が一定ではなく波を打つようになる。

遮断特性はフィルタの次数をあげることで急峻になっていくが、リプルは次数をあげて消える気配はない。

cut_sinc_cut_sinc

この点を考えてみる。

そもそもsinc関数をある範囲でバッサリカットするということは、矩形窓関数をかけることにあたる。したがって、周波数領域で矩形窓関数のフーリエ変換したものを理想的フィルタの形状に畳み込んでいることになる。

しだって、窓関数の次数をあげるとメインローブの幅が狭くなるので、通過と遮断の立ち上がりが急峻になる。

一方、リプルはサイドローブの影響。したがって、次数をあげてもサイドローブはなくならないので、リプルもなくならない。リプルを小さくするには、サイドローブが小さい窓関数でsinc関数をカットしてやればよいということになる。

例えば、ハミング窓を作用させることでリプルを小さくすることができる。ただし、遮断特性はなだらかになる。

cut_sinc_hamming
"""
Design filter by window
"""
# -*- coding: utf-8 -*-
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import sys, os
sys.path.append(os.pardir)  

fftpt = 1024
tap_num = 20
#sinc function
x = np.linspace(-tap_num,tap_num,fftpt)
s = np.sinc(x)
bc.showPlot(s, title='sinc func')
bc.showSpectrum(s)
#fft sinc in varing tap
tap_list = [3,5,10,15,20]
bin_per_tap = fftpt/(2*tap_num)
cut_sinc = np.zeros(fftpt, 'float')
f_list = np.linspace(-np.pi, np.pi, fftpt)
w_hamming = signal.hamming(fftpt)
for tap in tap_list:
    cut_sinc = s.copy()
    cent_bin = int(fftpt/2)
    cut_bin = int(tap * bin_per_tap)
    cut_sinc[:cent_bin - cut_bin] = 0
    cut_sinc[cent_bin + cut_bin:] = 0
    #bc.showPlot(cut_sinc, title='sinc func(cut)')
    cut_sinc = w_hamming * cut_sinc
    fftdt = np.fft.fft(cut_sinc, fftpt)/np.sqrt(fftpt)
    fftdt = np.fft.fftshift(fftdt)
    plt.plot(f_list, 20*np.log10(np.abs(fftdt)))
plt.legend(tap_list)
plt.title('fft sinc func (cut)')
plt.show()

窓関数によって周波数分解能とダイナミックレンジの性質は異なるのでその違いをみてみる。窓関数のフーリエ変換結果から窓関数の性質はみてとれるが、実際に信号に窓関数をかけたものをフーリエ変換してその性質を確認してみる。 ここでは以下の窓関数を比較してみる。

  • 矩形
  • hann
  • hamming
  • blackman
  • blackman-harris

比較に入る前に、窓関数の意義を考えるため、窓関数を作用させない場合(=矩形窓をかける場合)を見てみる。

フーリエ変換をするのに適当に信号を切り取るとその範囲がちょうど信号の周期になっているケースはほとんどないので、信号が不連続になり、その影響が周波数領域では周波数のもれとして現れる。

サンプリングレート:8192Hz
切り取り範囲(=fftポイント数):1024
周波数:513[Hz]

w_effect_rect_513

ちなみに、切り出した範囲が周波数の定数倍である場合は信号は不連続にはならず、したがって周波数もれも生じない。ピークがピンとたつ。

サンプリングレート:8192Hz
切り取り範囲(=fftポイント数):1024
周波数:512[Hz]

w_effect_rect_512

適当に信号を切り出すと不連続なものになってしまうが、 窓関数を作用させることで、信号を半ば無理やり連続的なものにし、周波数もれの影響を小さくすることが可能となる。

周波数分解能の違い

窓関数による周波数分解能の違いをみるために、周波数が隣接した2つの信号に窓関数かけフーリエ変換してみる。

サンプリングレート:8192Hz
切り取り範囲(=fftポイント数):1024
周波数:513[Hz]と523[Hz]

矩形
w_effect_res_rect
hann
w_effect_res_hann
hamming
w_effect_res_hamming
blackman
w_effect_res_blackman
blackman-harris
w_effect_res_blackman_harris

このフーリエ変換の結果の1binあたりの周波数は$8192/1024 = 8[\rm{Hz}]$である。

なにもしない(=矩形)場合がメインビームの幅が一番せまい。矩形窓のメインビームの幅の半分は 1binあたりの周波数 = 8[Hz]なので分離できている。

一方、hann窓やhamming窓は、メインビームの幅の半分が2x(1binあたりの周波数)=2x8[Hz] = 16[Hz]であるので、10[Hz]差の信号の分離ができていない。

ダイナミックレンジの違い

窓関数によるダイナミックレンジの違いをみるために、電力に差がある周波数の離れた2つの信号に窓関数かけフーリエ変換してみる。

サンプリングレート:8192Hz
切り取り範囲(=fftポイント数):1024
周波数:100[Hz]と900[Hz]

矩形
w_effect_dyn_rect
hann
w_effect_dyn_hann
hamming
w_effect_dyn_hamming
blackman
w_effect_dyn_blackman
blackman-harris w_effect_dyn_blackman_harris

矩形はメインビームは細いが、サイドローブが大きいので、900[Hz]の信号は確認できない。また、hamming窓もサイドローブが大きいので900[Hz]の信号は確認できない。

"""
compare window function effect
"""
# -*- coding: utf-8 -*-
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import sys, os


#parameter
freq = 100
freq_diff = 800
pw_diff = 1E-6
T=3000
fftpt=1024
N=fftpt
freq_per_bin = T/fftpt
#generate signal
sig = lambda f1, T: np.exp(2*np.pi*1j * f1 *np.arange(T)/T)
nos = lambda T : np.random.randn(T) + 1j * np.random.randn(T)
s = sig(freq, T) + pw_diff * sig(freq + freq_diff, T) + nos(T) * 0
#window function
w_rec = np.ones(N)
w_hann = signal.hann(N)
w_hamming = signal.hamming(N)
w_blackman = signal.blackman(N)
w_blackmanharris = signal.blackmanharris(N)
w_dic = {'rectangle':w_rec, 'hann':w_hann, 
    'hamming':w_hamming, 'blackman':w_blackman, 'blackman_harris':w_blackmanharris}
#display window function
f = lambda x : np.fft.fftshift(np.fft.fft(x, fftpt)/np.sqrt(fftpt))
#freq_list = 2 * np.pi * np.linspace(-0.5, 0.5, fftpt)
freq_list = T * np.linspace(-0.5, 0.5, fftpt)
w_list = []
for k in w_dic:
    w_list.append(k)
    s_w = s[:fftpt] * w_dic[k]
    fftdt = f(s_w)
    fftdt_log10 = 20*np.log10(np.abs(fftdt))
    plt.plot(freq_list, fftdt_log10)
    plt.legend([k])
    plt.ylabel("Amplitude[dB]")
    plt.xlabel("freq [Hz]")
    title = 'compare window function effect\n'
    title += 'sampling rate:%d[Hz] ' % T
    title += 'freq:%d[Hz] ' % freq
    title += 'fftpt:%d' % fftpt
    #title += 'resolution %.2f[Hz]' % freq_per_bin
    plt.title(title)
    #plt.xlim([40,160])
    plt.ylim([-150,40])
    plt.show()
plt.close('all')

カイザー窓はパラメータを変更することで柔軟に特性を変更できるので信号処理にもよく用いられる。

パラメータ$\beta(=\pi\alpha)$が0の場合は矩形窓。$\beta$を大きくすることで周波数分解能は悪くなるがダイナミックレンジが広くなる。つまり、$\beta$を変更するだけで周波数分解能とダイナミックレンジという両立しない2種類の窓関数の特性を連続的に推移できる。

また、カイザー窓はDPSSのよい近似にもなっている。DPSSとは、全電力に対するメインローブの電力比を最大化するという特徴をもっている。

kaiser_func
kaiser_fft
"""
compaire kaiser function for different beta
"""
# -*- coding: utf-8 -*-
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import sys, os


#beta
beta = {'ractangle 0':0, 'similar to a Hamming 5':5, 'similar to a Hann 6':6,
    'Similar to a Blackman 8.6':8.6}
w_list = []
N = 51
for k in beta:
    w_list.append(k)
    plt.plot(signal.kaiser(N, beta[k]))
plt.legend(w_list)
plt.ylabel("Amplitude")
plt.xlabel("sample")
plt.title('kaiser function')
plt.show()
plt.close('all')
#display frequency response of window function
fftpt=8192
f = lambda x : np.fft.fftshift(np.fft.fft(x, fftpt)/np.sqrt(fftpt))
freq_list = 2 * np.pi * np.linspace(-0.5, 0.5, fftpt)
#each
w_list = [] 
for k in beta:
    w_list.append(k)
    fftdt_lin = f(signal.kaiser(N, beta[k]))
    fftdt_log = 20*np.log10(np.abs(fftdt_lin))
    plt.plot(freq_list, fftdt_log)
plt.legend(w_list)
plt.ylabel("Amplitude[dB]")
plt.xlabel("bin")
plt.title('frequency response of window function')
plt.xlim([-1,1])
plt.ylim([-170,0])
plt.show()

窓関数の形状

以下の窓関数を考えてみる。
  • 矩形
  • han
  • hamming
  • blackman
  • blackman harris

各々の窓関数の形は以下のようになる。ただ、関数の形をみても性質はよくわからない。

win_all

窓関数の特徴をみるために、窓関数をフーリエ変換してみる。時間領域で対象関数に窓関数をかけるということは、周波数領域で畳み込みをするということであるため、窓関数をフーリエ変換した関数の形状をみることで、周波数の分解能やダイナミックレンジの性質を把握することができる。

窓関数の式

各々の窓関数の式は以下のようになる。

hann及びhammingは以下の式で表され

$a + b \cos(2\pi \frac{t}{N})$

hannは$a=b=1/2$、hamningは$a=0.54,b=0.46$としたもの。

blackmanは、

$a + b \cos(2\pi \frac{t}{N}) + \cos(4\pi \frac{t}{N})$

で$a=0.42,b=0.5,c=0.08$としたもの。

窓関数をフーリエ変換

上記の各々の窓関数をフーリエ変換してみる。

解析的に求めるための準備としてまず$\cos(2\pi n \frac{t}{N}) $のフーリエ変換を求めるておくとよい。

$\cos(2\pi n \frac{t}{N}) $のフーリエ変換は以下のようになる \[W_{n}(\omega) = (-1)^{n}\left(\frac{2 L^2 \omega}{L^2\omega^2 - 4 n^2 \pi^2}\right)\sin\left(\frac{N\omega}{2}\right) \]

展開式はこちらを参照。

例えば、hannやhammingのフーリエ変換は上記を利用して以下のように求められる。

\[a + b \cos(2\pi \frac{t}{N}) \xrightarrow{F} aW_{0} + bW_{1} \]

hannは$a=b=1/2$であるので

\[W(\omega) = \left(\frac{1}{\omega}-\frac{L^2 \omega}{L^2\omega^2 - 4 \pi^2}\right)\sin\left(\frac{N\omega}{2}\right) \]

となる。

グラフの形状は以下のようになり最初に0と交わる$\omega$は$\frac{4\pi}{N}$となる。$\sin$を見ると$\frac{2\pi}{N}$周期で0と交わるとはずだがそうななっていない。

$\frac{1}{\omega}-\frac{L^2 \omega}{L^2\omega^2 - 4 \pi^2}$の影響によりそうなっている。展開式はこちらの式(22)を参照。

各々の窓関数のフーリエ変換の形状をみてみる

矩形

フーリエ変換した式は$W_{0}=\frac{2}{\omega}\sin\left(\frac{N\omega}{2}\right)$となり、最初の節の位置は$\frac{2\pi}{N}$となる。

win_rec
hann

最初の節の位置は$\frac{4\pi}{N}$。矩形と比べると幅が広いので周波数分解能は落ちるがサイドローブが矩形と比べて小さい。 サイドロープがシャープに落ちているのでダイナミックレンジを広くとれる。 ただし、メインローブに近いサイドローブがあまり落ちていないので、 近接した周波数が含まれる信号に対しては向いていない。

win_hann
hamming

最初の節は$\frac{4\pi}{N}$。幅はhannと同じだが細い。hannと比較してメインローブのすぐ脇にでるサイドローブの大きさが 抑えられてるので周波数が近接した複数の信号に向いている。 ただ、サイドローブの減衰が鈍いので広いダイナミックレンジは望めない。

win_hamming
balckman

最初の節の位置は$\frac{6\pi}{N}$。幅がhannやhammingより広いがサイドローブがとても小さい。サイドローブの減衰がよいのでダイナミックレンジが広くとれ、 強い信号と弱い信号が周波数を隔てて存在する場合に有効。 ただ、メインローブの幅が広いので近接した信号の分離ができずに ひとつのピークになってしまう危険性がある。

win_blackman
balckman harris

blackmanと比べて少し太くなるがサイドローブがblackmanより小さくなる。

win_blackmanharris

上記の窓関数を重ね合わせて表示して、サイドローブの大きさやメインローブの幅の比較してみる。

win_all

窓関数の使い分け

窓関数は、サイドローブのレベルが低く、メインローブが細いほうが、望ましい。

  • サイドローブが低い=離れた周波数の信号の影響度が小さい → ダイナミックレンジが大きい
  • メインローブの帯域幅が狭い = 近い周波数に対する影響度が小さい → 周波数分解能が高い
ただ、両方は性質はトレードオフの関係にあるので、用途に合わせてどちらの特徴が重要であるかを見極めて窓関数を選ぶ必要がある。

こちらのサイトには以下のような指針が紹介されている。
  • 短時間の過渡信号 → 矩形
  • ピークの値を正確にしたい → フラット・トップ
  • 非常に近い帯域の信号をはっきり分離したい → カイザー
  • 汎用的、どういった信号か不明 → ハニング、最小 4 項ブラックマン・ハリス
  • とにかくダイナミックレンジを広くしたい → ドルフ・チェビシェフ、ガウス、DPSS でパラメータ調整
まだ、大きなポイント数でFFTでき周波数分解能を確保できる場合は、サイドローブを低く抑える性質をもった窓関数を選ぶ傾向にあると。 参考資料としてNIのこちらのホワイトペーパーがあげられている。
"""
Compare window function
"""
# -*- coding: utf-8 -*-
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import sys, os


N = 51
#window function
w_rec = np.ones(N)
w_hann = signal.hann(N)
w_hamming = signal.hamming(N)
w_blackman = signal.blackman(N)
w_blackmanharris = signal.blackmanharris(N)
w_dic = {'rectangle':w_rec, 'hann':w_hann, 
    'hamming':w_hamming, 'blackman':w_blackman, 'blackman_harris':w_blackmanharris}
#display window function
w_list = []
for k in w_dic:
    if k == 'rectangle':
        continue
    w_list.append(k)
    plt.plot(w_dic[k])
plt.legend(w_list)
plt.ylabel("Amplitude")
plt.xlabel("sample")
plt.title('compare window function shape')
plt.show()
plt.close('all')
#display frequency response of window function
fftpt=8192
f = lambda x : np.fft.fftshift(np.fft.fft(x, fftpt)/np.sqrt(fftpt))
freq_list = 2 * np.pi * np.linspace(-0.5, 0.5, fftpt)
#each 
for k in w_dic:
    fftdt_lin = f(w_dic[k])
    fftdt_log = 20*np.log10(np.abs(fftdt_lin))
    #fftdt = np.abs(fftdt_lin)
    plt.plot(freq_list, fftdt_log)
    plt.legend([k])
    plt.ylabel("Amplitude[dB]")
    plt.xlabel("bin")
    plt.title('frequency response of window function')
    plt.xlim([-1,1])
    plt.ylim([-170,0])
    plt.show()
#all
w_list = []
for k in w_dic:
    w_list.append(k)
    fftdt_lin = f(w_dic[k])
    fftdt_log = 20*np.log10(np.abs(fftdt_lin))
    plt.plot(freq_list, fftdt_log)
plt.legend(w_list)
plt.ylabel("Amplitude[dB]")
plt.xlabel("bin")
plt.title('compare frequency response of window function')
plt.xlim([-1,1])
plt.show()
plt.close('all')



GitHubを用いてチームで開発することを想定した手順を調べてみた。
(参照)

ローカルリポジトリと該当のリモートリポジトリを連携したのち、

1.ローカルリポジトリでリモートのmasterブランチに入る


git checkout master
コードを編集するときは必ず一度masterに戻ってからブランチを切り直すこと。


2.作業用のブランチをつくる


git checkout -b new-branch1
これでローカルにブランチが作成される。この時点ではローカルにのみnew-branch1が作成された状態で、リモートにはnew-branch1ブランチは作成されていない。

ちなみに
git branch -vv
でローカルのリポジトリとリモートのリポジトリの関係を確認することができる。
master 61d459a [origin/master]
* new-branch1 73f24c0

3.常にリモートのmasterに更新があるかを確認する


作業用のブランチを作成し作業を開始したのちもgit pullでリモートのmasterの更新を確認する。
git pullでmasterで更新があったときは、リベースしてmasterの更新をローカルで作業しているブランチに反映する。参照

4. 開発がおわったらリモートリポジトリへpush


開発が終わったらリモートリポジトリへpushする。現状ではローカルのnew-branch1ブランチはリモートのブランチと紐付いていない。直接masterにはpushしない。リモートリポジトリにブランチを作成し、そこにpushする。
git push -u origin new-branch1
-uオプションをつけるとリモートリポジトリに新たにブランチを作成しそこにpushしてくれる。
現状で-uオプションをつけずにpushすると、push先がないとエラーとなる。

git branch -vvで

master 61d459a [origin/master]
* new-branch1 73f24c0 [origin/new-branch1]

とリモートリポジトリに新たにブランチが作成されて、ローカルのブランチとひもづいていることがわかる。

ちなみに、ローカルブランチでコミット履歴はそのままリモートのブランチにも引き継がれるので、ローカルで開発中に個人的にコミットした内容は、pushまえに整理しておいたほうよい。

リモートにpushしたのちは、誰かのコードレビューをうけるので、レビューワが認識しやすいようにコミットのコメントやコミットの粒度を整理してからpushしたほうがよい。

5. pull requestを出す


Github上のリポジトリページからpull requestを出してレビューを受ける。
その後の手順は参照を参照



コミットは作業ログではない


細かすぎる作業メモのようなコミットにならないようなTipsが書かれている。
(参照)
そもそもはコミットするときから気をつける。
・不要なコミットは避ける
・amendを用いて作業内容を一つのコミットにまとめる
・rebaseでコミットログをまとめる


rebaseでコミットログをまとめる


git log —online
でコミットログを確認し、

git rebase -i HEAD~3
というようにまとめたいコミットログの数をHEAD~のあとに指定する。

VSのターミナルにviエディタが開き

pick cce1dc9 DB用クラス追加
pick 70b4379 タイポ修正
pick aebf72c  変数名修正

というように古いコミット順に表示される。

70b4379 タイポ修正及びpick aebf72c  変数名修正を
pick cce1dc9 DB用クラス追加にまとめてしまいたい場合

pick cce1dc9 DB用クラス追加
squash 70b4379 タイポ修正
squash aebf72c 変数名修正

とし、:x で実行する。

次に、コミットメッセージを入力するvi編集画面に切り替わるので
コミットメッセージを編集して :x を実行する。

そうすると「0b4379 タイポ修正」、「pick aebf72c 変数名修正」及び「cce1dc9 DB用クラス追加」のコミットが1つにまとまる。

ちなみにコミットハッシュは新しいものになる。(参照)


rebaseでコミットをまとめるとは、以前のコミットに新しいコミットをまとめていいく作業だからコミット日時が大きく戻ることになりかねない。

個人で作業しているときは問題ないが、複数人で同じブランチをさわっていて、リモートのブランチが更新されていた場合、実際は更新分を取り込んでいたとしても、pullしたときに日付が古いと言ってはじかれてしまう可能性はある。

だから、rebaseはコミットログをまとめるに万能ではないので、上記に書かれているように不要なコミットをしないように作業中から気をつけておくことが大切。

↑このページのトップヘ