Signal Processing : Temporal / Spectral Space (1)

Temporal / Spectral Space

 

필터링을 이해하려면 신호를 주파수 구성 요소로 보는 것이 좋다.

 

Fourier Series

 

푸리에 급수는 일련의 조화롭게 연결된 (모든 주파수는 가장 낮은 주파수의 배수) 정현파로, 일단 가중되고 합산되면 기본 주기 (신호에 존재하는 가장 낮은 주파수) 동안 임의의 신호와 동일하다. 또는 신호가 결정된 길이인 경우 해당 신호의 길이보다 1배 더 많다.

 

이를 설명하기 위해 듀티 사이클이 50%인 440Hz에서 구형파 신호를 분해하면, 신호의 정확히 절반이 양수이고 절반이 음수인 경우 초당 +1에서 -1까지 440번 진동하는 주기적 신호이다.

 

## load libraries
import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np
import cmath as cm
from scipy import signal as sps
from IPython.display import HTML

#This is the fundamental frequency. This is also the frequency of the first component.
fundamental = 440 
# Compute time and waveform data points
t = np.linspace(0,2/fundamental,1000);
x = sps.square(2*np.pi*440*t,0.5);
# show the waveform
fig, ax = plt.subplots(1,1,figsize=(12,4),sharex='col');

# Make the axes look nice
ax.axis([0,2.1/fundamental,-1.5,1.5]);
ax.set_xlabel('time [s]')
ax.set_ylabel('amplitude')
ax.set_xticks(np.arange(3)/fundamental)

# Draw the waveform
l1, = ax.plot(t,x,lw=4,color='k');

 

구형파에 대한 푸리에 급수 f(t)는 다음과 같이 계산된다.

 

 

이것은 시리즈가 모든 홀수 구성 요소로 구성되어 있음을 의미한다. (n = 1 , 3 , 5 , … , ∞) 여기서 각각을 1/n로 곱한다.

 

(구성 요소는 점점 작아짐) 구성 요소의 합은 다음과 같이 곱해진다. 4/π (즉, 1보다 약간 높음). 이 공식은 신호를 푸리에 계열로 분석적으로 분해할 수 있는 푸리에 변환에서 파생된다. 실제 세계에서 FFT (고속 푸리에 변환)와 같은 응용 알고리즘이 사용되며 이 알고리즘은 이산 (디지털) 푸리에 변환을 근사하는 데 매우 적합하다.

 

# Compute The Fourier series components
# Nframes determines how many frames the animation will have 
# each frame will add one more component, hence the number of
# frames is also the number of components to be computed
Nframes = 10;
# Set up empty line to accept the series' data
l2, = ax.plot([],[],lw=4,color='b');
# Set up empty lines to accept the components data
lines = [ax.plot([],[]) for _ in range(Nframes)];

# Add textual information on root mean square deviation
# and number of components
text = ax.text(0.0002,-1,"RMSD = {:.3}\n # of components = {:}".format(np.sqrt(np.mean(np.square(x))),0));

# Compute the components according to equation above
x_components = [4/np.pi*1/n*np.sin(2*np.pi*fundamental*n*t) for n in range(1,Nframes*2+1,2)];

# Define animation update function
def animate(i):
    # at frame 0, remove all components
    if i==0:
        for line in lines:
            line[0].set_data([],[]);
    #recompute the fourier series to the ith component
    x_series = sum(x_components[0:(i+1)],0);
    # plot series
    l2.set_data(t,x_series);
    # add new component to graph
    lines[i][0].set_data(t,x_components[i]);
    # update text
    x_rmsd = np.sqrt(np.mean(np.square(x-x_series)));
    text.set_text("RMSD = {:.3}\n # of components = {:>}".format(x_rmsd,i+1));
    
# Create animation object
ani = matplotlib.animation.FuncAnimation(fig, animate, frames=Nframes, interval = 500);
# Run animation
HTML(ani.to_jshtml())

 

구성 요소의 수가 증가할수록 시리즈는 원본에 점점 더 가까워진다. 평균 제곱근 편차 (맞춤 품질 의 지표, 낮을수록 좋음)는 로그 방식으로 감소한다.

 

Spectral representation

 

신호의 스펙터 (또는 스펙트럼)는 주파수 signature이다. 일반적으로 이것은 전력 스펙트럼 밀도 (PSD) 형식으로 표시된다. 주파수 대 전력 밀도 플롯이다. 주어진 주파수에서 스펙터의 진폭은 해당 주파수에 포함된 전력의 양을 알려준다.

 

# Initialize variables
# frequency of the periodic signals
f_signal = 4
# sampling rate
sr = 68

# timeseries vector of length 4*T (two periods(T), T = 1/f)
t = np.arange(0,4*1/f_signal,1/sr)

 

scipy 각 신호는 이론적인 표현이므로 한 번은 signal.periodogram 함수로 계산된 PSD를 사용하여 한 번은 왼쪽의 시간 공간에, 스펙트럼 공간에는 두 번 표시된다.

 

Dirac delta function (δ)

 

δ 신호 처리에 있어 매우 중요한 기능이다. 이것은 전자 시스템을 특성화하는데 도움이 되는 소위 임펄스 응답에 대한 입력이다. 또한 일반적인 신호의 스펙트럼에도 자주 나타난다.

 

디랙 델타 함수는 다음과 같이 정의된다.

 

 

즉, 진폭이 무한대이고 폭이 0이고 표면적이 1인 스파이크이다. 이것은 물론 디지털 세계에서 사용할 수 없으므로 다음과 같이 표현된다. 

 

 

n은 샘플 번호이고 sr은 샘플 속도이다. 이것은 진폭 sr 및 너비의 사각 펄스를 생성한다.

1/sr은

곡선 아래의 면적을 제공한다.

 

# Initialize figure
from scipy.signal import periodogram
fig, ax = plt.subplots(1,3,figsize=(12,4));
# Direct current(DC)
x = [0 if x !=0 else sr for x in t]
f = [0,32]
X = [1,1]
fp,Xp = periodogram(x,fs=sr,detrend=False)

ax[0].axis([0,1,0,1.1])
ax[0].plot(t,x)
ax[1].axis([-1,sr/2,10**-2,3])
ax[1].plot(f,X)
ax[1].set_yscale('log')
ax[2].axis([-1,sr/2,10**-2,3])
ax[2].plot(fp,Xp)
ax[2].set_yscale('log')

주기도의 진폭은 일반적으로 크기 차이의 차수에 관심이 있기 때문에 일반적으로 로그 스케일로 표현된다. 주파수 축 범위는 0에서 sr/2 이다. 이것은 Nyquist의 주파수 이하의 모든 것을 다루기 때문에 표시하기에 편리한 양의 데이터이다. 그것은 유일한 이유가 아니라, 그 이상을 사용하는 것과도 관련이 있다. scipy.signal.periodogram.

 

이론적 PSD는 모든 주파수에서 1의 일정한 진폭이다. 주기도가 진폭 1에서 2까지 약간 벗어났음을 알 수 있다. 이는 작업할 데이터 창이 매우 짧기 때문이다. δ의 범위는 −∞에서 +∞까지 있어햐 한다. 그러나 대부분 크기 차이의 순서에 관심이 있으므로 2의 요소가 그렇게 많이 떨어지지 않다.

 

https://neuro.inf.unibe.ch/AlgorithmsNeuroscience/Tutorial_files/Temporal_vs_Spectral.html