الفلك

فهم الرسم البياني للدورة الشهرية اللومبي القرمزي

فهم الرسم البياني للدورة الشهرية اللومبي القرمزي



We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

أحاول استخدام الرسم البياني الدوري لمعرفة متى تكون الإشارة دورية أم لا من خلال اتباع البرنامج التعليمي الخاص بالدورة الشهرية Lomb-scargle الفلكية هنا.

http://docs.astropy.org/en/stable/stats/lombscargle.html

لقد قمت بمحاكاة بعض البيانات ، أحدها عبارة عن جيبية (فترة = 200) وواحدة تمثل انحرافًا غاوسيًا (أي حدث عابر واحد). كان الأمل هو أن ينتقي الرسم البياني الفترة الزمنية للكائن الدوري ويعطي فترة للوقت الذي قد يعني وقوع حدث واحد في النافذة.

لسوء الحظ ، النتائج غير منطقية على الإطلاق. لقد أرفقت الكود في النهاية والأرقام التي تم إنشاؤها أدناه. نظرًا للضوضاء العشوائية في المحاكاة ، تختلف كل نتيجة وأقدم مثالين للنتائج أدناه. أستخدم نفس الطريقة الموضحة في الرابط أعلاه حيث نستخدم الوظيفةLombScargle.model ()على أفضل تردد مناسب منالتردد [argmax [الطاقة]]

الخط الأحمر هو الوظيفة الحقيقية التي قمت بمحاكاة البيانات منها. اللون الأخضر هو الأنسب من الرسم البياني. مخططات اليد اليمنى هي PSD من مخطط الدورة الشهرية.

مثال 1

هنا يمكننا أن نرى أفضل تردد مناسب للجيوب الأنفية (الجزء العلوي الأيمن) الذي تم انتقاؤه هو 0.105 (أي فترة 9-10 أيام) وهو ليس قريبًا من التردد 0.05 الذي أتوقعه لشيء مدته 200 يوم ، ومع ذلك ، عندما أقوم بتغذية أفضل تردد مناسب يبلغ 0.105 لمُجرب طراز لومب سارجل ، فإن منحنى دوري مطابق بشكل جيد يخرج بفترة 200 يوم؟

غير منطقي.

مثال 2

هنا قمت بتشغيل الكود مرة أخرى وهذه المرة تم تبديل النتائج؟ إنها تتناسب مع فترة زمنية كبيرة جدًا بحيث يمكنني القول بثقة أنها حدث عابر واحد ، لكن الملاءمة الجيبية فظيعة. أفضل تردد لا يزال 0.105 (الفترة = 10) ومع ذلك ، فإن نموذج اللومب القرمزي يغلف شيئًا يبدو أنه يحتوي على فترة 60 يومًا ، وهو خطأ؟

هل يمكنني الحصول على توضيح إذا كنت أفعل شيئًا خاطئًا؟ لقد تم إخباري أن الرسم البياني الدوري هو الأداة الفعلية لبيانات مثل هذه العينات غير المنتظمة حتى الآن ... تبدو النتائج مروعة في نصف الوقت.

لتوضيح أسئلتي

  1. هل يمكن تفسير كيف في الرسم الأول أن أفضل شكل متناسب يبلغ 0.105 والذي أقوم بإدخاله في مجرب نموذج lomb-Scargle الفلكي يخلق بطريقة ما نسخة جيبية متطابقة بشكل صحيح مع تردد 0.05؟ ما هو التفسير؟

  2. لماذا توجد 5 قمم قوية في الرسم البياني الأيمن العلوي للدورة الشهرية في المثال 1 عندما أتوقع 1 فقط؟ الوسطان قريبان من القيمة الحقيقية 0.05 (عند 0.045 و 0.055)

هذا هو الرمز المختصر المستخدم لمحاكاة البيانات ورسمها

استيراد numpy كـ np import matplotlib.pyplot كـ plt استيراد scipy.stats كما sc من astropy.stats استيراد LombScargle import math #simulate parameters data_range = [i for i in range (1،1001)] number_of_samples = 50 gauss_skew = sc.skewnorm. pdf skew = -10 period = 200 location = data_range [int (len (data_range) / 2)] y1 = [(2 * (1. + np.sin (2. * np. pi * x / period)) + np .random.normal (loc = 0.0، scale = 0.5)) لـ x في نطاق البيانات] أخطاء 1 = [np.random.normal (loc = 0.0، scale = 1) لـ x في نطاق البيانات] y2 = [(1000 * (gauss_skew ( x ، الانحراف ، loc = الموقع ، المقياس = 50)) + np.random.normal (loc = 0.0 ، مقياس = 1)) لـ x في نطاق البيانات] أخطاء 2 = [np.random.normal (loc = 0.0 ، مقياس = 1 ) من أجل x في data_range] sample_rate = int (len (data_range) / number_of_samples) # لتقليل حجم البيانات قليلاً y1 = y1 [:: sample_rate] y2 = y2 [:: sample_rate] أخطاء 1 = أخطاء 1 [:: sample_rate] أخطاء 2 = أخطاء 2 [:: sample_rate] x1 = x2 = data_range [:: sample_rate] الحقيقة 1 = [2 * (1. + np.sin (2. * np. pi * x / period)) لـ x في نطاق البيانات] الحقيقة 2 = [1000 * (جاو ss_skew (x، -10، loc = location، scale = 50)) لـ x في data_range] الحقيقة = [الحقيقة 1، الحقيقة 2] x = [x1، x2] y = [y1، y2] الأخطاء = [الأخطاء 1 ، الأخطاء 2] الشكل ، ax = plt.subplots (nrows = 2، ncols = 2) ax [0] [0]. errorbar (x1، y1، yerr = errors1، fmt = 'o') ax [0] [0] .set_xlabel (' time ') ax [1] [0] .errorbar (x2، y2، yerr = errors2، fmt =' o ') ax [1] [0] .set_xlabel (' time ') ax [0] [1] .set_xlabel ("التردد") ax [1] [1] .set_xlabel ("التردد") لـ i في النطاق (0،2): التردد ، الطاقة = LombScargle (x [i] ، y [i] ، الأخطاء [i]) .autopower () # احصل على أفضل تردد مناسب كما في البرنامج التعليمي best_frequency = التردد [np.argmax (power)] print ('best frequency:'، best_frequency) t_fit = np.linspace (x [i] [0]، math .floor (x [i] [- 1])، num = 50) #Fit the best fit Frequency #plot the best model on the best fit y_fit = LombScargle (x [i]، y [i]، errors [ i]). model (t_fit، best_frequency) ax [i] [0] .plot (t_fit، y_fit، 'g') ax [i] [0] .plot (data_range، truths [i]، 'r') # ارسم فأس PSD [i] [1]. plot (التردد ، الطاقة) ax [i] [1] .axvline (x = best_frequency ، color = 'black' ، ls = "-") plt.show ()

أعتقد أنني أرى ما يجري. أفضل ما يناسبك هو في الواقع تردد ضربات بين التردد الحقيقي (0.005) ومضاعفة تردد أخذ العينات (0.05). هذا بالفعل ينتج نموذجًا يمر عبر نقاط البيانات الخاصة بك ولكن نظرًا لأنك قمت برسم 50 نقطة فقط في النموذج ، لم تتمكن من رؤيته.

إذا قمت بتغيير هذا السطر t_fit = np.linspace (x [i] [0] ، math.floor (x [i] [- 1]) ، العدد = 1000)

حتى يظهر النموذج في نقاط أكثر يمكنك أن ترى أنه في الواقع بتردد أعلى بكثير.

الحيلة هي أن تعرف الإجابة (تقريبًا) قبل أن تبدأ.

إذا قمت بتحديد نطاق التردد الخاص بك في أمر Lomb-Scargle بحيث يكون أقل من تردد Nyquist:

التردد ، الطاقة = LombScargle (x [i] ، y [i] ، الأخطاء [i]). الطاقة التلقائية (nyquist_factor = 1.5)

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

في البيانات المتباعدة بشكل غير متساو ، الأمور ليست بهذه البساطة لأن تردد نيكويست غير محدد جيدًا.


ملاحظة: هذه المؤامرة تحتوي على خطأ¶

بعد نشر الكتاب ، أشار أحد القراء إلى أن هذه المؤامرة تحتوي على خطأ مطبعي. بدلاً من تمرير البيانات المزعجة إلى روتين Lomb-Scargle ، مررنا البيانات الأساسية غير المزعجة. تسبب هذا في تقدير مبالغ فيه لقوة Lomb-Scargle.

لهذا السبب ، أضفنا مؤتمرين إضافيين إلى هذا البرنامج النصي: الأول يعيد إنتاج المؤامرة الحالية بدون خطأ مطبعي. في ذلك ، نرى أنه بالنسبة للبيانات الصاخبة ، لم يتم الكشف عن الفترة لمجرد

30 نقطة خلال عشر فترات.

في المؤامرة الإضافية الثانية ، نقوم بزيادة خط الأساس وعدد النقاط بعامل عشرة. مع هذا التكوين ، يتم الكشف عن الذروة ، والجوانب النوعية للمناقشة أعلاه صحيحة.


