リサンプリングの方法としてゼロ点を挿入・サンプル点を間引く方法と、フーリエ変換を使う方法の二種類が考えられる。
ゼロ点挿入・サンプル点間引き
N倍にアップサンプリングすることを考える。
まず、N-1のゼロを元信号のサンプル点の間に挿入して(インターポレーション)、サンプリング周波数をあげる。
以下の図は、3倍にアップサンプリングする場合の例。元信号のサンプル点の間に2個ゼロ点を挿入している。
元信号
サンプル点の間にゼロ点を2個入れた状態

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

そこで、LPFで不要な繰り返しの部分をカットしてあげることで、元信号のサンプリング周波数をN倍した信号が得られる。
LPFで繰り返しの部分をカット
元信号とアップサンプリングした信号を重ねて表示してみる。ほぼ重なっている。アップサンプリングすることで曲線がなめらかになっている分、元信号が少しカクカクにみえる。

一方、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個のゼロ点を挿入する。
元信号をフーリエ変換した状態
中央にゼロ点を挿入しサンプル数を64から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()
コメント