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

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

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

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

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

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

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

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

一方、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 
"""
# -*- 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()