何かあれば GitHub のリポジトリに issue を作るか ryukau@gmail.com までお気軽にどうぞ。
Update: 2023-10-30
注意: 内容が怪しいです。
GottleibとShuの論文 “ON THE GIBBS PHENOMENON AND ITS RESOLUTION” を基にギブス現象を抑える方法を試します。
この文章で実装したコードをGitHubで見ることができます。実行するにはSciPy、pyFFTW、matplotlibが必要です。パラメータを変えて遊んでみてください。
GottleibとShuの論文で先行研究として紹介されていた方法です。ここではFIRフィルタとして実装します。
フィルタを適用する部分を実装します。
import numpy
import matplotlib.pyplot as pyplot
from pyfftw.interfaces.numpy_fft import rfft, irfft
from scipy.signal import *
from scipy.special import (factorial, eval_gegenbauer, gamma)
def spectralLowpass(sig, numSeries):
if numSeries >= len(sig):
return sig
= rfft(sig)
spec len(sig) - 1] = 0
spec[numSeries:return irfft(spec)
def additiveSaw(length, numSeries):
return spectralLowpass(numpy.linspace(-1, 1, length), numSeries)
def applyFilter(source, filterFunc, power=1):
= numpy.zeros(len(source))
filt = (len(source) - 1) / 2
half for i in range(len(filt)):
= (i - half) / half
eta = filterFunc(eta)**power
filt[i] /= numpy.sum(filt)
filt
return (convolve(source, filt, mode='full'), filt)
def plotGibbsSuppression(length, filterFunc, numSeries=16, power=1):
= additiveSaw(length, numSeries)
source = applyFilter(source, filterFunc, power)
suppressed, filt = numpy.append(source, numpy.zeros(length))
source # 結果をプロットする。ここでは省略。
フィルタの適用は plotGibbsSuppression
から行います。
power
の値を大きくするとフィルタのカットオフ周波数が高くなります。入力信号は低いほうから
numSeries
個の周波数成分を持つ鋸歯波です。
# Lanczosフィルタの適用例。
1024, lanczos, 16, 1) plotGibbsSuppression(
\[ \sigma_1(\eta) = 1 - \eta \]
実装します。
def fejer(eta):
return 1 - eta
フィルタを適用した結果です。画像のSignalは入力信号、Resultはフィルタを適用した結果です。一番下の図はFIRフィルタの係数です。
\[ \sigma_2(\eta) = \frac{\sin(\pi\eta)}{\pi\eta} \]
実装します。
def lanczos(eta):
return numpy.sinc(eta)
フィルタを適用した結果です。
\[ \sigma_3(\eta) = \frac{1}{2}(1 + \cos(\pi\eta)) \]
実装します。
def raisedCosine(eta):
return (1 + numpy.cos(numpy.pi * eta)) / 2
フィルタを適用した結果です。
\[ \sigma_4(\eta) = \sigma_3^4(\eta)(35 - 84\sigma_3(\eta) + 70\sigma_3^2(\eta) - 20\sigma_3^3(\eta)) \]
実装します。
def sharpendRaisedCosine(eta):
= raisedCosine(eta)
s return s**4 * (35 - 84 * s + 70 * s**2 - 20 * s**3)
フィルタを適用した結果です。
\[ \sigma_5(\eta) = e^{-\alpha\eta^p} \]
実装します。
def exponential(eta, alpha=2**10, p=4):
return numpy.exp(-alpha * eta**p)
フィルタを適用した結果です。
\[ \sigma_6(\eta) = 1 - \frac{(2p-1)!}{(p-1)!} \int_0^{\eta} [t(1-t)]^{p-1} dt \]
Maximaで展開します。 \(p\) は任意の正の実数です。
declare(eta, real, t, real);
p: 4;
I: integrate((t*(1 - t))**(4 - 1), t, 0, eta);
expand(1 - (2 * p - 1)! / (p - 1)! * I);
適当に \(p = 4\) として実装します。
def daubechies(eta):
return 1 - 210 * eta**4 + 504 * eta**5 - 420**eta**6 + 120 * eta**7
フィルタを適用した結果です。
GottleibとShuの論文で提案されていた方法です。この方法はギブス現象が発生する前の信号が解析的な形で表現できるときに使えます。解析的な形とは \(f(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \dots\) というような形のことです。
まずGegenbauer係数 \(\hat{g}_k^{\lambda}\) を求めます。
\[ \hat{g}_k^{\lambda} = \frac{1}{h_k^{\lambda}} \int_{-1}^{1} (1 - x^2)^{\lambda - \frac{1}{2}}f_N(x)C_k^{\lambda}(x)\,dx ,\quad 0 \le k \le m \]
\(h_k^{\lambda}\) は次の式で求められます。
\[ h_k^{\lambda} = \pi^{\frac{1}{2}}C_k^{\lambda}(x) \frac{\Gamma(\lambda + \frac{1}{2})}{\Gamma(\lambda)(k + \lambda)} \]
Gegenbauer係数 \(\hat{g}_k^{\lambda}\) を用いて出力信号 \(g^m(x)\) を合成します。
\[ g^m(x) = \sum_{k = 0}^m \hat{g}_k^{\lambda} C_k^{\lambda}(x) \]
\(g^m(x)\) は真の信号と比例関係にあるそうです。つまり出力信号の振幅を調整する必要があります。
\(\lambda = N / 4\) 、 \(m = N / 2\) として実装します。
def gottliebShu(sig, N, gam=0.25):
"""
:param sig: numpy.array input signal.
:param N: Number of spectral partial sum or number of overtone.
:param gam: Positive real constant.
"""
= gam * N # N が大きいと eval_gegenbauerでオーバーフローする。
lam
= numpy.linspace(-1.0, 1.0, len(sig))
x = numpy.arange(0, int(N / 2))
k = eval_gegenbauer(numpy.tile(numpy.vstack(k), len(sig)), lam, x)
ggnbr = numpy.transpose(ggnbr)
ggnbr_t
= numpy.sqrt(numpy.pi) \
h_k_lam * ggnbr_t[-1] \
* gamma(lam + 0.5) / gamma(lam) / (k + lam)
= numpy.sum(
g_hat 1 - x * x, numpy.abs(lam - 0.5)) * sig * ggnbr,
numpy.power(=1) / h_k_lam
axis
return (numpy.sum(g_hat * ggnbr_t, axis=1), g_hat)
適用した結果です。画像の一番下の図は \(\hat{g}_k^{\lambda}\) です。
\(C_1^{\lambda}\) の成分で完璧に再現できています。 \(C_1^{\lambda}\) の成分は \(y = ax\) の形の直線です。入力信号の鋸歯波も直線なので、これはどうやっても成功しそうです。
他の信号でも試してみます。
def additiveNoise(length, numSeries):
return spectralLowpass(numpy.random.uniform(-1, 1, length), numSeries)
def analyticSignal(length, numSeries):
= numpy.linspace(-1, 1, length)
sig = numpy.zeros(length)
noise = 128
order = numpy.random.random(order)
rand for i in range(order):
+= rand[i] * sig**(i + 1)
noise return spectralLowpass(noise, numSeries)
def plotGottleibShuNoise(length, signalFunc, numSeries=16):
= signalFunc(length, numSeries)
source = gottliebShu(source, numSeries)
suppressed, g_hat = spectralLowpass(suppressed, numSeries)
regibbsed # 結果をプロットする。ここでは省略。
1024, additiveNoise, 16)
plotGottleibShuNoise(1024, analyticSignal, 16) plotGottleibShuNoise(
additiveNoise
の結果です。 additiveNoise
は乱数で生成したノイズにローパスフィルタをかけた信号です。画像のFiltered
resultはギブス現象を抑えた結果にフィルタをかけて、周波数成分の数を入力信号と揃えた信号です。
additiveNoise
について入力信号と、ギブス現象を抑えた結果にフィルタをかけた信号の周波数の比較です。powerは周波数成分の大きさ、phaseは周波数成分の位相です。
additiveNoise
ではギブス現象の抑制に失敗しています。乱数で生成したノイズは不連続点を含むのでGegenbauer
Polynomialを使う方法をそのまま使うことはできません。
analyticSignal
の結果です。 analyticSignal
は解析的な要素のみで合成した信号です。
analyticSignal
について入力信号と、ギブス現象を抑えた結果にフィルタをかけた信号の周波数の比較です。
analyticSignal
では理論どおりにギブス現象の抑制に成功していると言えそうです。ただし今回選んだパラメータでは目に見える誤差が出ています。
\(m\)
の値を大きくすれば誤差を減らすことができます。
\(\hat{g}_k^{\lambda}\) を求める式の中にある \(f_N(x)\) は \([-1, 1]\) の区間で連続であることが求められています。つまり、この方法を現実の信号に適用するには不連続点を見つけて、連続な区間ごとに信号を分割する必要があります。
\(\lambda\) と \(m\) についてのパラメータ調整が必要です。
\(\lambda\) と \(m\) が大きいと Gegenbauer Polynomial の計算でオーバーフローすることがあります。これが問題になるときは任意精度で実装するしかなさそうです。
Gegenbauer polynomialを使う方法の \(m\) を増やしていくことで、Gegenbauer polynomialによる信号の加算合成ができます。ただし信号の両端の値が非常に大きくなるのでTukey窓などを掛け合わせて抑える必要があります。また \(\lambda = 1\) とすることで、ある程度はオーバーフローを避けることができます。