# Spectrogram_python

Spectrogram(시간, 주파수) 구하는 파이썬 코드에 대해 기술한다.

 

[50Hz, 20Hz 짜리 합성된 신호(z)]

- Sampling rate(Fs) : 1초 당 2000개 

#-*- coding:utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

Fs = 2000.
Ts = 1 / Fs
te = 1.0
t = np.arange(0.0, te, Ts)

# Signal x (20Hz) + Signal y (50Hz)
x = np.cos(2 * np.pi * 20 * t)
y = np.cos(2 * np.pi * 50 * t)

# Signal z
z = x + y

N = len(z)

k = np.arange(N)
T = N / Fs
freq = k / T
freq = freq[range(int(N/2))]

# FFT 적용
yfft = np.fft.fft(z)
yf = yfft / N
yf = yf[range(int(N/2))]

plt.rcParams["figure.figsize"] = (15,4)

# FFT 출력
plt.plot(freq, abs(yf), 'b')
plt.xlabel('Frequency')
plt.ylabel('Amplitude')
plt.xlim(0, Fs / 20)
plt.show()

FFT 출력

 

FFT 결과를 보면 주파수 성분 50Hz, 20Hz가 존재하는 것을 확인 할 수 있다.

 

하지만 이 그림에서는 어느 시점에 해당 주파수 성분이 존재하는지 알 수 없다.

따라서, STFT (고정된 시간 길이로 분할 및 FFT 수행)을 통해 해당 시점에 주파수성분을 알 수 있게 된다.

 

 

[신호(z)에 대한 Spectrogram 분석]

Spectgrogram 구하는 4가지 방법 

- pylab

- matplotlib.specgram

- scipy의 signal

- librosz의 stft

import pylab 
from scipy import signal


## Spectrogram type 1
f, tt, Sxx = signal.spectrogram(z, Fs)

plt.pcolormesh(tt, f, Sxx, shading='gouraud')
plt.title('Spectrogram type 1')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.ylim(0, 200)
plt.show()
    
## Spectrogram type 2
pp, ff, tt, pplot = pylab.specgram(
z,
NFFT=256,
Fs=2000,
detrend=pylab.detrend_none, 
window=pylab.window_hanning,
noverlap=int(128 * 0.5))
pylab.title('Spectrogram type 2')
pylab.ylim(0, 200)


## Spectrogram type 3
fig, (ax1, ax2) = plt.subplots(nrows=2)

ax1.plot(t, z)
ax1.set_xlabel('Time(s)')
ax1.set_ylabel('Amplitude')
fig.tight_layout()

Pxx, freqs, bins, im = ax2.specgram(z, NFFT=256, Fs=Fs, noverlap=int(128 * 0.5))
ax2.set_title('Spectrogram type 3')
ax2.set_xlabel('Time(s)')
ax2.set_ylabel('Frequency')
fig.tight_layout()

plt.show()

 

signal.spectrogram

- 해당 결과를 y축 0~200으로 확대시키면 다음과 같이 50Hz, 20Hz에서의 변화를 감지할 수 있다.

pylab.specgram

plt.specgram

원 z 신호(위)와 스펙트로그램(아래)

 

librosa.stft

- 파이썬 STFT를 사용한 Spectrogram

import librosa.display

# STFT power spectrum
D_octave = np.abs(librosa.stft(z))
librosa.display.specshow(librosa.amplitude_to_db(D_octave, ref=np.max), sr=Fs, y_axis='linear', x_axis='time')
plt.title('Spectrogram type 4')
plt.ylim(0, 200)
plt.show()