from pylab import*
import scaleogram as scg
axes = scg.plot_wav('cmor1-1.5', figsize=(14,3))
show()
>>> import pywt
>>> wavlist = pywt.wavelist(kind='continuous')
>>> wavlist
import numpy as np
import pywt
import matplotlib.pyplot as plt
vav='cmor1.5-1.0'
wav = pywt.ContinuousWavelet(vav)
# вывести диапазон, в котором будет оцениваться вейвлет
print("Непрерывный вейвлет будет оцениваться во всем диапазоне [{}, {}]".format( wav.lower_bound, wav.upper_bound))
width = wav.upper_bound - wav.lower_bound
scales = [1, 2, 3, 4, 10, 15]
max_len = int(np.max(scales)*width + 1)
t = np.arange(max_len)
fig, axes = plt.subplots(len(scales), 2, figsize=(12, 6))
for n, scale in enumerate(scales):
# Следующий код адаптирован из внутренних частей cwt
int_psi, x = pywt.integrate_wavelet(wav, precision=10)
step = x[1] - x[0]
j = np.floor(
np.arange(scale * width + 1) / (scale * step))
if np.max(j) >= np.size(int_psi):
j = np.delete(j, np.where((j >= np.size(int_psi)))[0])
j = j.astype(np.int)
# normalize int_psi для более простого построения
int_psi /= np.abs(int_psi).max()
# дискретные выборки интегрированного вейвлета
filt = int_psi[j][::-1]
# CWT состоит из свертки фильтра с сигналом в этой шкале
# Здесь мы строим это дискретное ядро свертки в каждом масштабе.
nt = len(filt)
t = np.linspace(-nt//2, nt//2, nt)
axes[n, 0].plot(t, filt.real, t, filt.imag)
axes[n, 0].set_xlim([-max_len//2, max_len//2])
axes[n, 0].set_ylim([-1, 1])
axes[n, 0].text(50, 0.35, 'scale = {}'.format(scale))
f = np.linspace(-np.pi, np.pi, max_len)
filt_fft = np.fft.fftshift(np.fft.fft(filt, n=max_len))
filt_fft /= np.abs(filt_fft).max()
axes[n, 1].plot(f, np.abs(filt_fft)**2)
axes[n, 1].set_xlim([-np.pi, np.pi])
axes[n, 1].set_ylim([0, 1])
axes[n, 1].set_xticks([-np.pi, 0, np.pi])
axes[n, 1].set_xticklabels([r'$-\pi$', '0', r'$\pi$'])
axes[n, 1].grid(True, axis='x')
axes[n, 1].text(np.pi/2, 0.5, 'scale = {}'.format(scale))
axes[n, 0].set_xlabel('time (samples)')
axes[n, 1].set_xlabel('frequency (radians)')
axes[0, 0].legend(['real', 'imaginary'], loc='upper left')
axes[0, 1].legend(['Power'], loc='upper left')
axes[0, 0].set_title('filter: wavelet - %s'%vav)
axes[0, 1].set_title(r'|FFT(filter)|$^2$')
plt.show()
>>> import pywt
>>> dt = 0.01 # 100 Hz sampling
>>> frequencies = pywt.scale2frequency('cmor1.5-1.0', [1, 2, 3, 4]) / dt
>>> frequencies
import pywt
dt = 0.25 # 4 Hz sampling
scale=range(1,4)
frequencies = pywt.scale2frequency('cmor1.0-0.5', scale) / dt
print(frequencies)
рядаpd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|")
from numpy import *
from scipy import *
import pandas as pd
from pylab import *
import pywt
def get_ave_values(xvalues, yvalues, n = 5):
signal_length = len(xvalues)
if signal_length % n == 0:
padding_length = 0
else:
padding_length = n - signal_length//n % n
xarr = array(xvalues)
yarr = array(yvalues)
xarr.resize(signal_length//n, n)
yarr.resize(signal_length//n, n)
xarr_reshaped = xarr.reshape((-1,n))
yarr_reshaped = yarr.reshape((-1,n))
x_ave = xarr_reshaped[:,0]
y_ave = nanmean(yarr_reshaped, axis=1)
return x_ave, y_ave
def plot_signal_plus_average(time, signal, average_over = 5):
fig, ax = subplots(figsize=(15, 3))
time_ave, signal_ave = get_ave_values(time, signal, average_over)
ax.plot(time, signal, label='Сигнал')
ax.plot(time_ave, signal_ave, label = 'Скользящее среднее сигнала (n={})'.format(5))
ax.set_xlim([time[0], time[-1]])
ax.set_ylabel('Амплитуда сигнала', fontsize=18)
ax.set_title('Сигнал + Скользящее среднее сигнала', fontsize=18)
ax.set_xlabel('Время', fontsize=18)
ax.legend()
show()
df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|")
N = df_nino.shape[0]
t0=1871
dt=0.25
time = arange(0, N) * dt + t0
signal = df_nino.values.squeeze()
scales = arange(1, 128)
plot_signal_plus_average(time, signal)
from numpy import *
from scipy import *
import pandas as pd
from pylab import *
import pywt
def get_ave_values(xvalues, yvalues, n = 5):
signal_length = len(xvalues)
if signal_length % n == 0:
padding_length = 0
else:
padding_length = n - signal_length//n % n
xarr = array(xvalues)
yarr = array(yvalues)
xarr.resize(signal_length//n, n)
yarr.resize(signal_length//n, n)
xarr_reshaped = xarr.reshape((-1,n))
yarr_reshaped = yarr.reshape((-1,n))
x_ave = xarr_reshaped[:,0]
y_ave = nanmean(yarr_reshaped, axis=1)
return x_ave, y_ave
def get_fft_values(y_values, T, N, f_s):
f_values = linspace(0.0, 1.0/(2.0*T), N//2)
fft_values_ = fft(y_values)
fft_values = 2.0/N * abs(fft_values_[0:N//2])
return f_values, fft_values
def plot_fft_plus_power(time, signal):
dt = time[1] - time[0]
N = len(signal)
fs = 1/dt
fig, ax = subplots(figsize=(15, 3))
variance = std(signal)**2
f_values, fft_values = get_fft_values(signal, dt, N, fs)
fft_power = variance * abs(fft_values) ** 2 # FFT power spectrum
ax.plot(f_values, fft_values, 'r-', label='FFT преобразование')
ax.plot(f_values, fft_power, 'k--', linewidth=1, label='Спектр мощности')
ax.set_xlabel('Частота[Герц / год]', fontsize=18)
ax.set_ylabel('Амплитуда', fontsize=18)
ax.legend()
show()
df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|")
N = df_nino.shape[0]
t0=1871
dt=0.25
time = arange(0, N) * dt + t0
signal = df_nino.values.squeeze()
scales = arange(1, 128)
plot_fft_plus_power(time, signal)
from numpy import *
import pandas as pd
from pylab import *
import pywt
def plot_wavelet(time, signal, scales,
waveletname = 'cmor1.0-0.4',
cmap = plt.cm.seismic,
title = 'Вейвлет-преобразование(Спектр мощности) сигнала',
ylabel = 'Период (год)',
xlabel = 'Время'):
dt = time[1] - time[0]
[coefficients, frequencies] = pywt.cwt(signal, scales, waveletname, dt)
power = (abs(coefficients)) ** 2
period = 1. / frequencies
levels = [2**-4 , 2**-3 , 2**-2 , 2**-1 , 2**0 , 2**1 , 2**2 , 2**3]
contourlevels = log2(levels)
fig, ax = subplots(figsize=(15, 10))
im = ax.contourf(time, log2(period), log2(power), contourlevels, extend='both',cmap=cmap)
ax.set_title(title, fontsize=20)
ax.set_ylabel(ylabel, fontsize=18)
ax.set_xlabel(xlabel, fontsize=18)
yticks = 2**arange(np.ceil(log2(period.min())), ceil(log2(period.max())))
ax.set_yticks(log2(yticks))
ax.set_yticklabels(yticks)
ax.invert_yaxis()
ylim = ax.get_ylim()
ax.set_ylim(ylim[0], -1)
cbar_ax = fig.add_axes([0.95, 0.5, 0.03, 0.25])
fig.colorbar(im, cax=cbar_ax, orientation="vertical")
show()
df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|")
N = df_nino.shape[0]
t0=1871
dt=0.25
time = arange(0, N) * dt + t0
signal = df_nino.values.squeeze()
scales = arange(1, 128)
plot_wavelet(time, signal, scales)
nino3 = "http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat"
url="http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat"
nino3 = pd.read_csv(url, sep = "|")
import pandas as pd
import pywt
from numpy import*
import scaleogram as scg
from pylab import*
url="sst_nino3.dat"
nino3 = pd.read_csv(url, sep = "|")
data = nino3.values.squeeze()
N = data.size
t0 = 1871; dt = 0.25
time = t0 + arange(len(data))*dt
wavelet = 'cmor1-0.5'
ax = scg.cws(time, data, scales=arange(1, 128), wavelet=wavelet,
figsize=(14, 3), cmap="jet", cbar=None, ylabel='Период (год)', xlabel="Время [Год]", yscale="log",
title='Вейвлет-преобразование временного ряда\n(спектр мощности)')
ticks = ax.set_yticks([2,4,8, 16,32])
ticks = ax.set_yticklabels([2,4,8, 16,32])
show()
К сожалению, не доступен сервер mySQL