👨‍👩‍👧‍👦 🧖 🕚 تحليل المويجات. الجزء 2 🐥 👨🏾‍🔬 🤯

مقدمة

يناقش هذا المنشور تحليل المويجات للسلسلة الزمنية. تتوافق الفكرة الرئيسية لتحويل المويجات مع تفاصيل العديد من السلاسل الزمنية التي تُظهر التطور في وقت خصائصها الرئيسية - متوسط القيمة والتباين والفترات والسعة ومراحل المكونات التوافقية. الغالبية العظمى من العمليات التي درست في مختلف مجالات المعرفة لديها الميزات المذكورة أعلاه.

الغرض من هذا المنشور هو وصف طريقة تحويل المويجات المستمر للسلسلة الزمنية باستخدام مكتبة PyWavelets ..

قليلا من التاريخ

الجيوفيزيائي D. Morlaix في أواخر 70s من القرن العشرين. لقد واجهت مشكلة تحليل الإشارات من أجهزة استشعار الزلازل التي تحتوي على مكون عالي التردد (النشاط الزلزالي) لفترة قصيرة من الوقت ومكونات التردد المنخفض (حالة هادئة من قشرة الأرض) لفترة طويلة. يتيح لك تحويل نافذة Fourier تحليل المكون عالي التردد أو المكون ذو التردد المنخفض ، ولكن ليس كلا المكونين في وقت واحد.

لذلك ، تم اقتراح طريقة تحليل يزيد فيها عرض وظيفة النافذة للترددات المنخفضة وينخفض للترددات العالية. تم الحصول على تحول النافذة الجديد نتيجة التمدد (الانضغاط) وإزاحة الوقت لدالة توليد واحدة (ما يسمى وظيفة التحجيم - وظيفة التحجيم ، والسكالت). كانت وظيفة التوليد هذه تسمى المويجات بواسطة د. مورليكس.

وافليت د. مورليكس

from pylab import* import scaleogram as scg axes = scg.plot_wav('cmor1-1.5', figsize=(14,3)) show()

تحويل المويجات المستمر (CWT)

تحويل المويجات أحادية المستوى:

coefs ، الترددات = pywt.cwt (البيانات ، المقاييس ، المويجات)

حيث:

البيانات: (array_like) إشارة الإدخال ؛

المقاييس (array_like) : - مقياس المويجات للاستخدام. يمكنك استخدام النسبة f = scale2frequency (scale ، المويجات) / sampling_period وتحديد التردد الفعلي f. هنا ، f في hertz عند أخذ العينات في فترة بالثواني ؛

المويجات (كائن المويجات أو الاسم) : - المويجات للاستخدام ؛

sampling_period (تعويم): - فترة أخذ العينات لترددات الإخراج (معلمة اختيارية). القيم المحسوبة لـ coefs مستقلة عن أخذ العينات_فترة (بمعنى ، المقاييس لا يتم تحجيمها بواسطة أخذ عينات الفترة) ؛

coefs (array_like) : - المويجات المستمرة - تحويل إشارة الدخل لمقياس معين و المويجات ؛

الترددات (array_like) : - إذا كانت وحدة فترة أخذ العينات هي الثانية وتم ضبطها ، فسيتم عرض الترددات بالهرتز. خلاف ذلك ، يفترض فترة أخذ العينات من 1.

ملاحظة: يعتمد حجم صفائف المعامل على طول صفيف الإدخال وطول المقاييس المحددة.

الحصول على قائمة بأسماء المويجات المتوفرة المتوافقة مع cwt:

 >>> import pywt >>> wavlist = pywt.wavelist(kind='continuous') >>> wavlist

