[Python] Return Dominant frequencies from S transform

Based on following link :S transform

i have implemented following spectral analyzer code (just call it analyzer)

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import irfft
from scipy.signal.windows import tukey
from scipy.signal import periodogram
from twistpy.utils import stransform,istransform
dt = 0.01 # Sampling interval (s)
tmin = 0 # Start time of the signal (s)
tmax = 2.93 # End time of the signal (s)

t = np.arange(tmin, tmax, dt) # Time axis
# print(t)
f1 = 10 # Frequency of first sinusoid
f2 = 10.6 # Frequency of second sinusoid

u = 30*np.sin(2 * np.pi * f1 * t) + 40 * np.cos(2 * np.pi * f2 * t) # Generate signal
L =10
K =len(u)-L+1

X =np.row_stack([u[i:i+L] for i in range(0,K)])
print(X.shape)
U, S, V = np.linalg.svd(X)
U1=np.concatenate((U[:,0],U[:,2]))
print(U1.shape)

[f1,pxx]=periodogram(u,fs=100)
# taper = tukey(len(t)) # Taper the signal to avoid edge artifacts
# u *= taper
# [f,pxx]=periodogram(u,fs=1/dt)
# plt.plot(f,pxx)
# plt.show()
u_stran_k1, f = stransform(u, k=11) # Compute the S-transform with k=1
f_Hz = f / dt # Convert frequency vector to Hz
f_Hz =np.array(f_Hz)
# print(f_Hz)
# print(f_Hz[np.where((f_Hz>9.9) & (f_Hz<11))] )
plt.subplot(1,3,1)
plt.plot(f1,pxx)
plt.xlabel('frequency')
plt.ylabel('Power spectral density')
# print(f_Hz[np.abs(u_stran_k1)>5])
# Plot the magnitude of the S-transform for k=1
plt.subplot(1,3,2)
plt.imshow(
np.abs(u_stran_k1),
origin="lower",
extent=[tmin, tmax, f_Hz[0], f_Hz[-1]],
aspect="auto",cmap='gray'
)
plt.subplot(1,3,3)
plt.hist( np.abs(u_stran_k1))
# print(max(np.abs(u_stran_k1).ravel()))
plt.show()


in this code dt =0.01 so fs =1/ds =100, also we have two frequency :

f1 = 10 # Frequency of first sinusoid
f2 = 10.6 # Frequency of second sinusoid


difference between them df =f2-f1=10.6-10=0.6 duration of the signal is 2.93 second, so that total number of sample is : 293 , we know that spectral resolution depends on number of samples and depend on time period of observation, in our case periodogram can resolve two signal if difference between frequencies is greater than 0.341296928327645, i would like to check if S transform works fine for this case and here is image :

my question is : from image itself or from np.abs(u_stran_k1) how can i determine frequencies? thanks in advance

Continue reading...