برنامج AFNI: 3dLombScargle

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

أهمية خاصة هو تطبيق هذه الوظيفة على
سلسلة زمنية للراحة ربما تكون قد خضعت للرقابة. النظرية وراء
تعود الرياضيات والخوارزميات الخاصة بذلك إلى مجموعات منفصلة ، بشكل أساسي
في مجال التطبيقات الفيزيائية الفلكية: Vaníček (1969 ، 1971) ،
لومب (1976) ، Scargle (1982) ، Press & amp Rybicki (1989). صرخ لهم.

هذا التنفيذ الخاص يرجع إلى Press & amp Rybicki (1989) ، بقلم
ترجمة تطبيق فورتران المنشور بشكل أساسي إلى لغة C ،
أثناء استخدام GSL لـ FFT ، بدلاً من realft () NR ، وصنع
عدة تعديلات على أساس ذلك.

تم إجراء تكييف Lomb-Scargle مع تغييرات طفيفة إلى حد ما هنا
PA Taylor (الإصدار 1.4 ، يونيو 2016).

+ الاستخدام:
أدخل سلسلة زمنية حجمية 4D (مجموعة بيانات BRIK / HEAD أو NIFTI)
بالإضافة إلى ملف اختياري 1D مكون من 0 و 1 يحدد النقاط
لفرض الرقابة (على سبيل المثال ، كل 0 يمثل نقطة / حجم للرقابة)
إذا لم يتم إدخال ملف 1D ، فسيقوم البرنامج بالتحقق من وحدات التخزين الموجودة
صفر بشكل موحد واعتبر هؤلاء خاضعين للرقابة.

الناتج هو مخطط دوري LS ، يصف المقادير الطيفية
حتى بعض "الحد الأقصى للتردد" - الحد الأقصى الافتراضي هنا هو ما
تردد نيكويست للسلسلة الزمنية * كان من الممكن أن يكون * بدونها
أي رقابة. (من المثير للاهتمام أن هذا التحليل يمكن أن يكون في الواقع
يتم تطبيقه بشكل شرعي في حالات تقدير محتوى التردد و gtNyquist.
رائع!)

سيكون طيف التردد في النطاق [df، f_N] ، حيث:
df = 1 / T ، و T هي المدة الإجمالية للسلسلة الزمنية غير الخاضعة للرقابة
f_N = 1 / dt ، و dt هو وقت أخذ العينات (أي TR)
والفاصل الزمني للترددات هو أيضًا df.
يجب أن تكون هذه النطاقات وأحجام الخطوات * مستقلة * عن الرقابة
وهي ملكية لطيفة لـ Lomb-Scargle-iness.

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ الإخراج:
1) PREFIX_time.1D: ملف 1D للنقاط الزمنية التي تم أخذ عينات منها (بوحدات
ثواني) من التحليل (وربما الخاضعة للرقابة)
مجموعة البيانات.
2) PREFIX_freq.1D: ملف 1D لنقاط عينات التردد (بالوحدات
من 1 / ثانية) من الدورة الشهرية الإخراج / الطيف
مجموعة البيانات.
3) PREFIX_amp + Orig: مجموعة بيانات حجمية تحتوي على مشتق LS
أو طيف الاتساع (يُسمى افتراضيًا "أمبير") أو أ
PREFIX_pow + طيف القدرة الأصلي (انظر "-out_pow_spec" ، المسمى "pow")
واحد لكل فوكسل.
يرجى ملاحظة أن سعة الإخراج والطاقة
الأطياف "من جانب واحد" ، لتمثيل
* نصف * سعة أو قوة تردد معين
(انظر الملاحظة التالية).

+ ملاحظة حول مسائل فورييه + بارسيفال (يرجى سامح المحرج
تنسيق):
في الصيغة المستخدمة هنا ، لسلسلة زمنية x [n] بطول N ،
قيمة مخطط الدورة الشهرية S [k] مرتبطة بقيمة السعة | X [k] |:
(1) S [k] = (| X [k] |) ** 2 ،
لكل k-th التوافقي.

ترتبط نظرية بارسيفال بالتقلبات الزمنية في السعات الطيفية ،
يذكر أن (بالنسبة إلى السلاسل الزمنية الحقيقية بمتوسط ​​صفر):
(2) sum_n = (1 / N) * sum_k <| X [k] | ** 2> ،
= (1 / N) * sum_k ،
حيث ن = 0،1. N-1 و k = 0،1. N-1 (NB: A [0] = 0 ، لمتوسط ​​الصفر
مسلسل). إن LHS هو أساسًا تباين السلاسل الزمنية
(مرات N-1). ما سبق مشتق من تحويل فورييه للرياضيات ، و
أطياف Lomb-Scargle تقريبية لـ Fourier ، لذلك ما سبق
يمكن توقع الاحتفاظ بها تقريبًا ، إذا سارت الأمور على ما يرام.

نتيجة أخرى مرتبطة بفورييه هي أنه بالنسبة لسلاسل زمنية حقيقية ومنفصلة ،
السعات الطيفية / قيم القدرة متماثلة ودورية في N.
لذلك ، | X [k] | = | X [-k] | = | X [N-k-1] | (في مجموعة صفرية القاعدة
عد)
التمييز بين الترددات المفهرسة الموجبة والسالبة
يمكن اعتبارها تدل على موجات تتحرك لليمين واليسار ، والتي
كلاهما يساهم في الطاقة الإجمالية لتردد معين.
المحصلة هي أنه يمكن للمرء أن يكتب معادلة بارسيفال على النحو التالي:
(3) sum_n = (2 / N) * sum_l <| X [l] | ** 2> ،
= (2 / N) * sum_l ،
حيث ن = 0،1. N-1 و l = 0،1. (N / 2) -1 (لاحظ عامل 2 الآن
تظهر في علاقات RHS). هذه التناظرات / الاعتبارات
هي السبب

يتم إخراج قيم التردد N / 2 هنا (نفترض
أن السلاسل الزمنية ذات القيمة الحقيقية فقط هي المدخلات) ، دون أي خسارة
معلومة.

بالإضافة إلى ذلك ، بهدف التعبير عن السعة الكلية
أو قوة تردد معين ، والتي قد يرغب العديد من الأشخاص في استخدامها
تقدير معلمات "التوصيلية الوظيفية" الطيفية مثل ALFF ،
fALFF ، RSFA ، إلخ (باستخدام ، على سبيل المثال ، 3dAmptoRSFC) ، لذلك نحن
لاحظ أن السعة * الكلية * أو القدرة لتردد معين ستكون كذلك
يكون:
أ [ك] = 2 * | س [ك] |
P [k] = 2 * S [k] = 2 * | X [k] | ** 2 = 0.5 * A [k] ** 2
بدلاً من جزء السفر الأيمن / الأيسر فقط. هذه الأنواع من
الكميات (A و P) يشار إليها أيضًا بأطياف "ثنائية الجانب". ال
يمكن بعد ذلك كتابة علاقة بارسيفال الناتجة:
(4) sum_n = (1 / (2N)) * sum_l ،
= (1 / N) * sum_l

،
حيث ن = 0،1. N-1 و l = 0،1. (N / 2] -1. بطريقة ما ، يبدو الأمر أسهل
لإخراج القيم من جانب واحد ، X و S ، بحيث يكون Parsevalian
تبدو قواعد الجمع أكثر تشابهًا.

مع وضع كل ذلك في الاعتبار ، يتم إخراج نتائج 3dLombScargle بتنسيق
يتبع. للسعات ، ما يلي تقريبا. علاقة بارسيفيليان
يجب أن تبقى بين السلسلة الزمنية "هولي" x [m] من M من النقاط و
سلسلة التردد Y [l] من L

نقاط M / 2 (حيث يقترب <| Y [l] |>
سعة فورييه <| X [l] |> كعدد النقاط الخاضعة للرقابة
إنقاص و M- & gtN):
(5) sum_m = (1 / L) * sum_l ،
حيث م = 0،1. M-1 و l = 0،1. L-1. بالنسبة إلى طيف القدرة T [l]
من L.

قيم M / 2 ، ثم:
(6) sum_m = (1 / L) * sum_l
لنفس نطاقات التلخيصات.

لذا ، يرجى مراعاة ذلك عند استخدام مخرجات هنا. 3dAmpToRSFC
لهذا الغرض عند حساب المعلمات الطيفية (من
السعة).

+ تشغيل:
-prefix PREFIX: اسم بادئة الإخراج لحجم البيانات ، ملف النقطة الزمنية 1D
وملف التردد 1D.
-inset FILE: سلسلة زمنية من وحدات التخزين ، مجموعة بيانات حجمية رباعية الأبعاد.

