2016-03-27 2 views
2

ここで私は、信号を生成することができます:私は生成したいビニングされたウィンドウのパワースペクトルを計算することによって、信号のスペクトログラムをプロットするにはどうすればよいですか?

time_step = 1/fs 
ps = np.abs(np.fft.fft(x))**2 
freqs = np.fft.fftfreq(x.size, time_step) 
idx = np.argsort(freqs) 
plt.plot(freqs[idx], 256*ps[idx]/max(ps[idx])) # set max to 256 for later image plotting purposes 
plt.xlabel('frequency') 
plt.ylabel('power') 
plt.show() 

enter image description here

次へ:私は、全体の信号のパワースペクトルをプロットすることができます。ここ

import numpy as np 
from matplotlib import pyplot as plt 
from numpy.lib import stride_tricks 
import seaborn as sns 
sns.set(style = "darkgrid") 

fs = 48000.0 
t = np.arange(0, 10, 1.0/fs) # 0 to 10 sec at 48k samples per second 
f0 = 1000 
phi = np.pi/2 # pi/2 

x = 0 # initial x 
f = [500, 100, 40, 1] #vector of frequencies 
A = [1, 0.5, 0.25, 0.1] #vector of amplitudes 
for i in range(0, len(f)): 
    x = x + A[i] * np.sin(2 * np.pi * f[i] * t + phi) #add waves 
x = x + max(x) # shift plot upwards 
plt.plot(t, x) 
plt.axis([0, .05, 0, max(x)]) 
plt.xlabel('time') 
plt.ylabel('amplitude') 
plt.show() 

enter image description here

を周波数(y軸)と時間(x軸)の画像として表されるスペクトログラムですが、私は新しいこの段階でwindow function(長方形、ハミング、ハニングなど)を使用する方法について混乱しています。これを行う適切な方法がありますので、自分の選択したウィンドウ機能を使用して時間内に信号を分割することができます。

答えて

1

この追加:

M = 5000 
overlap = 500 
unique = M - overlap 
han = np.hanning(M) 
f_border = 2*max(f) 

for i in range(0, x.shape[0], unique): 
    if i + M > x.shape[0]: 
     break 
    curr_x = x[i:i+M] 
    y = 10*np.log10(np.abs(np.fft.fft(curr_x*han))**2) 
    if i == 0: 
     freqs = np.fft.fftfreq(curr_x.size, time_step) 
     idx = np.argsort(freqs) 
     freqs = freqs[idx] 
     idx2 = np.where(np.logical_and(freqs > 0, freqs < f_border))[0] 
    y = y[idx][idx2][np.newaxis].T 
    try: 
     stereogram = np.hstack([stereogram, y]) 
    except NameError: 
     stereogram = y 

fig = plt.figure() 
ax = fig.add_subplot(111) 
ax.imshow(stereogram) 
yticks = ax.get_yticks()[1:-1] 
plt.yticks(yticks, (yticks * f_border/yticks[-1]).astype('str')) 
plt.ylabel('frequency') 
plt.xlabel('time') 
plt.show() 

enter image description here

関連する問題