0

I have a signal from a magnetic detector that I'm interested in analyzing, I've made signal decomposition using wavedec()

coeffs = pywt.wavedec(dane_K180_40['CH1[uV]'], 'coif5', level=5)

And I've received decomposition coefficients as follows:

cA1, cD5, cD4, cD3, cD2, cD1 = coeffs

These are ndarrays objects with various lengths. cD1 is (1519,) cD2 is (774,) and so on. Different length of arrays is my main obstacle. coefficients

My question:

I have to make DWT Scaleogram and I can't stress it enough that I've tried my best and couldn't do it.

What is the best approach? Using matpllotlib's imshow() as follows:

plt.imshow(np.abs([cD5, cD4, cD3, cD2, cD1]), cmap='bone', interpolation='none', aspect='auto')

gives me an error

TypeError: Image data of dtype object cannot be converted to float

I've tried to google it since I'm not an expert in python and I've tried to change the ndarrays to float.

What is the best for plotting scaleogram, matshow, pcolormesh? ;D

Muhammad Arslan
  • 380
  • 1
  • 4
  • 10
StrEtron
  • 1
  • 1

1 Answers1

0

Basically, each cDi array has half the amount of samples as the previous array (this is not the case for every mother wavelet!), so I create a 2D numpy array where the first element is the 'full' amount of samples, and for each subsequent level I repeat the samples 2^level times so that the end result is a rectangular block. You can pick whether you want the Y-axis plotted as a linear or as a logarithmic scale.

# Create signal
xc = np.linspace(0, t_n, num=N)
xd = np.linspace(0, t_n, num=32)
sig = np.sin(2*np.pi * 64 * xc[:32]) * (1 - xd)
composite_signal3 = np.concatenate([np.zeros(32), sig[:32], np.zeros(N-32-32)])

# Use the Daubechies wavelet
w = pywt.Wavelet('db1')
# Perform Wavelet transform up to log2(N) levels
lvls = ceil(log2(N))
coeffs = pywt.wavedec(composite_signal3, w, level=lvls)
# Each level of the WT will split the frequency band in two and apply a
# WT on the highest band. The lower band then gets split into two again,
# and a WT is applied on the higher band of that split. This repeats
# 'lvls' times.
#
# Since the amount of samples in each step decreases, we need to make
# sure that we repeat the samples 2^i times where i is the level so
# that at each level, we have the same amount of transformed samples
# as in the first level. This is only necessary because of plotting.
cc = np.abs(np.array([coeffs[-1]]))
for i in range(lvls - 1):
    cc = np.concatenate(np.abs([cc, np.array([np.repeat(coeffs[lvls - 1 - i], pow(2, i + 1))])]))

plt.figure()
plt.xlabel('Time (s)')
plt.ylabel('Frequency (Hz)')
plt.title('Discrete Wavelet Transform')
# X-axis has a linear scale (time)
x = np.linspace(start=0, stop=1, num=N//2)
# Y-axis has a logarithmic scale (frequency)
y = np.logspace(start=lvls-1, stop=0, num=lvls, base=2)
X, Y = np.meshgrid(x, y)
plt.pcolormesh(X, Y, cc)

use_log_scale = False

if use_log_scale:
    plt.yscale('log')
else:
    yticks = [pow(2, i) for i in range(lvls)]
    plt.yticks(yticks)

plt.tight_layout()
plt.show()
Oguz Meteer
  • 79
  • 9