["cgau1" ، "cgau2" ، "cgau3" ، "cgau4" ، "cgau5" ، "cgau6" ، "cgau7" ، "cgau8" ، "cmor" ، "fbsp" ، "gaus1" ، "gaus2" ، " gaus3 '،' gaus4 '،' gaus5 '،' gaus6 '،' gaus7 '،' gaus8 '،' mexh '،' morl '،' shan ']

لتحليل المويجات ، يعد اختيار المويجات الأم أحد المهام الرئيسية. لذلك ، قبل كل اختيار للمويجات ، من الضروري ببساطة استخدام البرنامج التالي الذي يسمح لك بتحديد الخصائص الديناميكية للمويجات:

قائمة 1

 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()

بعد ذلك ، سننظر في وظائف المويجات وخصائصها باستخدام البرنامج أعلاه:

المكسيكي قبعة المويجات mexh:

v a r p h i l e f t (t r i g h t) = f r a c 2 s q r t 3 s q r t [4] p i c d o t e x p^{- f r a c t^{2} 2} c d o t l e f t (1 - t^{2} r i g h t)

$\ varphi \ left (t \ right) = \ frac {2} {\ sqrt {3} \ sqrt [4] {\ pi}} \ cdot exp ^ {- \ frac {t ^ {2}} {2} } \ cdot \ left (1-t ^ {2} \ right)$

المويجات "morl" Morlaix:

v a r p h i l e f t (t r i g h t) = e x p^{- f r a c t^{2} 2} c d o t c o s (5 t)

$\ varphi \ left (t \ right) = exp ^ {- \ frac {t ^ {2}} {2}} \ cdot cos (5t)$

مركب Morlet المويجات "cmorB-C"

v a r p h i l e f t (t r i g h t) = f r a c 1 s q r t p i B e x p^{- f r a c t^{2} B} c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

حيث B هو عرض النطاق الترددي ؛ C هو التردد المركزي.

فيما يلي (بدون إشارة خاصة) يتم تعيين قيم B ، C بنقطة عائمة.

موجات غوسية "gausP"

v a r p h i l e f t (t r i g h t) = C c d o t e x p^{- t^{2}}

$\ varphi \ left (t \ right) = C \ cdot exp ^ {- t ^ {2}}$

حيث: P - يتم إدخال عدد صحيح من 1 إلى 8 في اسم المويجات ، على سبيل المثال ، "gaus8" ؛ C- المدمج في التطبيع ثابت.

موجات غاوسية متكاملة "cgauP"

v a r p h i l e f t (t r i g h t) = C c d o t e x p^{- t^{2}} c d o t e x p^{- j t}

$\ varphi \ left (t \ right) = C \ cdot exp ^ {- t ^ {2}} \ cdot exp ^ {- jt}$

حيث: P - يتم إدخال عدد صحيح من 1 إلى 8 في اسم المويجات ، على سبيل المثال ، "gaus8" ويتوافق مع ترتيب مشتق وظيفة المويجات ؛ C- المدمج في التطبيع ثابت.

شانون موجات "shanB-C"

v a r p h i l e f t (t r i g h t) = s q r t B c d o t f r a c s i n (p i B t) p i B T c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ sqrt {B} \ cdot \ frac {sin (\ pi B t)} {\ pi B T} \ cdot exp ^ {j 2 pi C t}$

حيث: B هو عرض النطاق الترددي ؛ C هو التردد المركزي.

CWT في PyWavelets ينطبق على البيانات المنفصلة - تلافحات مع عينات من التكامل المويجات. إذا كان المقياس صغيرًا جدًا ، فسنحصل على مرشح منفصل مع أخذ العينات غير الكافي ، مما يؤدي إلى تجانس ، كما هو موضح في الرسم البياني لـ "cmor1.5-1.0" للمويجات.

في العمود الأيسر ، تعرض الرسوم البيانية مرشحات منفصلة تستخدم في الإلتفاف لمختلف المقاييس. العمود الأيمن هو أطياف القدرة فورييه المقابلة لكل مرشح. من خلال المقياسين 1 و 2 للرسم البياني "cmor1.5-1.0" ، يمكن ملاحظة أن التنعيم يحدث بسبب انتهاك قيود Nyquist. يمكن أن تحدث الظاهرة المشار إليها في الموجات الأخرى عند اختيار C و B ، لذلك ، عند العمل مع wavelets ، يجب عليك استخدام البرنامج - القائمة 1.

لربط مقياس معين بتردد إشارة معين ، يجب أن تكون فترة أخذ العينات إشارة معروفة. باستخدام وظيفة pywt.scale2frequency () ، يمكنك تحويل قائمة المقاييس إلى الترددات المقابلة. يعتمد الاختيار الصحيح للمقاييس على المويجات المحددة ، لذلك يجب استخدام pywt.scale2frequency () للحصول على فكرة عن نطاق التردد المقابل للإشارة.

 >>> import pywt >>> dt = 0.01 # 100 Hz sampling >>> frequencies = pywt.scale2frequency('cmor1.5-1.0', [1, 2, 3, 4]) / dt >>> frequencies

صفيف ([100. ، 50. ، 33.33333333 ، 25.])

المويجات - تحليل السلاسل الزمنية باستخدام CWT PyWavelets

مجموعة بيانات El-Nino عبارة عن مجموعة بيانات سلسلة زمنية تستخدم لتتبع El Nino وتحتوي على قياسات فصلية لدرجات حرارة سطح البحر من 1871 إلى 1997. لفهم قوة مخطط المقياس ، دعنا نتخيله لمجموعة بيانات النينو إلى جانب بيانات السلاسل الزمنية الأصلية وتحويل فورييه.

لتحليل المويجات للسلسلة الزمنية ، يجب تنفيذ الخطوات التالية:

1. حدد المويجات الأم: حدد المويجات المورليت المعقدة "cmorB-C":

v a r p h i l e f t (t r i g h t) = f r a c 1 s q r t p i B e x p^{- f r a c t^{2} B} c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

النطاق الترددي - B والتردد المركزي - الذي سنحدد من خلاله في الخطوة التالية.

2. تحديد التردد المركزي عن طريق اتخاذ dt = 0.25 لسلسلة الوقت لدينا:

 import pywt dt = 0.25 # 4 Hz sampling scale=range(1,4) frequencies = pywt.scale2frequency('cmor1.0-0.5', scale) / dt print(frequencies)

[2. 1. 0.66666667]

3. نجد تحويل فورييه للوحة الأم المويجة cmor1.0-0.5 (نستخدم القائمة 1):
سيتم تقييم المويجات المستمرة على كامل النطاق [-8.0 ، 8.0]

4. حدد البيانات للسلسلة الزمنية:

 pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|")

البيانات هي قياسات درجة حرارة سطح البحر ربع السنوية من 1871 إلى 1997. لهذه البيانات: t0 = 1871 dt = 0.25

5. نبني سلسلة زمنية مع إشارة ومتوسطها المتحرك على مخطط واحد:

ادراج 2

 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)

