栏目分类:
子分类:
返回
名师互学网用户登录
快速导航关闭
当前搜索
当前分类
子分类
实用工具
热门搜索
名师互学网 > IT > 软件开发 > 后端开发 > Python

220430-三种短时傅立叶变换STFT的绘制及数据提取

Python 更新时间: 发布时间: IT归档 最新发布 模块sitemap 名妆网 法律咨询 聚返吧 英语巴士网 伯小乐 网商动力

220430-三种短时傅立叶变换STFT的绘制及数据提取

  • Figure

  • Code
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal
# - [How to plot Spectrogram using STFT in python?](https://stackoverflow.com/questions/43109708/how-to-plot-spectrogram-using-stft-in-python/48205933#48205933)

sampling_rate = 8000.0
fft_size = 1024
step = fft_size/16
time = 2
t = np.arange(0, time, 1/sampling_rate)
sweep = signal.chirp(t, f0=100, t1=time, f1=0.8*sampling_rate/2, method='logarithmic')

fig, ax = plt.subplots(1,3, figsize=(12, 3))
ax = ax.ravel()

# Figure 1
[spectrum, freqs, t, im] = ax[0].specgram(sweep, fft_size, sampling_rate, noverlap=1024-step, cmap='jet')
ax[0].set_xlabel('Time (s)')
ax[0].set_ylabel('Freq. (Hz)')
ax[0].set_title('By plt.specgram()')

# Figure 2
# plt.imshow(im.get_array(), cmap='jet', aspect='auto')
ax[1].pcolormesh(t, freqs[::-1], im.get_array(), cmap='jet', shading='gouraud')
ax[1].set_xlabel('Time (s)')
ax[1].set_ylabel('Freq. (Hz)')
ax[1].set_title('By im.get_array()')

# Figure 3
ax[2].pcolormesh(t, freqs, np.log(spectrum), cmap='jet', shading='gouraud')
ax[2].set_xlabel('Time (s)')
ax[2].set_ylabel('Freq. (Hz)')
ax[2].set_title('By np.log(spectrum)')

plt.tight_layout()
plt.savefig('output.png', dpi=300)

# Compare values
print(im.get_array()[:4, :4])
print('n')
print(np.log(spectrum)[:4, :4])
  • Comparison
[[-204.33209898 -219.84393176 -214.14525629 -203.59444001]
 [-201.31048748 -216.82498615 -210.87029865 -200.53882072]
 [-201.27672382 -216.79880707 -210.16099174 -200.40558268]
 [-201.22097657 -216.75575367 -209.18803361 -200.19223749]]
 
[[-22.84165273 -20.81034315 -21.35552053 -27.64431125]
 [-21.8274879  -20.04740739 -20.59310002 -24.67213702]
 [-21.15707948 -19.84181708 -20.38858459 -23.23503521]
 [-20.4209166  -19.50932131 -20.05665799 -22.21228083]]
  • STFT vs. CWT

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal
# - [How to plot Spectrogram using STFT in python?](https://stackoverflow.com/questions/43109708/how-to-plot-spectrogram-using-stft-in-python/48205933#48205933)

fs = 1024
fft_size = 1024
step = fft_size/16
time=2
tt = np.arange(0, time, 1/fs) # 
# sweep = signal.chirp(t, f0=100, t1=time, f1=0.8*sampling_rate/2, method='logarithmic')
sweep = signal.chirp(tt, f0=30, t1=time, f1=500, method='quadratic')
print(len(sweep))

# fig, ax = plt.subplots(1,3,figsize=(24,6))
fig = plt.figure(figsize=(12,6))

# ax[0]
ax1 = plt.subplot(211)
ax1.plot(sweep)
ax1.set_title('Raw signal')

ax2 = plt.subplot(234)
[spectrum, freqs, t, im] = ax2.specgram(sweep, Fs=fs, cmap='jet')
ax2.set_xlim(0,2)
ax2.set_title('By specgram')

ax3 = plt.subplot(235)
f1, t1, Zxx = signal.stft(sweep, fs=fs)
ax3.pcolormesh(t1, f1, np.log(np.abs(Zxx)), shading='gouraud', cmap='jet') # vmin=0,  vmax=np.max(sweep),
ax3.set_title('By STFT')

ax4 = plt.subplot(236)
import pywt
wavename = 'cmor3-3'
totalscal = 256
fc = pywt.central_frequency(wavename)
cparam = 2 * fc * totalscal
scales = cparam / np.arange(totalscal, 1, -1)
[cwtmatr, frequencies] = pywt.cwt(sweep, scales, wavename, 1.0 / fs)
ax4.pcolormesh(tt, frequencies, np.abs(cwtmatr), shading='gouraud', cmap='jet') 
ax4.set_title('By CWT')
ax4.set_xlim(0,2)

plt.tight_layout()
  • Reference
    How to plot Spectrogram using STFT in python?
转载请注明:文章转载自 www.mshxw.com
本文地址:https://www.mshxw.com/it/846331.html
我们一直用心在做
关于我们 文章归档 网站地图 联系我们

版权所有 (c)2021-2022 MSHXW.COM

ICP备案号:晋ICP备2021003244-6号