Continuous Wavelet Transform with Scipy.signal (Py

2019-04-01 00:01发布

问题:

I search to draw a time-frequency signal with a discrete temporal signal (sampling step = 0.001sec). I use Python and the library Scipy.signal. I use the function cwt(data, wavelet, widths), which returns a matrix, to do a continuous wavelet transform, with the complex morlet wavelet (or gabor wavelet). Unfortunately, there is not a lot of documentations of this use. The best which I found are: - this for Matlab (I try to find the same scale-time result) but I have naturally not access to the same fonctions, - And this which explain what is continuous wavelet transform, without details of wavelet parameters.

First step: Obtain a scale-translation signal. In doubt, I associated directly the array “widths” with the array of the possible different scales. Because, I don’t understood what is parameter width if it’s not scale. Perhaps, you would tell me “it’s the width of your current wavelet”! But, even now, I'm not sure how link width with scale… In the Morlet documentation in Scipy, it seems that the link could be: "s: Scaling factor, windowed from -s*2*pi to +s*2*pi", so, I thought that width = 4*pi*scale (width=width of the window). But when I draw the wavelets, more scale increases, more the visual width of the wavelet decreases...

My second problem is to find and draw the equivalent with frequency. In literature, I find this formula: Fa = Fc / (s*delta), where Fa is the final frequency, Fc the center frequency of a wavelet in Hz, s the scale and delta the sampling period. So, ok for scale (if I find the link with the width) and delta (=0.001sec), but it’s more complicated with center frequency of the wavelet. In scipy documentation, I find that: “The fundamental frequency of this wavelet [morlet wavelet] in Hz is given by f = 2*s*w*r / M, where r is the sampling rate [s is here Scaling factor, windowed from -s*2*pi to +s*2*pi. Default is 1; w the width; and M the length of the wavelet].” I think it’s the center frequency, is it?

Thank you

Here my remanied code for cwt():

def MyCWT(data, wavelet, scales):

output = zeros([len(scales), len(data)], dtype=complex)

for ind, scale in enumerate(scales):

    window = scale*4*pi*10#Number of points to define correctly the wavelet
    waveletLength = min(window, len(data))#Number of points of the wavelet
    wavelet_data = wavelet(waveletLength, s=scale)#Need to precise w parameter???

    #To see the wavelets:
    plot(wavelet_data)
    xlabel('time (10^-3 sec)')
    ylabel('amplitude')
    title('Morlet Wavelet for scale='+str(scale)+'\nwidth='+str(window))
    show()

    #Concolution to calculate the current line for the current scale:
    z = convolve(data, wavelet_data, mode='same')

    i = 0
    for complexVal in z:
        output[ind][i] = complex(complexVal.real, complexVal.imag)         
        i+=1

return output

回答1:

The widths parameter is an array of width sizes to which the wavelet is stretched to before convolving the wavelet with the data.

You should choose a range starting with a value slightly smaller than your expected signal width, up to slightly larger. The more values you supply, the slower the calculation but the higher the resolution.

Check out the documentation or the referenced paper Bioinformatics (2006) 22 (17): 2059-2065. doi: 10.1093/bioinformatics/btl355 for more information.