尝试探索时频分析以获取实验信号的有价值信息。
我需要获得这样一个图,它将可视化频谱随时间的变化-对于所有信号持续时间-以提取关于当前信号的一些有价值的信息。
尝试使用wavelets来实现此目的。
这是一个link with example,我用它来计算系数和可视化频谱图。
我收到的频谱图并不反映我所拥有的数据,似乎我对可视化技术或小波计算有一些误解。
第一次
1.技术问题:如何将原始信号和小波谱图显示在同一张图上,并共用x轴?(以便能够方便地浏览数据)
1.我看不到信号第一部分的共振,但我可以使用时域信号相对容易地做到这一点。为什么?
1.左上角和右上角圆角区域中的数据是什么?
1.也许我选择了错误的scales
或levels
(请参阅提供的代码中的数组,我已经玩过这些数字,但我不知道如何选择它们的值?)
1.如何在左侧刻度上显示频率-而不是句点?
1.如何添加色彩Map表以了解当前的功率系数值?
1.也许有一些工具/方法-自动获得levels
和scales
值?或者一些一般性的建议/教程来了解如何获得平滑的时频分析数据可视化?
范例重现和范例数据。
# wavelet application for time-frequency analysis
# example used https://ataspinar.com/2018/12/21/a-guide-for-using-the-wavelet-transform-in-machine-learning/
import pywt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib widget
### loading data
# initial data for analysis
tr_link = 'https://drive.google.com/uc?export=download&id=1FEQvfqQ8BE1rS8BuTqQeRsTBJx_kiH15'
# input dataframe - transmission data - experimental results with the corresponding theoretical data and uncertainty data
tr_df = pd.read_csv(tr_link, index_col=0)
#data with the actual resonance positions
ladder_df = pd.DataFrame({'E': {0: 40.29368109287083, 1: 55.462104682940094, 2: 66.19723693698184, 3: 81.22602772049886, 4: 101.05019001782576, 5: 140.79384791944875, 6: 157.22007983931803, 7: 171.43325808950274, 8: 217.16542649339505, 9: 221.41056096126604, 10: 243.58428941894317, 11: 271.0699959595112, 12: 285.04407065596706, 13: 300.3253362518625, 14: 309.7129973708237, 15: 314.68310009948743, 16: 323.26453920844216, 17: 342.54164584343937, 18: 356.5478102534998, 19: 382.23962909125567, 20: 397.8070178951269, 21: 421.2843581368074, 22: 433.36971676875646, 23: 440.66505067627634, 24: 466.7307530820385, 25: 477.03060262148057, 26: 482.96167606482015, 27: 497.424475651996, 28: 515.5927570438026, 29: 544.9552781653138, 30: 549.3628110144758, 31: 594.9779029029645, 32: 613.1079845026825, 33: 634.7001884990783, 34: 667.2508925896644, 35: 690.0964618342111, 36: 706.2573081949196, 37: 723.455953863036, 38: 753.9685262219448, 39: 764.1023564395969, 40: 783.3375513683885, 41: 803.4441302741657, 42: 814.9006657436429, 43: 829.3065932509774, 44: 858.2810554651402, 45: 883.026051087727, 46: 913.8980609316312, 47: 933.3841005436072, 48: 936.9012403488496, 49: 960.07254411068, 50: 970.1599371526584, 51: 988.6230483408726, 52: 994.1094565980452}, 'tof': {0: 0.00040406988890069487, 1: 0.0003449018934825362, 2: 0.00031598105788645324, 3: 0.00028557849361103294, 4: 0.00025638211406266864, 5: 0.00021771059568433755, 6: 0.00020620234645780873, 7: 0.00019761038050053123, 8: 0.00017594559752920423, 9: 0.00017428275970206476, 10: 0.00016631594547793704, 11: 0.0001578317819181072, 12: 0.0001539969269073221, 13: 0.00015011363900030704, 14: 0.00014787189376091746, 15: 0.00014672587278245796, 16: 0.00014480971199822796, 17: 0.00014077095253061672, 18: 0.00013804430306554194, 19: 0.00013343808731047422, 20: 0.00013086684344748386, 21: 0.0001272621057195108, 22: 0.00012552178822158457, 23: 0.00012450607469089692, 24: 0.00012107366944339282, 25: 0.00011979555386093621, 26: 0.00011907818543729406, 27: 0.00011738300517409577, 28: 0.00011535543399251944, 29: 0.00011229554167237238, 30: 0.00011185753091262304, 31: 0.00010761419597738966, 32: 0.0001060606849751197, 33: 0.00010429807317815873, 34: 0.00010180440741783211, 35: 0.00010016063185845586, 36: 9.904631835791354e-05, 37: 9.790170126465565e-05, 38: 9.596823414123544e-05, 39: 9.535185526580462e-05, 40: 9.421496658465192e-05, 41: 9.307049217268802e-05, 42: 9.243740952453541e-05, 43: 9.166004036892226e-05, 44: 9.015621997184851e-05, 45: 8.893095729395043e-05, 46: 8.747264517648235e-05, 47: 8.65896588575137e-05, 48: 8.643322569777624e-05, 49: 8.542420745803129e-05, 50: 8.499627822618716e-05, 51: 8.423006517753538e-05, 52: 8.400650521127758e-05}} )
### preprocessing steps to get data in a proper format
# rescaling the tof
tr_df['tof'] = tr_df['tof']/1e6 # to get time in seconds
# calculating the step betweeb the points byt hte TOF
tr_df['d_tof'] = tr_df['tof'] - tr_df['tof'].shift(1)
N = tr_df.shape[0] #number of points in the signal
# time step calculation - dt
min_timebin = abs(tr_df['d_tof'].min())
max_timebin = abs(tr_df['d_tof'].max())
std_timebin = abs(tr_df['d_tof'].std())
timebin_avg = abs(tr_df['d_tof'].mean())
print(min_timebin, timebin_avg, max_timebin, std_timebin)
# let's assume that it's the mean value because we have some discretization & computational errors
timebin = round(abs(tr_df['d_tof'].mean())*1e7, 0)/1e7 # average rounded step by tof between each point
t0 = tr_df['tof'].min()
dt = timebin
print('dt = ', dt)
print('t0 = ', t0)
# time array
time_ar = np.arange(0, N) * dt + t0
# signal data array
signal_ = tr_df['exp_trans'].to_numpy()
# plotting the signal
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(tr_df['tof'], tr_df['exp_trans'], label='signal')
#plotting the positions of the peaks
for index, row in ladder_df.iterrows():
ax.axvline(x = row['tof'], color = 'r', linestyle = 'dashed', linewidth = 1.0, alpha=0.5)
ax.set_xlim([time_ar[0], time_ar[-1]])
ax.set_ylabel('Signal Amplitude', fontsize=10)
ax.set_title('Signal', fontsize=10)
ax.set_xlabel('Time', fontsize=10)
ax.legend()
ax.grid()
plt.show()
# for wavelet plotting
cmap = plt.cm.seismic
# wavelet calculation
scales = np.arange(1, 64) # what are these scales??
[coefficients, frequencies] = pywt.cwt(signal_, scales, 'cmor', dt)
pwr = (abs(coefficients)) ** 2
periods = 1. / frequencies
levels = [0.015625, 0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32] # what are these levels?
contourlevels = np.log2(levels)
fig, ax = plt.subplots(figsize=(8, 5))
im = ax.contourf(time_ar, np.log2(periods), np.log2(pwr), contourlevels, extend='both', cmap=cmap)
plt.show()
P.S.实际的任务是识别噪声信号中的共振。假设真实的和虚假共振(或随机噪声)在频谱或某些参数上不同,这些参数可以帮助识别噪声数据中的共振或使用该信息来自动化识别过程。(我有许多实验信号需要评估--具有已知的真实的共振位置,因此下一步是使用ML技术自动识别实验数据中的共振)。
1条答案
按热度按时间sauutmhj1#
1.使用matplotlib的twinx创建共享图形x轴的双轴。
1.这是影响锥,描述为here。此区域(左上角和右上角)中的数据不可信,因为它们可能受到边缘效应的影响。
1.尺度定义了你如何“拉伸和压缩”小波来分析你的输入数据。你可以使用pywt的scale2freq把尺度和频率联系起来。
im = plt.contourf(time_ar, scales, np.log2(pwr), contourlevels, extend='both', cmap=cmap)
1.您指的是colorbar:简单地将
plt.colorbar()
加在plt.show()
之前。1.检查here“选择cwt秤”部分。