如何通過計算裝倉窗口的功率譜來繪製信號的頻譜圖？

Question

在這裏，我可以生成一個信號：

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() 

在這裏，我可以繪製整個信號的功率譜：

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() 

接下來我要生成頻譜圖，表示爲頻率（y軸）和時間（x軸）的圖像，但我是新的進行傅立葉分析，並在這個階段對如何使用window function（矩形，漢明，漢寧等）感到困惑。有沒有適當的方法來做到這一點，以便我可以使用我選擇的窗口函數及時分解信號？

Answer 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() 

Answer 2

，或者您可以使用matplotlib.pyplot.specgram看到：http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.specgram