時間周波数解析の基礎の基礎, FFT, STFT, 窓関数などについて

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

t1= 10
Ts = 0.01
fs = 1/Ts
t = np.arange(0,t1+Ts,Ts)

f0=1.00
f1=10

y = signal.chirp(t,f0=f0,t1=t1, f1=f1,method="linear")

FONTSIZE = 18
LINEWIDTH = 1

plt.figure()
plt.plot(t, y, linewidth=LINEWIDTH)
plt.grid(True)
plt.xlabel("Time, s", fontsize=FONTSIZE)
plt.ylabel("Amplitude", fontsize=FONTSIZE)
plt.xticks(fontsize=FONTSIZE)
plt.yticks(fontsize=FONTSIZE)

"""


ここでは、チャープ信号をFFTする
"""

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

t1= 10
Ts = 0.01
fs = 1/Ts
t = np.arange(0,t1+Ts,Ts)

f0=1.00
f1=10

y = signal.chirp(t,f0=f0,t1=t1, f1=f1,method="linear")

Y = np.fft.rfft(y)                     #rは実数のみを表す
f_range = np.fft.rfftfreq(len(y), d=Ts)  # frequency axis
Amp = np.abs(Y)

FONTSIZE = 18
LINEWIDTH = 1

# ---- Fig1: time signal ----
plt.figure()
plt.plot(t, y, linewidth=LINEWIDTH)
plt.grid(True)
plt.xlabel("Time, s", fontsize=FONTSIZE)
plt.ylabel("Amplitude", fontsize=FONTSIZE)
plt.xticks(fontsize=FONTSIZE)
plt.yticks(fontsize=FONTSIZE)


plt.figure()
plt.plot(f_range, Amp, linewidth=LINEWIDTH)
plt.grid(True)
plt.xlabel("Frequency, Hz", fontsize=FONTSIZE)
plt.ylabel("Amplitude", fontsize=FONTSIZE)
plt.xticks(fontsize=FONTSIZE)
plt.yticks(fontsize=FONTSIZE)


plt.show()

"""

ここでは、チャープ信号を矩形窓を使って短時間フーリエ変換を行う
"""

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

t1= 10
Ts = 0.01
fs = 1/Ts
t = np.arange(0,t1+Ts,Ts)

f0=1.00
f1=10

y = signal.chirp(t,f0=f0,t1=t1, f1=f1,method="linear")

Y = np.fft.rfft(y)                   
f_range = np.fft.rfftfreq(len(y), d=Ts)  # frequency axis
Amp = np.abs(Y)


win_size = 256
inc = 1
nfft = 2048
noverlap = win_size-inc

# Rectangle window
win_rec = signal.windows.boxcar(win_size)
"""
detrend str or function or False, optional
Specifies how to detrend each segment. If detrend is a string, it is passed as the type argument to the detrend function. If it is a function, it takes a segment and returns a detrended segment. If detrend is False, no detrending is done. Defaults to False.

return_onesidedbool, optional
If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.

boundarystr or None, optional
Specifies whether the input signal is extended at both ends, and how to generate the new values, in order to center the first windowed segment on the first input point. This has the benefit of enabling reconstruction of the first input point when the employed window function starts at zero. Valid options are ['even', 'odd', 'constant', 'zeros', None]. Defaults to ‘zeros’, for zero padding extension. I.e. [1, 2, 3, 4] is extended to [0, 1, 2, 3, 4, 0] for nperseg=3.

paddedbool, optional
Specifies whether the input signal is zero-padded at the end to make the signal fit exactly into an integer number of window segments, so that all of the signal is included in the output. Defaults to True. Padding occurs after boundary extension, if boundary is not None, and padded is True, as is the default.


"""
f,tt,zxx = signal.stft(y,fs=fs,window=win_rec,nperseg=win_size,noverlap=noverlap,nfft=nfft,boundary=None)



FONTSIZE = 18
LINEWIDTH = 1

# ---- Fig1: time signal ----
plt.figure()
plt.plot(t, y, linewidth=LINEWIDTH)
plt.grid(True)
plt.xlabel("Time, s", fontsize=FONTSIZE)
plt.ylabel("Amplitude", fontsize=FONTSIZE)
plt.xticks(fontsize=FONTSIZE)
plt.yticks(fontsize=FONTSIZE)


plt.figure()
plt.plot(f_range, Amp, linewidth=LINEWIDTH)
plt.grid(True)
plt.xlabel("Frequency, Hz", fontsize=FONTSIZE)
plt.ylabel("Amplitude", fontsize=FONTSIZE)
plt.xticks(fontsize=FONTSIZE)
plt.yticks(fontsize=FONTSIZE)


# ---- Fig3: Spectrogram (Rectangle) ----
plt.figure()
plt.pcolormesh(tt, f, 20 * np.log10(np.abs(zxx) + 1e-12), shading="auto")
cb = plt.colorbar()
cb.set_label("Magnitude, dB", fontsize=FONTSIZE)
plt.xlabel("Time, s", fontsize=FONTSIZE)
plt.ylabel("Frequency, Hz", fontsize=FONTSIZE)
plt.xticks(fontsize=FONTSIZE)
plt.yticks(fontsize=FONTSIZE)

plt.show()

"""
ここでは、STFT結果から振動周波数を得る
"""

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

t1= 10
Ts = 0.01
fs = 1/Ts
t = np.arange(0,t1+Ts,Ts)
t_ref = np.arange(0,10,0.01)
y_ref = 9/10*t_ref+1
f0=1.00
f1=10

y = signal.chirp(t,f0=f0,t1=t1, f1=f1,method="linear")

Y = np.fft.rfft(y)                     
f_range = np.fft.rfftfreq(len(y), d=Ts)  
Amp = np.abs(Y)


win_size = 256
inc = 1
nfft = 2048
noverlap = win_size-inc

# Rectangle window
win_rec = signal.windows.boxcar(win_size)

f,tt,zxx = signal.stft(y,fs=fs,window=win_rec,nperseg=win_size,noverlap=noverlap,nfft=nfft,boundary=None)
"""zxx.shape = (周波数ビン数, 時間フレーム数)"""
P_rec = np.abs(zxx) ** 2
kmax_rec = np.argmax(P_rec, axis=0)   # 各時間フレームごとに 最大周波数ビンを求めている
f_peak_rec = f[kmax_rec]


FONTSIZE = 18
LINEWIDTH = 1

plt.figure()
plt.plot(tt,f_peak_rec,label="peak frequency")
plt.plot(t_ref,y_ref,label="chirp")
plt.grid()
plt.legend()
plt.xlabel("Time, s",fontsize=FONTSIZE)
plt.ylabel("Frequency, Hz",fontsize=FONTSIZE)
plt.xticks(fontsize=FONTSIZE)
plt.yticks(fontsize=FONTSIZE)



plt.show()

win_rec = signal.windows.boxcar(win_size)

win_rec = signal.windows.blackman(win_size)