日に日に分からんことが増えていく…

φ(..)メモメモ

信号解析のフィルタリング

Python3

import scipy as sp
print(sp.__version__)

>>> 0.19.1

信号の作成

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

N = 1024            # サンプル数
dt = 0.001          # サンプリング周期 [s]
f1, f2, f3 = 10, 60, 300 # 周波数 [Hz]

t = np.arange(0, N*dt, dt) # 時間 [s]
x = 3*np.sin(2*np.pi*f1*t) + 0.3*np.sin(2*np.pi*f2*t) + 0.2*np.sin(2*np.pi*f3*t) # 信号

xn = x + np.random.randn(len(t)) * 0.3

フーリエ解析 (周波数の導出)

np.fft.fftfreq で周波数を抽出
np.fft.fft 高速フーリエ変換

F = np.fft.fft(x) # 変換結果
freq = np.fft.fftfreq(N, d=dt) # 周波数

Amp = np.abs(F/(N/2)) # 振幅

freq =  np.fft.fftfreq(len(x), d=dt)

fig, ax = plt.subplots()
ax.plot(freq[1:int(N/2)], Amp[1:int(N/2)])
ax.set_xlabel("Freqency [Hz]")
ax.set_ylabel("Amplitude")
ax.grid()
plt.show()

信号のフィルタリング

signal.lfilter はただの式 scipy のドキュメント参照。
a=1 だけなので、FIR によりフィルタリング。いい感じに low filtering などをしてくれる係数を出してくれる。
signal.firwin はフィルター用の窓関数を抽出。

プロット関数

def plot_show(x, y):
    sns.set()
    plt.figure(figsize=(14, 4))
    plt.plot(x, y)
    plt.margins(0)
    plt.show()

# ナイキスト周波数

1/(dt*2)

>> 500.0

フィルタリング

pass_zero= 1 の場合は、high pass or band pass

high pass

filter1 = signal.firwin(numtaps=255, cutoff=40, nyq=1/(dt*2), pass_zero=0)
y1 = signal.lfilter(filter1, 1, xn)

F1 = np.fft.fft(y1)
Amp1 = np.abs(F1/(N/2))

plot_show(t, y1)

low pass

filter1 = signal.firwin(numtaps=255, cutoff=60, nyq=1/(dt*2))
y1 = signal.lfilter(filter1, 1, xn)

F1 = np.fft.fft(y1)
Amp1 = np.abs(F1/(N/2))

plot_show(t, y1)

band pass

filter1 = signal.firwin(numtaps=255, cutoff=[40, 70], nyq=1/(dt*2), pass_zero=0)
y1 = signal.lfilter(filter1, 1, xn)

F1 = np.fft.fft(y1)
Amp1 = np.abs(F1/(N/2))

plot_show(t, y1)

band eliminate

filter1 = signal.firwin(numtaps=101, cutoff=[100, 200], nyq=1/(dt*2), pass_zero=1)
y1 = signal.lfilter(filter1, 1, xn)

F1 = np.fft.fft(y1)
Amp1 = np.abs(F1/(N/2))

plot_show(t, y1)

DWT での denoise

def maddest(d, axis=None):
    """
    Mean Absolute Deviation
    """
    return np.mean(np.absolute(d - np.mean(d, axis)), axis)


def denoise_signal( x, wavelet='db4', level=1):
    """
    1. Adapted from waveletSmooth function found here:
    http://connor-johnson.com/2016/01/24/using-pywavelets-to-remove-high-frequency-noise/
    2. Threshold equation and using hard mode in threshold as mentioned
    in section '3.2 denoising based on optimized singular values' from paper by Tomas Vantuch:
    http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf
    """
    
    # Decompose to get the wavelet coefficients
    coeff = pywt.wavedec( x, wavelet, mode="per" )
    
    # Calculate sigma for threshold as defined in http://dspace.vsb.cz/bitstream/handle/10084/133114/VAN431_FEI_P1807_1801V001_2018.pdf
    # As noted by @harshit92 MAD referred to in the paper is Mean Absolute Deviation not Median Absolute Deviation
    sigma = (1/0.6745) * maddest( coeff[-level] )

    # Calculte the univeral threshold
    uthresh = sigma * np.sqrt( 2*np.log( len( x ) ) )
    coeff[1:] = ( pywt.threshold( i, value=uthresh, mode='hard' ) for i in coeff[1:] )
    
    # Reconstruct the signal using the thresholded coefficients
    return pywt.waverec( coeff, wavelet, mode='per' )


def denoise_signal_raw( x, wavelet='db4', level=1):
    
    # Decompose to get the wavelet coefficients
    coeff = pywt.wavedec( x, wavelet, mode="per" )
    return pywt.waverec( coeff[: level + 1], wavelet, mode='per' )

CWT から scaleograms の作成

def plot_wavelet(time, signal, scales, 
                 waveletname = 'cmor', 
                 cmap = plt.cm.seismic, 
                 title = 'Wavelet Transform (Power Spectrum) of signal', 
                 ylabel = 'Frequence (Hz)', 
                 xlabel = 'Time'):
#     from IPython.core.debugger import Pdb; Pdb().set_trace()
    [coefficients, frequencies] = pywt.cwt(signal, scales, waveletname, dt)
    power = (abs(coefficients)) ** 2
    
    z = np.log2(power)
#     z = (z - z.mean()) / z.std()
#     z= np.clip(z, -1, 1)
    
    period = 1. / frequencies
    
    levels = np.linspace(-10, 10, 99)
    print(z.std(), z.mean())
    
    fig, ax = plt.subplots(figsize=(15, 10))
    print(frequencies)
    im = ax.contourf(time, np.log2(frequencies), z, extend='both', levels=levels, cmap=cmap)
    
    ax.set_title(title, fontsize=20)
    ax.set_ylabel(ylabel, fontsize=18)
    ax.set_xlabel(xlabel, fontsize=18)
    
    yticks = 2**np.arange(np.ceil(np.log2(frequencies.min())), np.ceil(np.log2(frequencies.max())))
    print(yticks)
#     yticks = frequencies

    ax.set_yticks(np.log2(yticks))
#     ax.set_yticks(yticks)
    
    ax.set_yticklabels(yticks)
    ax.invert_yaxis()
    ylim = ax.get_ylim()
    ax.set_ylim(ylim[1], ylim[0])
    
    cbar_ax = fig.add_axes([0.95, 0.5, 0.03, 0.25])
    fig.colorbar(im, cax=cbar_ax, orientation="vertical")
    plt.show()