6. نقوم بتنفيذ تحويل فورييه وطيف الطيف من السلاسل الزمنية:

قائمة 3

 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)

7. نحدد المقاييس: المقاييس = arange (1 ، 128) ؛ المستويات = [2 ** - 4 ، 2 ** - 3 ، 2 ** - 2 ، 2 ** - 1 ، 2 ** 0 ، 2 ** 1 ، 2 ** 2 ، 2 ** 3].

8. بناء جدول الرسم البياني:

قائمة 4

 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)

يوضح الرسم البياني أن معظم قدرة الطيف تتركز على مدى فترة تتراوح من 2-8 سنوات ، وهذا يتوافق مع تردد يتراوح بين 0.125 - 0.5 هرتز. في طيف فورييه ، يركز الأزياء أيضًا على قيم التردد هذه. ومع ذلك ، فإن تحويل المويجات يعطينا أيضًا معلومات مؤقتة ، لكن تحويل فورييه لا يوفر لنا.

على سبيل المثال ، على الرسم البياني للمقياس ، نرى أنه قبل عام 1920 كان هناك العديد من التقلبات ، بينما لم يكن هناك الكثير بين عامي 1960 و 1990. يمكننا أيضا أن نرى أنه مع مرور الوقت هناك تحول من فترات أقصر إلى فترات أطول.

باستخدام وحدة المقياس

Scaleogram هو أداة ملائمة لتحليل بيانات 1D مع تحويل المويجات المستمر استنادًا إلى مكتبة PyWavelets. تم تصميم الوحدة لتوفير أداة موثوقة لتحليل البيانات بسرعة.

Scaleogram لديه الميزات التالية:

توقيع دعوة للمبتدئين بسيطة
محاور قابلة للقراءة وتكامل matplotlib النقي
العديد من الخيارات للتكبير ، مرشح الطيف ، تطبيق colorbar ، إلخ ...
دعم لوحدات التردد والتردد وفقا للعلامة
السرعة ، يستخدم خوارزمية سجل * N للتحولات
قابلية: اختبار مع python2.7 و python3.7
رسائل خطأ مفصلة والوثائق مع أمثلة

ومع ذلك ، في أمثلة الاستخدام ، لا يتم تسجيل الوصول إلى البيانات بشكل صحيح:

 nino3 = "http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat"

يظهر التحذير التالي:

nino3 = pd.read_table (nino3)
مما يؤدي إلى تحذير:
تحذير (من وحدة التحذيرات):
ملف "C: /Users/User/Desktop/wavelet/pr.py" ، السطر 7
nino3 = pd.read_table (nino3)
FutureWarning: read_table مهملة ، استخدم read_csv بدلاً من ذلك ، ويمر sep = '\ t'

يجب كتابة الوصول إلى البيانات مثل هذا:

 url="http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|")

لإظهار استخدام الرسم البياني ومقارنة النتائج مع نتائج المثال السابق ، نستخدم نفس البيانات - على القياسات الفصلية لدرجة حرارة سطح البحر من 1871 إلى 1997. لهذه البيانات: t0 = 1871 dt = 0.25

قائمة 4

 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()

إذا قارنا مخطط المقياس بالمسح الضوئي الناتج ، إذن ، باستثناء لوحة الألوان ، تكون متطابقة ، وبالتالي تظهر نفس الديناميات في السلسلة الزمنية.
القائمة 4 هي أبسط من القائمة 3 ولديها أداء أفضل.

عندما لا يسمح طيف الترددات وطيف القدرة وفقًا لـ فورييه بالحصول على معلومات إضافية حول ديناميكيات السلاسل الزمنية ، يُفضَّل استخدام scalogram.

النتائج

ويرد مثال لتحليل المويجات لسلسلة زمنية عن طريق مكتبة PyWavelets.

تحليل المويجات. الجزء 2

مقدمة

تحويل المويجات المستمر (CWT)

المويجات - تحليل السلاسل الزمنية باستخدام CWT PyWavelets

باستخدام وحدة المقياس

النتائج

More articles: