カテゴリ: ひとり学習

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

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

ゼロ点を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')

↑このページのトップヘ