-censor_1D C1D: صف واحد أو عمود من 1 ثانية (احتفاظ) وأصفار (خاضع للرقابة)
يصف أي مجلدات من FILE يتم حفظها في ملف
أخذ العينات والتي تخضع للرقابة ، على التوالي. ال
يجب أن يكون طول قائمة الأرقام من
نفس طول عدد المجلدات في FILE.
إذا لم يتم إدخالها ، فسيقوم البرنامج بالبحث عن مجموعات فرعية
من جميع الأصفار وافترض أن هذه تخضع للرقابة.
-censor_str CSTR: سلسلة محددات نمط AFNI لوحدات التخزين * للاحتفاظ بها
التحليل. مثل:
'[0..4,7,10..$]'
لماذا نشير إليها على أنها "سلسلة الرقيب" عندما تكون كذلك
حقًا قائمة المجلدات التي يجب الاحتفاظ بها. حسنًا ، لقد صنعت
الشعور في ذلك الوقت. يمكن لمؤرخي المستقبل مبارزة مع
حبر عنها.

- قناع قناع: اختياري ، قناع حجم لتحليل بالإضافة إلى أي
voxel بقيم صفرية بشكل موحد عبر الزمن
تنتج طيفًا صفريًا.

-out_pow_spec: قم بالتبديل لإخراج طيف السعة للتكرار
بدلا من الدورة الشهرية. في الصيغة المستخدمة
هنا ، لسلسلة زمنية بطول N ، طيف القدرة
ترتبط القيمة S بقيمة السعة X على النحو التالي:
S = (X) ** 2. (بدون هذا الخيار ، يكون الإخراج الافتراضي هو
طيف السعة.)

-nyq_mult N2: يمكن أن تتضمن مخططات الدورة الشهرية L-S ترددات أعلى مما
عادةً ما يتم اعتباره Nyquist (تم تعريفه هنا
مثل:
f_N = 0.5 * (عدد العينات) / (إجمالي الفاصل الزمني)
بشكل افتراضي ، سيكون الحد الأقصى للتردد هو ما
f_N * كان من الممكن أن يكون إذا لم يتم فرض رقابة على النقاط
حدث. (هذا يجعل من السهل مقارنة أطياف L-S
عبر مجموعة بنفس بروتوكول الفحص ، حتى لو
هناك اختلافات طفيفة في الرقابة حسب الموضوع.)
القيم المقبولة هي & gt0. (بالنسبة لأولئك الذين يقرؤون ملف
أوراق الخوارزمية ، هذا يحدد معلمة "hifac".)
إذا لم يكن لديك سبب وجيه لتغيير هذا ،
لا تغيره!
-nifti: قم بالتبديل إلى ملف وحدة التخزين * .nii.gz
(التنسيق الافتراضي هو BRIK / HEAD).

+ مثال:
3dLombScargle -prefix LSout -inset TimeSeries.nii.gz
-mask mask.nii.gz -censor_1D censor_list.txt


____________________________________________________________________________
تم إنشاء هذه الصفحة تلقائيًا في الأربعاء 23 يونيو 15:19:40 EDT 2021


نتائج

الجاموس الأفريقي

كان LSP لمواقع Cilla مميزًا للمشي العشوائي غير الدوري (بدون ذروة) ، في حين اقترح LSP لمواقع Pepper دورة يومية (الشكل 3). عند تطبيق نمط Pepper للملاحظات المفقودة على السلاسل الزمنية لموقع Cilla ، استرجعنا ، كما هو متوقع ، إشارة دورية يومية (الشكل 3). يوضح هذا أن النمط الدوري تم إنشاؤه عن طريق الارتباط التلقائي بالطريقة التي كانت بها الملاحظات مفقودة ، وليس سلوك الحيوانات. ال ص- أكدت قيمة 0.29 من 150 محاكاة في اختبار النموذج الفارغ المطبق على مواقع Pepper أن النمط الدوري المرصود في سجل تتبع Pepper كان مصطنعًا.

غادر: رسم بياني لبيانات تتبع الجاموس سيلا وفقًا لجدول أخذ العينات الأصلي. مركز: رسم بياني لبيانات Cilla عند إعادة أخذ العينات لتقليد جدول أخذ عينات Pepper. حق: الرسم البياني لبيانات تتبع Pepper. تشير الأسهم إلى الذروة المقابلة لفترة اليوم الواحد

عند تطبيق LSP على سجل سرعة Cilla (مقابل سجل الموقع سابقًا) ، وجدنا ، كما هو متوقع ، دورية يومية كبيرة (ص-value & lt0.01 من 150 محاكاة ملف إضافي 2: الشكل C5). سيلا نشاط وهكذا كانت دورية على الرغم من لها استخدام الفضاء لم يكن (الشكل 3). يوضح هذا كيف يتم فصل الأنماط الدورية لاستخدام الفضاء عن دورية في مستويات النشاط في هذا النوع.

ولوح طيور القطرس

أظهر جميع أفراد طيور القطرس الأربعة تواترًا كبيرًا (جميع ص- القيم & lt0.005 من 1000 محاكاة لكل منهما) مع فترة تقديرية تتراوح من 22 إلى 35 يومًا. في الفرد الذي قام برحلات أطول ، ووصل إلى المياه على بعد 1800 كم من مستعمرة تكاثره مع دورية تقديرية بـ 29 يومًا ، كان الارتباط بمرور القمر واضحًا جدًا (الشكل 4). خلال ثلاث دورات خاضعة للمراقبة ، غادر الطائر المستعمرة قبل وقت قصير من القمر الجديد وبدأ رحلة العودة بعد وقت قصير من اكتمال القمر ، على الأرجح حتى يتمكن هو ورفيقه من الغذاء خلال فترات مع بعض ضوء القمر.

تتبع البيانات من طائر قطرس لوح لوح. غادر: تآمر بمقياس لون يشير إلى مرحلة القمر. الألوان الزرقاء قريبة من القمر الجديد والألوان الحمراء للقمر الكامل. حق: الرسوم البيانية. المخطط الدوري الأسود مأخوذ من بيانات التتبع ، أما المخطط الدوري الأحمر فهو من جدول أخذ العينات. يشير وجود ذروة في الرسم البياني الأسود لفترة دورة قمرية واحدة (خط عمودي) ، وغيابها في الرسم البياني الأحمر ، إلى أن النمط الدوري لا ينتج عن جدول أخذ العينات

الذئاب

كانت دورية ليوم واحد ذات دلالة إحصائية لجميع الأفراد الثمانية (جميع ص- القيم & lt0.05 من 1000 محاكاة لكل شكل 5 ، اللوحة اليمنى) ، حتى لو لم يكن النمط الدوري لاستخدام الفضاء واضحًا بشكل مباشر في بيانات التتبع الأولية (الشكل 5 ، اللوحة اليسرى). بعبارة أخرى ، أظهرت تلك الذئاب الثمانية ميلًا كبيرًا للقيام بدوريات في نطاقها الأصلي على طول الطرق المتكررة يوميًا ، لكن هذا النمط تم حجبه بسبب ضوضاء الخلفية المهمة (المكونات العشوائية غير الدورية للمسارات).

تتبع البيانات من الذئاب المأهولة. غادر: رسم لفرد واحد (أماديو) حسب توقيت اليوم ، لإظهار صعوبة الكشف البصري عن النمط الدوري بسبب المكون العشوائي المهم في عملية الحركة. حق: مخطط دوري متعدد الأفراد من أماديو و 7 ذئاب أخرى ، يُظهر دورية ليوم واحد وسلسلة التوافقيات المرتبطة بها (الأسهم الحمراء)


قابلية رصد المكونات الطيفية التي تتجاوز حد نيكويست في إشارات العينات غير المنتظمة

يعد تحديد مكون الإشارة بتردد يتجاوز حد نيكويست مشكلة صعبة في نظرية الإشارة وكذلك في بعض مجالات التطبيقات المحددة مثل علم الفلك والعلوم الحيوية. نتيجة لنظرية أخذ العينات المعروفة لإشارة عينات موحدة هي أن المكون الطيفي فوق حد نيكويست يتم تسميته إلى نطاق تردد أقل ، مما يجعل التمييز بين المكونات الحقيقية والمستعارة أمرًا مستحيلًا. ومع ذلك ، فإن أخذ العينات غير المنتظم يوفر إمكانية لتقليل المكونات ذات الأسماء المستعارة وكشف المعلومات فوق حد Nyquist. في هذا البحث ، نقدم تحليلًا نظريًا لتخفيض المكونات ذات الأسماء المستعارة في المخطط الدوري اللامعلمي لخطتي أخذ عينات: نمط أخذ العينات العشوائي ونمط أخذ العينات الناتج عن تعديل تردد النبض المتكامل (IPFM) ، والأخير مقبول على نطاق واسع باعتباره معدل ضربات القلب. نموذج التوقيت. تم اقتراح صيغة عامة تتعلق بتباين انحرافات التوقيت عن مخطط منتظم مع قمع مكون مستعار. يتم توضيح العلاقة المشتقة من خلال الرسوم البيانية Lomb-Scargle المطبقة على بيانات محاكاة. البيانات التجريبية المقدمة التي تتكون من إشارة التنفس المستمدة من مخطط القلب الكهربائي وإشارة معدل ضربات القلب تدعم أيضًا إمكانية اكتشاف الترددات فوق حد نيكويست في الحالة المعروفة باسم التعرج القلبي.

1 المقدمة

يعد التعرج ظاهرة معروفة في نظرية الإشارة تحدث إذا تم أخذ عينات من إشارة الوقت المستمر بتردد

أقل من ضعف الحد الأقصى للتردد

من إشارة العينة ، أي

. انخفض تردد أخذ العينات إلى النصف ،

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

، ما لم يتم تطبيق مرشح مضاد للتشويش. ومع ذلك ، لا تنطبق هذه القواعد بشكل كامل على البيانات المأخوذة من عينات غير موحدة.

يمكن استغلال أخذ العينات غير المنتظم عن قصد في تصميم النظام ، أو يمكن أن يعكس شروط تقييدية محددة في قياس البيانات واستكشافها في بعض المجالات العلمية ، مثل علم الفلك وعلم الزلازل وعلم الوراثة والقياسات الفيزيائية الحيوية وعلم وظائف القلب. لا يمكن الحصول على البيانات في لحظات زمنية محددة ، فهي مفقودة أو تالفة. في حالات أخرى ، تملي الطبيعة أوقات أخذ العينات غير المنتظمة بشكل أساسي ، كما هو الحال في تحليل تقلب معدل ضربات القلب. يجعل أخذ العينات غير المنتظم تحليل الإشارة أكثر تعقيدًا ، لأن الغالبية العظمى من تحليل البيانات وإجراءات معالجة الإشارات الرقمية قد تم تطويرها للإشارات التي تم أخذ عينات منها بشكل موحد. من ناحية أخرى ، تعد الإشارات التي تم أخذ عينات منها غير منتظمة أكثر إفادة ويمكن اكتشاف المكونات التي تتجاوز حد نيكويست.

يحتاج مفهوم تردد نيكويست إلى بعض التوضيح عند النظر في الإشارات غير المنتظمة. لا يوجد تعريف مقبول جيدًا إذا كان تكرار أخذ العينات غير ثابت [1]. يعرف بعض المؤلفين تردد Nyquist في الحالة غير الموحدة على أنها

هو أصغر فاصل زمني بين البيانات ، بينما يستنتج البعض الآخر تردد نيكويست من متوسط ​​أوقات أخذ العينات

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

التحليل الطيفي للبيانات المأخوذة من العينات غير المنتظمة له تاريخ قديم نسبيًا بسبب الأهمية النظرية والعملية لهذه المشكلة. تم اقتراح طرق مختلفة للتعامل مع العينات غير الموحدة. يمكن الاطلاع على مراجعة شاملة لأحدث ما توصلت إليه التقنية في [1]. سنقتصر نظرنا على الأساليب المستخدمة بشكل شائع في البحث عن بيانات معدل ضربات القلب المستخدمة في هذه الورقة كبيانات توضيحية لتوضيح الظواهر التي تم تحليلها. إن أبسط طريقة للتعامل مع البيانات التي تم أخذ عينات منها بشكل غير منتظم هو تجاهل عدم انتظام أخذ العينات ، ونقل العينات إلى شبكة منتظمة ، ومعاملتها كإشارة تم أخذ عينات منها بشكل موحد. ينتج عن هذا النهج تقليل السعات الطيفية ويتم ترجمة عدم انتظام أخذ العينات إلى ضوضاء طيفية. لا يمكن تمييز المكونات الموجودة فوق تردد Nyquist عن الترددات ذات الأسماء المستعارة. تعمل الطرق القائمة على الاستيفاء على إضعاف الترددات العالية مما يجعلها غير مناسبة لتحليل المكونات فوق تردد Nyquist. تُعرف الطريقة المستخدمة على نطاق واسع والتي لا تفترض إعادة تشكيل الإشارة باسم مخطط Lomb-Scargle ، وهو في الأساس المربعات الصغرى الملائمة للبيانات المقاسة بواسطة الجيوب الأنفية [2-6].

مخطط لومب-سكارغل (LS) [2] هو طريقة مصممة خصيصًا للعينات التي يتم الحصول عليها بشكل غير منتظم. تم اقتراح هذه الطريقة في الأصل لتحليل السلاسل الزمنية الفلكية ووجدت أيضًا تطبيقًا ناجحًا في تحليل بيانات القلب والأوعية الدموية ، مثل تلك الموجودة في تحليلات تقلب معدل ضربات القلب [7 ، 8] ، والإيقاعات اليومية [9] ، والتسلسل الجيني تحليل. تكون نتائج LP في العديد من المواقف متطابقة تقريبًا مع مخطط الدورة الشهرية الكلاسيكي (شوستر) المعمم على العينات غير المنتظمة [10] وبالتالي يمكن استنتاجها من FT للبيانات المأخوذة من العينات غير المنتظمة. يتم استخدام هذا النهج في [6] حيث يتم تقديم العديد من خصائص FT التي تم أخذ عينات منها بشكل غير منتظم ، على الرغم من الإشارة إليها بشكل غير دقيق على أنها خصائص الكثافة الطيفية لقدرة لومب.

على عكس ورقتنا البحثية ، تم إجراء التحليلات [6] دون إيلاء اهتمام خاص لتقدير تأثير التعرّف ، أي القدرة على الكشف عن مكونات أعلى من حد نيكويست. تدرس هذه الورقة فائدة أخذ العينات غير المنتظم لكل من تقليل المكونات ذات الأسماء المستعارة وكشف المكونات فوق معدل Nyquist. سنقدم تحليلًا نظريًا لحالة النموذج ، عندما يتم أخذ عينات من إشارة جيبية وفقًا لنمط أخذ عينات غير منتظم محدد مسبقًا. تم اقتراح صيغة عامة تتعلق بتباين انحرافات التوقيت عن مخطط منتظم مع قمع مكون مستعار. يتم توضيح الملاحظة المثبتة نظريًا للمكونات التي تتجاوز حد Nyquist من خلال أمثلة المحاكاة وبيانات القلب والجهاز التنفسي التي تم الحصول عليها من تجربة تم ترتيبها خصيصًا لهذا الغرض.

الخطوط العريضة للورقة على النحو التالي. يلخص القسم 2 تعريفات الرسوم البيانية المستخدمة في التحليل الطيفي للبيانات غير المنتظمة. تقدم الأقسام اللاحقة التحليل الطيفي النظري للعينة الجيبية المأخوذة من مخططين لأخذ العينات. في القسم 3 من الورقة ، سنشتق معادلات الكثافة الطيفية للقدرة (PSDs) للإشارات الجيبية غير المنتظمة التي تم أخذ عينات منها بشرط أن تخضع أوقات أخذ العينات لنموذج عشوائي. يوفر القسم 4 علاجًا للطيف الجيبي عندما يتم إنشاء نمط أخذ العينات بواسطة تعديل تردد النبضة المتكامل (IPFM). يتم اشتقاق الصيغة الكمية التي تصف تقليل المكون ذي الاسم المستعار في القسم 5. يتم تقديم التحليلات التوضيحية للبيانات المحاكاة في القسم 6. التحليلات الطيفية للبيانات التجريبية - التنفس المشتق من مخطط القلب الكهربائي (ECG-) وإشارات تغير معدل ضربات القلب (HRV) - توضيح أداء طريقة LS في حالة ما يسمى بالتعريف القلبي ، عندما يتم استنباط مكونات أعلى من تردد Nyquist بسبب عدم انتظام ضربات القلب في الجهاز التنفسي.

2. مخططات دورية لإشارات العينات غير المنتظمة

يمكن تعميم مخطط الدورة الشهرية التقليدي المستند إلى الحجم التربيعي لتحويل فورييه أيضًا للإشارات غير المنتظمة التي تم أخذ عينات منها بالشكل [4]: ​​& # 13

هي عينات بيانات مأخوذة في الوقت الحالي

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

. صاغ لومب [2] مخطط الدورة الشهرية عن طريق ملاءمة المربعات الصغرى للبيانات ذات القيمة الحقيقية المرصودة لنموذج النموذج: & # 13

& # 13 متغير تأخير الوقت الزائد

& # 13 في النموذج (2) من أجل تعاميد شرط الجيب وجيب التمام ، والتي تنشأ في صياغة المعادلات العادية لمسألة التركيب. ثم يتم حساب LS PSD كـ [4] & # 13

& # 13 بشرط أن تكون الحدود التربيعية لجيب الجيب وجيب التمام في (4) متساوية ومتقاربة

، نهج LS (4) للدورة الكلاسيكية (1). بشكل عام ، (1) و (4) ليسا متكافئين. يُفضل LS بسبب الخصائص الإحصائية المفيدة لاختبار الأهمية للقمم الطيفية [4 ، 11]. من ناحية أخرى ، (1) أكثر قابلية للتتبع من الناحية الحسابية.

3. نمط أخذ العينات العشوائي

ستتخذ الإشارة التي تم تحليلها شكل الجيب المعقد مع الاتساع

يتم أخذ عينات منها في اللحظات الزمنية

& # 13 من المفترض أن تكون لحظات أخذ العينات هي الأرقام العشوائية التي تتبع النموذج التالي: & # 13

هي فترة أخذ العينات المتوسطة ،

عبارة عن سلسلة من المتغيرات العشوائية المستقلة الموزعة بشكل متماثل. يمكن استنتاج متوسط ​​سلوك (1) ، صالح تقريبًا أيضًا للطيف LS ، من خلال تطبيق التوقع على تحويل فورييه لـ (5): & # 13

& # 13 بعد الاستبدال (5) ، (6) في (7) وإعادة الترتيب ، نحصل على & # 13

& # 13 في (8) ، يمكننا التعرف على FT لعينة جيبية بشكل موحد يمكن التعبير عنها بواسطة نواة Dirichlet الدورية: & # 13

& # 13 يمكن إعادة كتابة متوسط ​​FT (8) من حيث الوظيفة المميزة التي تم إزاحتها بواسطة التردد

& # 13 حيث الوظيفة المميزة

من عدم انتظام التوقيت

, تتميز بدالة كثافة الاحتمال

& # 13 لذلك ، تخفف الوظيفة المميزة الذروات الدورية - الأسماء المستعارة أو الصور الطيفية في (11) ، بينما تحتفظ بالقمة الحقيقية عند التردد

لا يفسر التوقع الناتج لـ FT (11) بدقة التأثير العشوائي لعدم انتظام أخذ العينات. لذلك ، فإن توقع حجم مربع FT ، أي PSD بتعريفه ، يتم تقييمه ليكون أكثر اكتمالًا على الرغم من طيف إشارة عينات غير منتظمة: & # 13

& # 13 حد الجمع في (13) بعد التعويض (5) هو & # 13

، لدينا بسبب الاستقلال المفترض: & # 13

& # 13 حيث افترضنا دالة احتمالية متناظرة في (12) وبالتالي وظيفة مميزة ذات قيمة حقيقية.

، التوقع (15) هو 1 ويمكن توسيعه رسميًا كـ & # 13

في (14) ومع مراعاة (15) ، (16) ، يبسط الجمع (14) إلى & # 13

& # 13 يمكن التعرف على مصطلح التجميع في (17) على أنه FT لعينات جيبية بشكل موحد مضروبة في النافذة المثلثية. من أجل تسهيل قراءة الصيغة ، سيتم الإشارة إلى FT كـ & # 13

& # 13 تصبح الصيغة النهائية بعد ذلك & # 13

& # 13 بخلاف (11) ، يتضمن تعبير PSD (19) ، إلى جانب القمم ذات الصلة بالجيوب الأنفية ، جزءًا من الطيف السلس المتعلق بعشوائية لحظات التوقيت. حيث

، يُظهر هذا الجزء الأملس واديًا بالقرب من التردد الحقيقي للجيوب الأنفية المعقدة. تم توضيح تكوين PSD من المصطلحات الواردة في (19) في الشكل 1.


4. نمط أخذ العينات المولدة من النموذج IPFM

نموذج تعديل تردد النبض المتكامل (IPFM) مقبول على نطاق واسع كنموذج توقيت العقدة الجيبية الأذينية في نمذجة تقلب معدل ضربات القلب. إنه نموذج للتكامل والنار ، يستخدم أيضًا لوصف نشاط الخلايا العصبية. يفترض نموذج IPFM أن إشارة التعديل تتكامل وأن نبضة الزناد تتولد عندما تصل الإشارة المتكاملة إلى عتبة (الشكل 2). أوقات حدوث الإيقاع (تمثل أوقات أخذ العينات)

التي تم إنشاؤها بواسطة نموذج IPFM يمكن تعريفها على أنها [12] & # 13

هو متوسط ​​فترة القلب ، و

هي إشارة تعديل بلا أبعاد يتم تفسيرها على أنها التغيير الجزئي لمعدل ضربات القلب اللحظي بالنسبة إلى


يمكن كتابة نمط أخذ العينات كسلسلة من نبضات ديراك: & # 13

& # 13 يُشار إلى تحويل فورييه لـ (21) باسم طيف الأعداد (SPC) [12]: & # 13

& # 13 طيف إشارة العينة

يرتبط بالطيف الحقيقي

من خلال التفاف مجال التردد: & # 13

& # 13 طيف التهم تم تحليله على نطاق واسع في العديد من الأعمال [12-14]. سنلخص هنا الحقائق المطلوبة لتقييم تأثير التعرج الطيفي. يمكن إعادة كتابة SPC المحددة بواسطة (22) من حيث إشارة التعديل كـ & # 13

& # 13 ومن المثير للاهتمام ، طيف إشارة التعديل

مدرج مباشرة في (24). بجانب نبضة ديراك الموجودة في أصل التردد ، هناك عدد لا حصر له من مصطلحات تشكيل التردد (FM) الملتف مع

حاضرون. تقع الموجات الحاملة لهذه المصطلحات المشكّلة FM التي تنتج قممًا في (24) عند مضاعفات متوسط ​​تردد أخذ العينات ، والتي تشرح أصل المكونات ذات الأسماء المستعارة في الطيف (23). بالنسبة لإشارة التعديل الجيبية: & # 13

هو تردد التشكيل ، يمكن توسيع الطيف (24) باستخدام وظائف Bessel من النوع الأول [15]. السعات الطيفية للحوامل تتبع دالة بيسيل من الدرجة 0

, the argument of which is given by the gradually increased index of the frequency modulation:

is the mean sampling frequency corresponding to the first carrier

5. Aliasing Reduction Index

For quantification of the aliased components amplitudes, we will restrict our considerations to the situation when the input frequency is in the range

. Unlike the derivation presented in Section 3, a real-valued sinusoid, comprising of positive and negative frequencies, will be considered in this treatment. It is well known that the spectrum of such a uniformly sampled signal is mirrored around the Nyquist frequency. If the frequency of the sampled real-valued signal exceeds the Nyquist limit, this frequency is folded over the Nyquist frequency and appears at a wrong position, in the range

. Specifically, if an input signal has frequency

, the aliased component appears at the frequency

. If the spectrum is plotted in a range above

, have equal amplitudes and there is no possibility to distinguish them. For nonuniformly sampled signal, this statement does not hold, and differences can be observed. A component with higher amplitude is considered to be the true component and other component at mirrored frequency as its alias. We introduce a quantity to describe the observed difference between the true and aliased components and name this quantity as the aliasing reduction index (ARI):

is the spectral amplitude at the true (input) frequency

is the spectral amplitude at the corresponding aliased frequency

. Next, we will be dealing with an evaluation of the proposed ARI for both considered sampling models—the random sampling pattern and the IPFM-generated sampling pattern.

The reduction of an aliased component in the situation with the random sampling according to (6) is deduced directly from (11). For a true component, the argument in (11) is

, while the component aliased into

is originated from the negative input frequency

. The ARI can be then written in terms of the characteristic function as

The result (28) depends on a particular probability distribution function which describes the random sampling jitter. For small-to-moderate values of ARI, a more general version of (28) can be derived using power series expansion of the characteristic function. Taking into account the presumed symmetry of the probability density function and neglecting the high-order terms in the expansion, the approximated characteristic function is expressed by the variance of the timing deviations:

in (28), we obtain an approximate formula for the ARI in term of the relative timing standard deviation

which is independent of a particular probability distribution type.

The proposed index can be evaluated also for the IPFM-generated sampling pattern. The spectrum (23) of a real-valued sinusoid, nonuniformly sampled in IPFM-generated time instants, is composed from SPC shifted in the frequency direction by

. A dominant aliased term is associated with the 1st FM carrier, that appears, after

. Its amplitude is determined by the 0th-order Bessel function of the first kind for

. The resulting formula for the ARI estimation, analogical to (28), becomes

ARI can be approximated by

In order to facilitate a comparison of the obtained formula (33) with the ARI derived for the random sampling model, we will rewrite (33) by using the standard deviation of the timing instants. Integration of the IPFM model formula (20) gives

The deviations of the sampling times from a regular scheme is thus sinusoidal, with amplitude

, and the standard deviation estimated as the RMS value of a continuous-time sinusoid is

Combining (33) with (35) and (26) results in the formula:

Surprisingly, (36) has exactly the same form as (30), despite the essentially different sampling pattern.

6. Simulation Example

In order to illustrate a mechanism of the aliasing in nonuniformly sampled signal and show an extent to which an aliased component reduction can be achieved, we have performed a spectral analysis of sequences synthesized according to both described sampling models. Since this work is inspired by analysis of the cardiac aliasing effect, we have arranged the signal parameters—duration, frequency—as to fit the experimental data presented in subsequent section. Each of the simulated signals consists of 120 samples of real-valued, unit-amplitude sinusoid. The frequencies are expressed as the fractions of the mean sampling rate

. The mean sampling rate is considered as the reciprocal value of the mean sampling interval. An input frequency exceeding the Nyquist limit was set to 0.65

which is aliased to the frequency 0.35

. Conversely, an input frequency 0.35

below the Nyquist limit 0.5

and an informationless peak in the spectrum can be observed when a periodogram is plotted above the Nyquist limit. Both spectral peaks appear as identical when the sampling is uniform. The nonuniform sampling, however, introduces differences between the true component and the unwanted—aliased or imaged components. This effect can be observed in the spectra plotted in the range

instead of a conventionally used range

Figures 3 and 4 show LS periodograms for both scenarios (below and above Nyquist limit) for the sampling times randomly deviated from regular positions. Spectral amplitudes are plotted as the square roots of the LS periodograms and scaled to directly reflect amplitudes of the generated sinusoids. The sampling irregularity follows the Gaussian distribution with the variance set according to a required irregularity level. Three levels were used in simulation: low, medium, high. The low irregularity is represented by sampling jitter with the standard deviation 5% of the mean sampling time

. Such a degree of the sampling nonuniformity caused only a minor spectral amplitude difference (about 5%). The medium level is represented by two-fold reduction of alias power (3 dB), that is, ARI

0.707. The formula derived in Section 5 allowed to estimate

required in simulation as

. A maximum deviation of a sampling instant is theoretically unlimited therefore, we incorporate a practical limit as the irregularity level which ensures a monotonic increase of the sampling times. This practical maximum is used as the high level of the irregularity. Since Gaussian variables are not confident to finite interval, determination of the high level deviations must be made in probabilistic sense. The choice

ensures that randomly chosen proximate sampling times

, fulfill the monotonicity condition

with probability 0.99. Presented figures show a gradual destruction of the aliased component (Figure 4) and spectral images (Figure 3) as the sampling irregularity level increases. In the high-level irregularity, the aliasing is effectively destroyed—the aliased peak is buried in a noise caused by sampling irregularity. Notice that a noise-like spectrum increase is evident as the sampling irregularity increases. A valley near the true frequency in noise PSD, which may be expected by examination of (19), is not evident, because it is overlapped by the smooth noise spectrum coming from the negative frequency.


(a)
(b)
(c)
(a)
(b)
(c) Spectral amplitudes of the LS spectra—nonuniformly sampled signal with the random sampling pattern. The input frequency 0.35

(a)
(b)
(c)
(a)
(b)
(c) Spectral amplitudes of the LS spectra—nonuniformly sampled signal with the random sampling pattern. The input frequency 0.65

A similar set of simulations (Figures 5 and 6) was performed for the sinusoid sampled by the pattern generated by the IPFM model driven by a sinusoidal modulating signal. The sampling time instants were obtained by means of numerical solution of the IPFM equation (34) using MATLAB function fzero. The frequency of the sampled sinusoid was set equal to the frequency of the modulating signal, a condition occurred in heart rate modulation by the respiration signal. The low and medium irregularity levels were established in the same way as in the random sampling pattern. The maximum irregularity level is considered when the normalized modulation signal amplitude

reaches the value1. Due to this restriction, the achieved attenuation cannot be as high as in the random sampling pattern scenarios. Notice that numbers of spurious peaks can be found in the spectrum which are explained by a complex structure of the SPC (24).


(a)
(b)
(c)
(a)
(b)
(c) Spectral amplitudes of periodograms—nonuniformly sampled signal with IPFM generated sampling pattern. Input frequency 0.35

(a)
(b)
(c)
(a)
(b)
(c) Spectral amplitudes of periodograms—nonuniformly sampled signal with IPFM generated sampling pattern. Input frequency

7. Cardiorespiratory Data Analysis

Beat-to-beat variations of consecutive heart periods or the instantaneous heart rate, known as the heart rate variability (HRV), become intensively studied over the last decades. Analysis of the HRV provides a noninvasive insight into the cardiovascular neural control and a spectral analysis of HRV signals is a widely used method for assessment of the autonomous nervous system. Generally, the power spectra of HRV can be divided into the high-frequency (HF, 0.15–0.4 Hz), low-frequency (LF, 0.04–0.15 Hz), and very-low-frequency range (VLF, 0.003–0.04 Hz). The LF power is modulated by both sympathetic and parasympathetic controls, while the HF power mainly reflects the parasympathetic influence linked to the respiration. Therefore, a concurrent measurement of the respiration signal is helpful for HRV interpretation. To obtain the respiration activity signal, an additional specialized instrument is not necessary because the respiration can be estimated from modulation of an electrocardiogram amplitude. Such a procedure is referred to as ECG-derived respiration [16–18]. Notice that a raw ECG signal itself includes spectral components related to heart rate oscillations and respiration [19], but they are mixed in the spectrum in a complex way [20], and thus, appropriate signals must be derived from measured ECG. Both the HRV as well as the respiration signal derivation incorporate heart beat detection, and the beat occurrence times define sampling instants, irregularly spaced by nature. Although a heart rate (HR) is essentially a discrete-time signal, computed from beat-to-beat occurrence times, it can be considered as a hypothetically continuous signal representing the autonomous nervous system activity which continually modulates the heart rate. Unlike an artificially designed system, a physiological system does not include antialiasing filter. For example, 0.4 Hz frequency considered as the highest limit of the standardized HF range requires sampling at a mean heart rate 0.8 Hz, and thus, the mean heart rate of at least 48 bpm is required for reliable spectral analysis. Therefore, a bradycardia or an extended HF range caused by the elevated respiration rate can lead to the aliasing. For this kind of “physiological aliasing,” the term cardiac aliasing has been naturalized. High frequency components modulating the heart rate have been actually observed in neonates, transplanted heart patients, and animals [21, 22]. Improper signal processing or ignoring rarely occurred conditions can then result in incorrect interpretation of a spectral analysis.

To test an ability of different methods to detect frequencies beyond the Nyquist rate in a real-world measured data, we have performed an experiment exploiting respiratory sinus arrhythmia. The measured subject was asked to breathe at an elevated rate in synchrony with a moving pattern displayed on the computer screen at the frequency of 0.8 Hz. An electrocardiogram (ECG) was recorded and two signals were extracted: the respiration and the instantaneous heart rate (HR). The respiration signal was obtained as the amplitude in bandpass-filtered ECG curve taken at maxima of QRS complexes [23]. The HR was computed from beat-to-beat time distances of the QRS maxima. Both signals are thus sampled at a variable rate dictated by the heart rate. The mean heart rate of 74 bpm (1.23 Hz) was sufficiently low to observe an aliasing.

In Figure 7(a), the respiration signal is analyzed as a sequence of values equidistantly spaced at multiples of the mean heart period. A meaningful spectral range is thus up to 0.615 Hz provided that the mean heart rate was 1.23 Hz. The respiration frequency 0.8 Hz is aliased near 0.4 Hz and spreads over a wider spectral range due to wandering of the mean heart rate. The respiration sequence was then interpolated to 4 Hz sampling frequency by means of the cubic spline interpolation. A small peak at the true respiratory frequency can be detected but the aliased component near 0.4 Hz largely exceeds the true peak in its amplitude—Figure 7(b). The Lomb-Scargle periodogram manifests a clearly visible peak at 0.8 Hz, Figure 7(c), albeit the components at aliased frequency region do not disappear completely.


(a)
(b)
(c)
(a)
(b)
(c) Spectral analysis of the ECG-derived respiration signal: (a)equidistantly repositioned QRS maxima samples with frequency axis scaled according to the mean heart period, (b) spline interpolated unevenly sampled QRS maxima, (c) the Lomb-Scargle periodogram of unevenly sampled QRS maxima.

The results of application of different methods for HRV spectrum computation are summarized in Figure 8. Repositioning of unevenly spaced samples, similarly as it was observed in Figure 7(a), yields a symmetrical spectrum that is unable to detect components above the mean heart rate halved. The cubic spline interpolation, although frequently used in HRV analysis, noticeably attenuates high-frequency components above the mean Nyquist frequency. The simplest interpolation type, nearest neighbor method, gives surprisingly better outcome, comparable to the Berger resampling method [24]. As in the case of the respiration signal analysis, only the last two methods presented, the Lomb-Scargle periodogram and the SPC, are capable to show a higher amplitude at the true position 0.8 Hz than near the aliased frequency (about 0.4 Hz). The apparent increase of the amplitudes at higher frequencies in Figure 8(f) can be explained by low-pass filtering effect of the integrator in the IPFM model that does not affect the spectrum of counts. In the low frequency range, all methods give similar results.


(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f) HRV spectra obtained by means of different methods: (a)equidistantly repositioned HR samples with frequency axis scaled according to the mean heart period unevenly sampled heart rate sequence resampled to evenly spaced signal- (b) cubic-spline interpolation, (c) nearest neighbour interpolation, (d) Berger method spectra computed without resampling-, (e) Lomb-Scargle periodogram, (f) spectrum of counts.

8. Conclusion

The fact that a spectrum of nonuniformly sampled signal conveys useful information above the Nyquist limit has been pointed out in several works [25, 26]. In this work, we have presented a theoretical treatment which quantitatively explains this phenomenon. The analyses of simulated and experimental data illustrate a possibility to identify components beyond the Nyquist limit in real-world problems.

Although the spectrum of a nonuinformly sampled signal is not alias free, unlike a uniform sampling situation, the aliased spectral amplitudes are attenuated when compared to the true components amplitudes. Analogically, components below the Nyquist limit can be observed also as images above this limit. The analytic treatment of these unwanted components constitutes a key contribution of this paper. We have studied the aliasing effect in two sampling schemes: random deviations from regular times and sampling times generated by the IPFM model. The former was chosen as a basic model which allowed to explain a mechanism of aliasing suppression. The later is inspired by a model used to describe generation of events in biological systems, such as the sinus node or the neuronal firing.

The two essentially distinctive sampling models required different way of analysis. The random model analysis relies on terms common in probability/statistics. We have shown that the characteristic function of timing jitter plays a key role in composition of the resulting spectrum. The characteristic function modulates amplitudes of aliased or imaged spectral lines. Besides the line portion of the spectrum, a smooth noise portion shaped by a term of characteristic function is present due to a random nature of the sampling process, despite the sampled signal is deterministic. The IPFM sampling scheme uses concepts originated in communication systems theory and presented analysis is purely deterministic. Consequently, the smooth spectrum is not present on the other hand, spurious peaks frequently occur due to a complex nature of the spectrum of the modulated signal. In the both situations, the relative reduction of aliased components was evaluated by a quantity named in this paper as the aliasing reduction index. Interestingly, the derived formula approximating this index for small-to-moderate sampling irregularity is identical for the both sampling models. The formula which we have found seems to be quite universal and relates a decrease of unwanted aliased components with the standard deviation of timing irregularity measured as the fraction of the mean sampling period. Validity of the formula in more general sampling schemes, such as those with dependent sampling irregularities and additive random sampling, is intended to be studied in a future work.

Mitigation of the aliasing effect by means of nonuniform sampling not only attracts scientific interest but also has practical implications. We have demonstrated the possibility to detect frequencies above the Nyquist limit in the ECG-derived respiration signal and the heart rate signal. The condition of aliasing, so-called cardiac aliasing, was elicited by increasing the respiration rate, and signals derived from recorded electrocardiogram were analyzed by different methods. Conventionally used interpolation/resampling methods were found to be unsuitable in the condition of cardiac aliasing. The Lomb-Scarglre method, the classical periodogram used for nonuniform samples and the spectrum of counts are all capable to uncover components above the Nyquist frequency. Since real-world measurements are usually contaminated by noise, a signal to be detected needs not to be a pure sinusoid, or it can be random and even nonstationary, development of a universal method able to resolve between true and aliased spectral components, with properly assigned significance measure, is another challenging task intended for future work.

مراجع

  1. P. Babu and P. Stoica, “Spectral analysis of nonuniformly sampled data𠅊 review,” Digital Signal Processing, vol. 20, no. 2, pp. 359–378, 2010. View at: Publisher Site | Google Scholar
  2. N. R. Lomb, “Least-squares frequency analysis of unequally spaced data,” الفيزياء الفلكية وعلوم الفضاء, vol. 39, no. 2, pp. 447–462, 1976. View at: Publisher Site | Google Scholar
  3. P. Stoica, J. Li, and H. He, “Spectral analysis of nonuniformly sampled data: a new approach versus the periodogram,” IEEE Transactions on Signal Processing, vol. 57, no. 3, pp. 843–858, 2009. View at: Publisher Site | Google Scholar
  4. J. D. Scargle, “Studies in astronomical time series analysis. II—statistical aspects of spectral analysis of unevenly spaced data,” Astrophysical Journal, vol. 263, pp. 835–853, 1982. View at: Publisher Site | Google Scholar
  5. P. Vaníპk, “Approximate spectral analysis by least-squares fit-Successive spectral analysis,” الفيزياء الفلكية وعلوم الفضاء, vol. 4, no. 4, pp. 387–391, 1969. View at: Publisher Site | Google Scholar
  6. P. Laguna, G. B. Moody, and R. G. Mark, “Power spectral density of unevenly sampled data by least-square analysis: performance and application to heart rate signals,” IEEE Transactions on Biomedical Engineering, vol. 45, no. 6, pp. 698–715, 1998. View at: Publisher Site | Google Scholar
  7. “Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology,” Circulation, vol. 93, no. 5, pp. 1043–1065, 1996. View at: Google Scholar
  8. G. G. Berntson, J. T. Bigger, D. L. Eckberg et al., “Heart rate variability: origins, methods, and interpretive caveats,” Psychophysiology, vol. 34, no. 6, pp. 623–648, 1997. View at: Google Scholar
  9. M. Teplan, L. Molლn, and M. Zeman, “Spectral analysis of cardiovascular parameters of rats under irregular light-dark regime,” in Proceedings of the 8th International Conference on Measurement, pp. 343–346, Smolenice, Slovak Republic, 2011. View at: Google Scholar
  10. V. V. Vityazev, “Time series analysis of unequally spaced data: intercomparison between estimators of the power spectrum,” in Proceedings of the Astronomical Data Analysis Software and Systems VI ASP Conference Series, vol. 125, pp. 166–169, 1997. View at: Google Scholar
  11. A. Schwarzenberg-Czerny, “The distribution of empirical periodograms: Lomb-Scargle and PDM spectra,” Monthly Notices of the Royal Astronomical Society, vol. 301, no. 3, pp. 831–840, 1998. View at: Google Scholar
  12. J. Mateo and P. Laguna, “Improved heart rate variability signal analysis from the beat occurrence times according to the IPFM model,” IEEE Transactions on Biomedical Engineering, vol. 47, no. 8, pp. 985–996, 2000. View at: Publisher Site | Google Scholar
  13. M. Brennan, M. Palaniswami, and P. Kamen, “Distortion properties of the interval spectrum of IPFM generated heartbeats for heart rate variability analysis,” IEEE Transactions on Biomedical Engineering, vol. 48, no. 11, pp. 1251–1264, 2001. View at: Publisher Site | Google Scholar
  14. J. Púčik, O. Ondráპk, E. Cocherová, and A. Sultan, “Spectrum of counts computation for HRV analysis,” in Proceedings of 19th International Conference Radioelektronika (RADIOELEKTRONIKA '09), pp. 255–258, April 2009. View at: Publisher Site | Google Scholar
  15. E. J. Bayly, “Spectral analysis of pulse frequency modulation in the nervous systems,” IEEE Transactions on Biomedical Engineering, vol. 15, no. 4, pp. 257–265, 1968. View at: Google Scholar
  16. G. Moody, R. Mark, A. Zoccola, and S. Mantero, “Derivation of respiratory signals from multi-lead ECGs,” IEEE Computer Society, vol. 12, pp. 113–116, 1985. View at: Google Scholar
  17. G. D. Clifford, F. Azuaje, P. E. McSharry et al., Advanced Methods and Tools for ECG Data Analysis, Artech House, 2006.
  18. D. Widjaja, J. Taelman, S. Vandeput et al., “ECG-derived respiration: comparison and new measures for respiratory variability,” in Proceedings of the Computing in Cardiology (CinC '10), pp. 149–152, September 2010. View at: Google Scholar
  19. A. Gersten, O. Gersten, A. Ronen, and Y. Cassuto, “The RR interval spectrum, the ECG signal and aliasing,” . Eprint, http://arxiv.org/abs/physics/9911017v1. View at: Google Scholar
  20. J. Šurda, S. Lovás, J. Púčik, and M. Jus, “Spectral properties of ECG signal,” in Proceedings of the of the 17th International Conference Radioelektronika, pp. 537–541, Brno, Czech Republic, 2007. View at: Google Scholar
  21. E. Toledo, I. Pinhas, D. Aravot, and S. Akselrod, “Very high frequency oscillations in the heart rate and blood pressure of heart transplant patients,” Medical and Biological Engineering and Computing, vol. 41, no. 4, pp. 432–438, 2003. View at: Google Scholar
  22. H. A. Campbell, J. Z. Klepacki, and S. Egginton, “A new method in applying power spectral statistics to examine cardio-respiratory interactions in fish,” Journal of Theoretical Biology, vol. 241, no. 2, pp. 410–419, 2006. View at: Publisher Site | Google Scholar
  23. J. Púčik, M. Uhrík, A. Sultan, and A. Šurda, “Experimental setup for cardio-respiratory interaction study,” in Proceedings of the 8th Czech-Slovak Conference, Trends in Biomedical Engineering, pp. 126–129, Bratislava, Slovakia, 2009. View at: Google Scholar
  24. R. Berger, S. Akselrod, D. Gordon, and R. Cohen, “An efficient algorithm for spectral analysis of heart rate variability,” IEEE Transactions on Bio-Medical Engineering, vol. 33, no. 9, pp. 900–904, 1986. View at: Google Scholar
  25. W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling et al., Numerical Recipes in Fortran 77: The Art of Scientific Computing, vol. 1, Cambridge University Press, 1992.
  26. G. L. Bretthorst, “Nonuniform sampling: bandwidth and aliasing,” Concepts in Magnetic Resonance A, vol. 32, no. 6, pp. 417–435, 2008. View at: Publisher Site | Google Scholar

Copyright

Copyright © 2012 Jozef Púčik et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Note: This Plot Contains an Error¶

After the book was in press, a reader pointed out that this plot contains a typo. Instead of passing the noisy data to the Lomb-Scargle routine, we had passed the underlying, non-noisy data. This caused an over-estimate of the Lomb-Scargle power.

Because of this, we add two extra plots to this script: the first reproduces the current plot without the typo. In it, we see that for the noisy data, the period is not detected for just

30 points within ten periods.

In the second additional plot, we increase the baseline and the number of points by a factor of ten. With this configuration, the peak is detected, and the qualitative aspects of the above discussion hold true.


1. Introduction

The functioning of eukaryotic cells is controlled by accurate timing of biological cycles, such as cell cycles and circadian rhythms. These are composed of an echelon of molecular events and checkpoints. At the transcription level, these events can be quantitatively observed by measuring the concentration of messenger RNA (mRNA), which is transcribed from DNA and serves as the template for synthesizing protein. To achieve this goal, in the microarray experiments, high-throughput gene chips are exploited to measure genome-wide gene expressions sequentially at discrete time points. These time series data have three characteristics. Firstly, most data sets are of small sample size, usually not more than 50 data points. Large sample sizes are not financially affordable due to high cost of gene chips. Also the cell cultures lose their synchronization and render meaningless data after a period of time. Secondly, the data are usually evenly sampled and have many time points missing. Thirdly, most data sets are customarily corrupted by experimental noise and the produced uncertainty should be addressed in a stochastic framework.

Extensive genome-wide time course microarray experiments have been conducted on organisms such as Saccharomyces cerevisiae (budding yeast) [1], human Hela [2], and Drosophila melanogaster (fruit fly) [3]. Budding yeast in [1] has served as the predominant data source for various statistical methods in search of periodically expressed genes, mainly due to its pioneering publication and relatively larger sample size compared with its peers. By assuming the signal in the cell cycle to be a simple sinusoid, Spellman et al. [1] and Whitfield et al. [2] performed a Fourier transformation on the data sampled with different synchronization methods, while Giurcaneanu [4] explored the stochastic complexity of the detection mechanism of periodically expressed genes by means of generalized Gaussian distributions. Ahdesmäki et al. [5] implemented a robust periodicity testing procedure also based on the non-Gaussian noise assumption. Alternatively, Luan and Li [6] employed guide genes and constructed cubic B-spline-based periodic functions for modeling, while Lu et al. [7] employed up to three harmonics to fit the data and proposed a periodic normal mixture model. Power spectral density estimation schemes have also been employed. Wichert et al. [8] applied the traditional periodogram on various data sets. Bowles et al. [9] compared Capon and robust Capon methods in terms of their ability to identify a predetermined frequency using evenly sampled data sets, under the assumption of a known period. Lichtenberg et al. [10] compared [167] while proposing a new score by combining the periodicity and regulation magnitude. The majority of these works dealt with evenly sampled data. When missing data points were present, either the vacancies were filled by interpolation in time domain, or the genes were discarded if there were more than 30% data samples missing.

Biological experiments generally output unequally spaced measurements. The major reasons are experimental constraints and event-driven observation. The rate of measurement is directly proportional to the occurrence of events. Therefore, an analysis based on unevenly sampled data is practically desired and technically more challenging. While providing modern spectral estimation methods for stationary processes with complete and evenly sampled data [11], the signal processing literature has witnessed an increased interest in analyzing unevenly sampled data sets, especially in astronomy, in the last decades. The harmonics exploited in discrete Fourier transform (DFT) are no longer orthogonal for uneven sampling. However, Lomb [12] and Scargle [13] demonstrated that a phase shift suffices to make the sine and cosine terms orthogonal. The Lomb-Scargle scheme has been exploited in analyzing the budding yeast data set by Glynn et al. [14]. Schwarzenberg-Czerny [15] employed one-way analysis of variance (AoV) and formulated an AoV periodogram as a method to detect sharp periodicities. However, it relies on an infeasible biological assumption, that is, the observation duration covers many cycles. Along with this line of research, Ahdesmäki et al. [16] proposed to use robust regression techniques, while Stoica and Sandgren [17] updated the traditional Capon method to cope with the irregularly sampled data. Wang et al. [18] reported a novel technique, referred to as the missing-data amplitude and phase estimation (MAPES) approach, which estimates the missing data and spectra iteratively through the expectation maximization (EM) algorithm. In general, Capon and MAPES methods possess a better spectral resolution than Lomb-Scargle periodogram. In this paper, we propose to analyze the performance of three of the most representative spectral estimation methods: Lomb-Scargle periodogram, Capon method, and the MAPES technique in the presence of missing samples and irregularly spaced samples. The following questions are to be answered in this study: do technically more sophisticated schemes, such as MAPES, achieve a better performance on real biological data sets than on simpler schemes? Is the efficiency sacrificed in using these advanced methods justifiable?

The remainder of this paper is structured as follows. In Section 2, we introduce the three spectral analysis methods, that is, Lomb-Scargle, Capon and MAPES. Hypothesis tests for periodicity detection and the corresponding -values are also formulated. The multiple testing correction is discussed. Section 3 presents simulation results. The performances of the three schemes are compared based on published cell-cycle and noncell-cycle genes of the Saccharomyces cerevisiae (budding yeast). Then the spectral analysis for the data set of Drosophila melanogaster (fruit fly) is performed, and a list of 149 genes are presented as cycle-related genes. The synchronization effects are also considered. Concluding remarks and future works constitute the last section, and full results are provided online in the supplementary materials [19].


Time-uncertain spectral analysis

Now let’s consider age uncertainties.

The GeoChronR approach to quantifying and propagating those uncertainties is to leverage the power of ensembles. Here we will illustrate this with MTM. Let us repeat the previous analysis by looping over the 1,000 age realizations output by the tuning algorithm HMM-Match. That means computing 1000 spectra. We’ve already seen that nuspectral is too slow for an ensemble job, so let’s leave it out. Also, since REDFIT is a tapered version of Lomb-Scargle, and comes with more features, let’s focus here on comparing an REDFIT and MTM in ensemble mode.

This took a few seconds with 1000 ensemble members, and the resulting output is close to 10 Mb. Not an issue for modern computers, but you can see why people weren’t doing this in the 70’s, even if they wanted to. Now, let’s plot the result:

To this we can add the periods identified as significant owing to MTM’s harmonic ratio test. geoChronR deals with it by computing at each frequency the fraction of ensemble members that exhibit a significant peak. One simple criterion for gauging the level of support for such peaks given age uncertainties is to pick out those periodicities that are identified as significant more than 50% of the time.

One could of course, impose a more stringent criterion (e.g. 80%, 90%, etc). To do that, specify which(mtm.ens.spec$sig.freq>0.8) or which(mtm.ens.spec$sig.freq>0.9) . That relatively close peaks appear significant may be an artifact of neglecting the multiple hypothesis problem (cf Vaughan et al, 2011).

So what do we find? Consistent with previous investigations (e.g. Mudelsee et al 2009), the effect of age uncertainties is felt more at high frequencies, with the effect of broadening peaks, or lumping them together. Again, the peaks that rise above the power-law background are the

40kyr peaks, which is not surprising. There is significant energy around the precessional peaks, but it is rather spread out, and hard to tie to particular harmonics. No doubt a record with higher-resolution and/or tighter chronology would show sharper precessional peaks.

Do we get the same answer in REDFIT? Let’s find out.

again, it took a bit of time. let’s plot the result:

The result is quite a bit smoother here, perhaps because of the choice of parameters.


Making sense of the lomb-scargle periodogram - Astronomy

Context. Although timing variations in close binary systems have been studied for a long time, their underlying causes are still unclear. A possible explanation is the so-called Applegate mechanism, where a strong, variable magnetic field can periodically change the gravitational quadrupole moment of a stellar component, thus causing observable period changes. One of the systems exhibiting such strong orbital variations is the RS CVn binary HR 1099, whose activity cycle has been studied by various authors via photospheric and chromospheric activity indicators, resulting in contradicting periods.
Aims: We aim at independently determining the magnetic activity cycle of HR 1099 using archival X-ray data to allow for a comparison to orbital period variations.
Methods: Archival X-ray data from 80 different observations of HR 1099 acquired with 12 different X-ray facilities and covering almost four decades were used to determine X-ray fluxes in the energy range of 2-10 keV via spectral fitting and flux conversion. Via the Lomb-Scargle periodogram we analyze the resulting long-term X-ray light curve to search for periodicities.
Results: We do not detect any statistically significant periodicities within the X-ray data. An analysis of optical data of HR 1099 shows that the derivation of such periods is strongly dependent on the time coverage of available data, since the observed optical variations strongly deviate from a pure sine wave. We argue that this offers an explanation as to why other authors derive such a wide range of activity cycle periods based on optical data. We furthermore show that X-ray and optical variations are correlated in the sense that the star tends to be optically fainter when it is X-ray bright.
Conclusions: We conclude that our analysis constitutes, to our knowledge, the longest stellar X-ray activity light curve acquired to date, yet the still rather sparse sampling of the X-ray data, along with stochastic flaring activity, does not allow for the independent determination of an X-ray activity cycle.


شاهد الفيديو: تحديد التبويض للدورة المنتضمة و الغير منتضمة بدقة. الأيام المناسبة للحمل (أغسطس 2022).