قرار است تصور حرکتی1 را از سیگنال EEG تفکیک کنیم. با این روش می توانیم از سیگنال حاصله جهت استفاده در به حرکت درآوردن بازوهای مکانیکی استفاده کنیم. برای مثال می‌توان با ایجاد تصور حرکت دست‌ها و تفکیک و انتقال این سیگنال به فردی که از ناحیه دست قطع عضو گردیده یا نرون‌های عصبی آن از کار افتاده، دست مصنوعی و یا دست خود فرد را به حرکت درآوریم.
۱. مقدمه
در این گزارش، روش CSP به همراه سه تعمیم آن یعنی ACSP، ACCSP و SUTCCSP به طور کامل شرح داده می‌شود. الگو های فضایی مشترک روشی محبوب در جداسازی بین کلاس‌ها است. اعمال این روش به داده‌ها پیش از کلاس‌بندی، عملکرد را بالا تر می‌برد. این روش ماهیتی چند کاناله دارد و اغلب برای چند کانال سیگنال به کار می‌رود. یکی از مشکلات این روش عدم دخالت اطلاعات مربوط به فاز سیگنال‌ها در آن است. سیگنال‌های جهان ما اکثراً غیر‌دایره ای هستند. برای مثال در سیگنال EEG ، به دلیل اختلاف توان بین کانال‌ها و همبستگی بین آنها این موضوع رویت می‌شود. سه تعمیمی که مورد بحث قرار می‌گیرند، سعی در اضافه کردن اطلاعات فاز و جدا سازی بیشتر و بهتر بین کلاس‌ها را دارند.
۲. کارهای مرتبط
روش الگو های فضایی مشترک 2 روشی محبوب و پر کاربرد در جدا سازی بین داده‌های دو کلاس، مخصوصاً در تفکیک فعالیت‌های ذهنی در واسط‌های مغز کامپیوتر 3 مبتنی بر داده‌های EEG 4 است. عملکرد بالا در کلاس‌بندی و فیلترهای فضایی مفید به دست آمده از این روش، این الگوریتم را مشهور ساخته است [1] و [2].
این الگوریتم برای اولین بار برای تشخیص غیر طبیعی بودن 5 سیگنال EEG به کار گرفته شد [3] و سپس در جدا سازی بین کلاس‌های مربوط به واسط‌های مغز کامپیوتر و الگو های حرکتی 6 نیز مورد استفاده قرار گرفت [1] و [2]. به دلیل ماهیت چند کاناله بودن، این روش کاربرد فراوانی برای سیگنال های EEG دارد.
روش CSP، با اعمال فیلترهای فضایی به ورودی‌ها، واریانس سیگنال‌ها را در کلاس اول ماکزیمم و به طور همزمان در کلاس دیگر مینیمم می کند و سپس از سیگنال‌های فیلتر شده ویژگی‌های کلاس اول را می‌سازد. این کار به طور عکس نیز اتفاق می افتد. یک سری از فیلترهای فضایی، واریانس سیگنال را در کلاس دوم ماکزیمم و به طور همزمان در کلاس اول مینیمم می‌کنند. سپس با استفاده از این سیگنال‌های فیلتر شده، ویژگی‌های کلاس دوم ایجاد می‌شود. برای به دست آوردن این فیلترهای فضایی، CSP مجموع ماتریس‌های کواریانس دو کلاس را تبدیل به ماتریس همانی 7 میکند. این کار با اعمال تبدیلی به نام سفید کنندگی8 انجام می‌شود. با این کار، افزایش واریانس یک کلاس، واریانس کلاس دیگر را کاهش می‌دهد.
تعمیم‌های زیادی از این روش ارائه شده است که هر کدام سعی در بهبود کارایی و یا رفع محدودیتی در الگوریتم اصلی دارند. برای مثال در [5] حالت چند کلاسه‌ی CSP و در [6] انتخاب فیلتر طیفی به طور خود کار9 مطرح شده است. در این گزارش، به غیر از روش CSP ، سه تعمیم آن، یعنی ACSP 10 ، ACCSP 11 و SUTCCSP12 نیز شبیه سازی می شود.
این سه تعمیم، به غیر از اطلاعات دامنه‌ی سیگنال‌ها، اطلاعات فاز سیگنال‌ها را نیز در تفکیک کلاس‌ها دخیل می‌کنند. چرا که بسیاری از سیگنال‌های دنیای واقعی، سیگنال‌های غیر دایره ای 13 هستند. به عبارت دیگر، بین سیگنال‌های دو کانال، همبستگی 14 وجود دارد. روش ACSP فقط اطلاعات فاز سیگنال ها را در الگوریتم CSP دخیل میکند. ولی روش‌های ACCSP و SUTCCSP به غیر از این کار، از آمار درجه دوم مربوط به مجموعه اعداد مختلط15 استفاده می کنند. مباحث آماری که به ما می گوید، برای اعداد متغیر های تصادفی مختلف به غیر از کواریانس، به شبه کواریانس 16 هم نیاز داریم. شبه کواریانس داده های غیر دایره ای صفر نیستند! روش ACCSP ماتریسی حاصل از ترکیب کواریانس و شبه کواریانس را برای قطری‌سازی 17 استفاده میکند. در حالی که روش SUTCCSP ، به طور همزمان ماتریس های کواریانس و شبه کواریانس را قطری می سازد.
۳. الگو های فضایی مشترک برای سیگنال های حقیقی و تحلیلی
الگوریتم CSP با روش CSP مختلط یا همان ACSP18 تفاوتی ندارد. ولی در روش ACSP سیگنال های ورودی پیش از اعمال به الگوریتم، به صورت تحلیلی19 در می آیند. پس ابتدا فرض کنید که سیگنال ورودی، سیگنالی حقیقی است.
ماتریس داده های { Z }_{ a }\in { R }^{ N\times T } و { Z }_{ b } را در نظر بگیرید. { Z }_{ a } و { Z }_{ b } دو ماتریس با میانگین صفر هستند. N تعداد کانال ها و T تعداد نمونه ها در هر کانال است. در حقیقت N کانال داریم که سیگنال هر کانال T نمونه دارد و میانگین صفر است. ماتریس کواریانس { Z }_{ a } و { Z }_{ b } را به صورت زیر تعریف می کنیم:
که در آن منظور از E\left[ \cdot \right]،عملگر امید ریاضی و { \left( \cdot \right) }^{ H } بیانگر ترانهاده ی مزدوج مختلط یا همان هرمیتین است. ماتریس های کواریانس، ماتریس های N\times N هستند.
در بسیاری از مقالات به جای استفاده از عملگر امید ریاضی از تخمین لحظه ای آن استفاده می کنند و هم چنین، کواریانس را نرمالیزه می کنند تا مقادیر آن بین صفر تا یک باشد. یعنی:
که در آن منظور از tr\left( \cdot \right) ، جمع المان های قطری ماتریس است. می دانیم که محاسبه ی امید ریاضی با توجه به نیاز به تابع چگالی داده ها، کار مشکلی است. به همین دلیل اغلب آن را تخمین می زنند. نرمالیزه کردن ماتریس کواریانس به وسیله ی tr\left( \cdot \right) آن ماتریس انجام می شود. چرا که تریس در ماتریس کواریانس یک داده، واریانس کل آن داده است.
در مسائل شناسایی آماری الگو، اصولاً در داده های آموزشی، بیش از یک داده برای هر کلاس وجود دارد. پس در این جا ( و در تمامی تعمیم های CSP ) هم برای هر کدام از داده ها، ماتریس کواریانس مربوطه را به دست آورده و سپس میانگین ماتریس های مربوط به هر کلاس را به دست می آوریم. پس دو ماتریس { \bar { C } }_{ a } و { \bar { C } }_{ b } که N\times N ،هرمیتین و تقریباً همیشه مثبت هستند به دست می آید. اگر داده های ما مختلط باشند، ماتریس های کواریانس دارای مقادیر حقیقی روی قطر اصلی و مقادیر مختلط در سایر المان های ماتریس هستند.
حال یک ماتریس مرکب به صورت زیر تعریف می کنیم:
و { C }_{ c }ترکیبی از ماتریس های کواریانس هر دو کلاس است. می توان این ماتریس را به صورت زیر تبدیل کرد:
و { U }_{ c } ماتریس بردار های ویژه است که هر ستون آن را یک بردار ویژه تشکیل می دهد. بردار های ویژه ای که هر کدام مربوط به یکی از مقادیر ویژه ی موجود در ماتریس قطری { \Lambda }_{ c } ( ماتریس مقادیر ویژه ) هستند. در ضمن { U }_{ c } ، یک ماتریس واحد است. یعنی رابطه ی { U }_{ c }{ U }_{ c }^{ H } = I وجود دارد.
به تبدیل G =\sqrt { { \Lambda }_{ c }^{ -1 } } { U }_{ c }^{ H } ، تبدیل سفید کننده می گویند. چرا که:
یعنی این تبدیل، کواریانس مرکب را قطری ( ماتریس همانی ) می کند. پس می توان نوشت:
که در آن I نشان دهنده ی ماتریس همانی است.
حال اگر{ S }_{ a }=G{ \bar { C } }_{ a }{ G }^{ H } و { S }_{b }=G{ \bar { C } }_{ b }{ G }^{ H } را تعریف کنیم و B یک ماتریس متعامد نرمال باشد، آنگاه داریم:
چرا که می دانیم{ B }^{ H }B=I.
یعنی می توان گفت که ماتریس های { S }_{ a } و { S }_{ b } ماتریس بردار های ویژه ی مشترکی دارند. به طوری که:
حال با توجه به رابطه ی (10) ، اگر ماتریس مقادیر ویژه ی { \Lambda }_{ a } را به صورت نزولی مرتب کنیم، به طور هم زمان { \Lambda }_{ b } به صورت صعودی مرتب می شود. با این کار، بردار ویژه ای از ماتریس B که مربوط به بزرگ ترین مقدار ویژه ی { \Lambda }_{ a } می شود،مربوط به کوچک ترین مقدار ویژه ی { \Lambda }_{ b } است. یعنی یک بردار ویژه از B (مربوط به بزرگترین مقدار ویژه ی { \Lambda }_{a } )، به طور همزمان واریانس کلاس a را ماکزیمم و واریانس کلاس b را مینیمم می کند. پس این فیلتر فضایی را به صورت زیر تعریف می کنیم:
و با توجه به این تعریف، رابطه ی (9) و (10) داریم:
حال ماتریس داده ها را روی این فیلتر فضایی تصویر می کنیم تا داده های جدیدی به دست آوریم. در حقیقت داده های ورودی را توسط فیلتر های فضایی فیلتر می کنیم. یعنی:
هر کدام از ستون های ماتریس W=\left( { w }_{ 1 },{ w }_{ 2 },...,{ w }_{ N } \right)را یک فیلتر فضایی می نامند.
برای جدا سازی بین دو کلاس، از واریانس سیگنال های فیلتر شده در رابطه ی (12) استفاده می کنیم. می دانیم که هر سطر از ماتریس V ، یک سیگنال جدید است. به منظور کاهش بعد فضای ویژگی، m سطر اول V ( که وابسته به m مقدار ویژه ی بزرگ { \Lambda }_{ a } است، یعنی بیشترین واریانس را در کلاس a ایجاد می کنند ) و m سطر آخر V ( که وابسته به m مقدار ویژه ی بزرگ { \Lambda }_{ b } است یعنی بیشترین واریانس را در کلاس b ایجاد می کنند ) را انتخاب می کنیم. به عبارت دیگر، سطر های ماتریس V={ \left( { v }_{ 1 },...,{ v }_{ m },{ v }_{ N-m+1 },...,{ v }_{ N } \right) }^{ T } ، اختلاف واریانس بین دو کلاس را ماکزیمم می کنند. چرا که m سطر اول بیشترین واریانس در کلاس a و کمترین واریانس در کلاس b را دارند و m سطر آخر، دارای بیشترین واریانس در کلاس b و کمترین واریانس در کلاس a هستند. ویژگی های مورد نظر ما از رابطه ی زیر به دست می آید:
که در این رابطه p ، ویژگی مربوط به سطر pاُم ماتریس و منظور از var\left( { v }_{ p } \right) ، واریانس سطر pاُم ماتریس V است. در این رابطه، به منظور نرمالیزاسیون ، واریانس هر سطر ( سیگنال فیلتر شده ی هر کانال ) را به مجموع واریانس سطر ها تقسیم می کنیم. همچنین از لگاریتم برای تخمین توزیع نرمال داده ها استفاده می کنیم [4] .
حال بردار ویژگی مربوط به هر کلاس را از روی داده های آزمایشی به دست می آوریم و از آن ها برای آموزش کلاسیفایر مورد نظر استفاده می کنیم.
۴. تبدیل سیگنال حقیقی به سیگنال تحلیلی
می دانیم که سیگنال‌هایی که در طبیعت وجود دارند، داده هایی با مقادیر حقیقی هستند. پس در الگوریتم CSP ، عملگر { \left( \cdot \right) }^{ H } به عملگر ترانهاده یا { \left( \cdot \right) }^{ T } تبدیل می شود.
در 2010 ، فالزون روش ACSP پیشنهاد کرد [12] . در این روش، سیگنال های ورودی پیش از اعمال به الگوریتم CSP ، تبدیل به فُرم تحلیلی خود می شوند. پس ورودی Z ، ماتریسی با مقادیر مختلط است. در این جا هم هدف قطری سازی ماتریس کواریانس مرکب است. با این تفاوت که در این حالت، ماتریس های کواریانس، ماتریس بردار های ویژه و ضرایب فیلتر به دست آمده، همگی مختلط هستند. اما ماتریس مقادیر ویژه ی کواریانس مرکب \left( { \Lambda }_{ C } \right)مقادیر حقیقی دارد [8] .
در حقیقت، روش ACSP از اطلاعات فاز ورودی نیز برای متمایز کردن کلاس ها استفاده می کند. چرا که سیگنال های ورودی به فرم تحلیلی خود در آمده و اطلاعات فاز را نیز در بر دارد.
اگر فرض کنیم که x\left( t \right)سیگنال ورودی باشد، سیگنال تحلیلی مربوطه، z\left( t \right)، به صورت زیر تعریف می شود:
و \bar { x } \left( t \right) ، تبدیل هیلبرت سیگنال ورودی است:
که منظور از P.V. ، مقدار اصلی کوشی ( روشی برای مقدار دادن به انتگرال های ناسره ) است.
در خصوص استفاده از تبدیل هیلبرت باید به این موضوع توجه داشت که این تبدیل فقط برای سیگنال هایی با باند فرکانسی باریک استفاده می شود و نمی تواند برای سایر سیگنال ها، اطلاعات کاملی راجع به محتوای فرکانسی آن ها ارائه کند. همچنین زمانی که تغییرات دامنه سریع تر از تغییرات فاز است (مثلاً سیگنال هایی که میانگین صفر ندارند)، تبدیل هیلبرت، فرکانس لحظه ای منفی ایجاد می کند [9] . به همین جهت، برای رفع محدودیت پهنای باند تبدیل هیلبرت، باید از تبدیل فوریه و یا EMD در کنار ACSP استفاده کرد [5] .
طبق [5] ، المان های بردار ویژگی به دست آمده از ACSP برای هر کلاس را می توان از روابط زیر به دست آورد:
۵. الگو های فضایی مشترک مختلط افزوده (ACCSP)
برای یک بردار میانگین صفر مختلط مانند z ، ماتریس کواریانس به صورت \complement =E\left[ z{ z }^{ H } \right] محاسبه می شود. مباحث آماری مربوط به اعداد مختلط، به طور کامل تعمیمی از مباحث آماری اعداد حقیقی نیستند! پس ماتریس کواریانس \complement به طور کامل آمار درجه دوم بردار z را توصیف نمی کند!! پس برای برداری مانند z که مقادیر مختلط دارد، به ماتریس شبه کواریانس یعنی P=E\left[ z{ z }^{ T } \right] نیز احتیاج داریم.
متغیر تصادفی z={ z }_{ r }+{ jz }_{ i } را در نظر بگیرید. کواریانس این ماتریس برابر E\left[ z{ z }^{ \ast } \right] =E\left[ { \left| z \right| }^{ 2 } \right] =E\left[ { z }_{ r }^{ 2 }+{ z }_{ i }^{ 2 } \right] می شود. این در حالیست که شبه کواریانس z که به صورت E\left[ z{ z }^{ \ast } \right] =E\left[ { z }_{ r }^{ 2 }-{ z }_{ i }^{ 2 }+2j{ z }_{ r }{ z }_{ i } \right] =E\left[ { z }_{ r }^{ 2 } \right] -E\left[ { z }_{ i }^{ 2 } \right] +2jE\left[ { z }_{ r }{ z }_{ i } \right] تعریف می شود، هنگامی بی اثر است که به طور همزمان، { z }_{ r } و { z }_{ r } با یک دیگر ناهمبسته بوده و واریانس آن ها برابر باشد. به سیگنال هایی که شبه کواریانس صفر دارند، دایره ا ی مرتبه دوم یا مناسب می گویند. با این حال، حتی اگر داده ای در دنیای واقعی، دایره ا ی باشد، به دلیل مشاهده ی سیگنال در پنجره های کوچک، نویز های نا همسانگر و توان متفاوت کانال های مختلف، شبه کواریانس برابر صفر نیست [5] . شکل 1، نمایشی هندسی از توضیح داده های دایره ا ی و غیر دایره ا ی را نشان می دهد. میزان دایره ا ی بودن بستگی به میزان همبستگی بین قسمت حقیقی و موهومی داده ها دارد.همبستگی کمتر، دایره ا یبودن بیشتر و همبستگی بیشتر، غیر دایره ا یبودن بیشتری را ایجاد می کند [5] .
شکل 1: نمایش هندسی داده ها از نظر دایره ای و غیر دایره ای بودن. نمودار پراکندگی قسمت حقیقی در برابر قسمت موهومی دو نوع داده نمایش داده شده است. (a) داده ی دایره ای که همبستگی میان قسمت حقیقی و موهومی آن صفر است (b) داده ی غیر دایره ای که همبستگی میان قسمت حقیقی و موهومی آن 0.9 است.
از این رو، برای داده های مختلط، باید خواص آماری مشترک z و { z }^{ \ast } مورد بررسی قرار گیرد. پس فُرم افزوده ای از متغیر مختلط z به صورت \hat { z } ={ \begin{bmatrix} { z }^{ T } & { z }^{ H } \end{bmatrix} }^{ T }، تعریف می شود. ماتریس کواریانس این متغیر افزوده ( ماتریس افزوده ) به صورت زیر است:
این ماتریس کواریانس افزوده، خصوصیات آماری مرتبه ی دوم کامل کواریانس و شبه کواریانس را در بر دارد.
روش ACSP ( و به طور قطع روش CSP )، غیر دایره ای بودن داده ها را در نظر نمی گیرد. به همین دلیل روش الگوهای فضایی مشترک مختلط افزوده (ACCSP) مطرح می شود [5] . در این روش به جای ماتریس کواریانس، فیلتر های فضایی از طریق ماتریس کواریانس افزوده به دست می آید. ماتریس داده های تبدیل شده به فرم تحلیلی { Z }_{ b }\in { \complement }^{ N\times T } و { Z }_{ a } را داریم ( منظور از \complement ، مجموعه اعداد مختلط است ). این ماتریس ها را به فرم افزوده ی خود { \hat { Z } }_{ a } و { \hat { Z } }_{ b } تبدیل کرده و کواریانس آنها را به دست می آوریم:
ادامه ی روند، مانند الگوریتم اولیه ی CSP است. ماتریس کواریانس مرکب افزوده تعریف می شود. تبدیل سفید کننده روی آن اعمال می‌شود و ماتریس بردار های ویژه ی مشترک بین کلاس a و b به دست می آید. ماتریس فیلتر فضایی افزوده،{ \hat { W } }^{ H }، از رابطه ی (11) به دست می آید. حال برای هر داده ی ورودی Z ، فرم افزوده ی ورودی،\hat { Z }، را به دست می آوریم و سیگنال ها را با استفاده از رابطه ی زیر، فیلتر می کنیم:
ماتریس\hat { V }مقادیر مختلط دارد و مانند قبل، m سطر اول و m سطر آخر آن انتخاب می شود و از رابطه ی (13) و یا (17) و (18) برای استخراج ویژگی از این 2m سیگنال، استفاده می شود.
۶. الگو های فضایی مشترک مختلط افزوده با تبدیل ناهمبستگی قوی (SUTCCSP)
در این قسمت تبدیل ناهمبستگی قوی (SUT) ) را به کار گرفته می شود. این تبدیل، تعمیمی از تبدیل سفید کنندگی (7) برای متغیر های تصادفی مختلط غیر دایره ای ا ست. در این تعمیم از CSP که به طور مخفف SUTCCSP نامیده می شود، هردوی ماتریس های کواریانس و شبه کواریانس به طور همزمان قطری می شوند.
ما در SUT به دنبال ماتریس Q می گردیم به طوری که:
و C ماتریس کواریانس، P ماتریس شبه کواریانس و Λ ماتریسی قطری است. به عبارت دیگر، SUT به طور همزمان، ماتریس کواریانس و شبه کواریانس را قطری می کند.
فرض کنید ماتریس داده های تبدیل شده به فرم تحلیلی { Z }_{ a },{ Z }_{ b }\epsilon { C }^{ N\times T } را در اختیار داریم. این بار علاوه بر تعریف ماتریس کواریانس مرکب، ماتریس شبه کواریانس مرکب را نیز تعریف می کنیم. به طوری که:
تبدیل سفید کنندگی را به { C }_{ c } اعمال می کنیم. می دانیم که { U }_{ c } و { \Lambda }_{ c } به ترتیب، ماتریس بردار های ویژه و ماتریس مقادیر ویژه ی { C }_{ c } هستند. همان طور که بیان شد نتیجه ی اعمال این تبدیل، رابطه ی(7)است.
حال فرض کنید تبدیل سفید کننده را به صورت زیر به شبه کواریانس اعمال کنیم:
ماتریس { \bar { P } } ، متقارن است. چرا که رابطه ی { \bar { P } }={ \bar { P } }^{ T } برقرار است. اکنون از قضیه ای به نام فاکتورگیری تاکاگی ( که به آن تجزیه به مقادیر تکین متقارن یا SSVD نیز می گویند ) استفاده می کنیم [13] . این قضیه به ما می گوید که اگر ماتریسی ( مانند {\bar{ P }} )مربعی و متقارن باشد، می توان آن را به صورت زیر تجزیه کرد:
که \Lambda ماتریسی است که مقادیر آن، جذر مقادیر تکین ماتریس {\bar{ P }} است. به طور کامل تر، می دانیم که طبق تجزیه ی SVD می‌توان ماتریس {\bar{ P }} را به صورت \bar { P } =U\Lambda { V }^{ H } تجزیه کرد. قابل ذکر است که ماتریس {\bar{ P }} می تواند مربعی نباشد، U و V ماتریس های واحد و \Lambda هم ماتریسی با ابعاد {\bar{ P }} است که عناصر آن جذر مقادیر تکین {\bar{ P }}هستند. اگر {\bar{ P }} ماتریسی مربعی باشد ( که هست! )، \Lambda ماتریسی قطری می شود.
همان طور که بیان شد، داریم:
ابتدای دو طرف رابطه را { U }^{ H } در ضرب می کنیم :
در رابطه ی بالا از خاصیت واحد بودن ماتریس U استفاده کردیم. حال انتهای دو طرف رابطه را در { U }^{ \ast } ضرب می کنیم :
حال با تعریف D={ U }^{ H }{ V }^{ \ast } ، به رابطه ی زیر می رسیم:
که D یک ماتریس واحد است. چرا که از ضرب دو ماتریس واحد حاصل شده است.
حال به رابطه ی (28) باز می گردیم:
اگر Y=U{ D }^{ \frac { 1 }{ 2 } } را تعریف کنیم، به عبارت زیر می رسیم:
به رابطه ی (30) تبدیل تاکاگی یا SSVD می گویند. باز هم مشخص است که Y یک ماتریس واحد است. چرا که از ضرب دو ماتریس واحد به دست می آید.
همان طور که پیش تر گفته شد، به دنبال ماتریس Q می گردیم که در رابطه ی (23) صدق کند. اگر Q={ Y }^{ H }G را تعریف کنیم، آنگاه:
پس ملاحظه می شود که با تعریف Q={ Y }^{ H }G ، ماتریس مورد نظر در رابطه ی (23) را یافتیم. برای کواریانس می توان نوشت:
اگر { S }_{ a }=Q{ C }_{ a }{ Q }^{ H } و { S }_{ b }=Q{ C }_{ b }{ Q }^{ H } را تعریف کنیم، مانند الگوریتم اصلی CSP ، می توان گفت که ماتریس های { S }_{ a } و { S }_{ a } ماتریس بردار های ویژه ی مشترکی دارند. به طوری که:
به همین ترتیب، می خواهیم برای شبه کواریانس دو کلاس هم، ماتریس بردار های ویژه ی مشترکی تخمین بزنیم. پس:
ابتدا و انتهای عبارت بالا را در { \Lambda }^{ -\frac { 1 }{ 2 } } و { \Lambda }^{ -\frac { 1 }{ 2 } T } ضرب می کنیم:
در رابطه ی بالا، \Lambda یک ماتریس قطری با عناصر حقیقی است. پس رابطه ی { \Lambda }^{ -\frac { 1 }{ 2 } }={ \Lambda }^{ -\frac { 1 }{ 2 } T } برقرار است. برای یافتن ماتریس بردار های ویژه ی مشترک شبه کواریانس ها، به صورت زیر عمل می کنیم:
در حقیقت با اعمال تبدیل SUT ، به طور هم زمان ماتریس های کواریانس و شبه کواریانس را به ماتریس همانی تبدیل کردیم ( رابطه ی (32 )و (34)). با این کار، اختلاف واریانس بین دو کلاس ماکزیمم خواهد شد. این بار به جای یک دسته فیلتر فضایی، برای هر کدام از کواریانس و شبه کواریانس، یک دسته فیلتر فضایی تعریف می کنیم:
و به همین ترتیب، سیگنال ها به صورت زیر فیلتر می شوند:
تا بتوانیم المان های بردار های ویژگی را طبق (17) و (18) یا (13) به دست آوریم.
زمانی که ماتریس کواریانس مرکب را با اعمال تبدیل سفید کنندگی، به ماتریس همانی تبدیل می کنیم، مجموع ماتریس بردار های ویژه کلاس a و b همانی می شود. این کار باعث می شود که سطر اول ماتریس فیلتر شده ی V ، بیشترین واریانس را در کلاس a و کمترین واریانس را در کلاسbداشته باشد و همین طور سطر آخر این ماتریس فیلتر شده بیشترین واریانس را در کلاس a و کمترین واریانس را در کلاس b دارد. حال در این جا با اعمال تبدیل SUT به غیر از ماتریس کواریانس مرکب، شبه کواریانس مرکب نیز به ماتریسی همانی تبدیل می شود. پس طبیعی است که سیگنال های فیلتر شده به وسیله ی فیلتر های فضایی شبه کواریانس نیز، خاصیت های V را در مورد اختلاف واریانس ها در کلاس های مختلف داشته باشند.
می دانیم که تبدیل تاکاگی منجر به رابطه ی (23) شد. نگاهی به این رابطه می اندازیم:
می دانیم کهQZ مقادیری مختلط دارد. پس مقادیر ماتریس \Lambda به صورت زیر بیان می شوند:
این در حالیست که ماتریس \Lambdaای که از تبدیل تاکاگی به دست می آید، مقادیری حقیقی دارد. این موضوع نشان می دهد که مقادیر ماتریس \Lambda که به صورت \lambda =E\left[ { z }_{ r }^{ 2 } \right] -E\left[ { z }_{ i }^{ 2 } \right] هستند، اطلاعات اضافی راجع به اختلاف توان بین قسمت های حقیقی و موهومی منبع ورودی به ما می دهند. پس روش SUTCCSP بر خلاف روش CSP ، تفاوت بین توان قسمت های حقیقی و موهومی ورودی را نیز در نظر می گیرد.
۷. رابطه ی بین CSP و ACCSP
در این قسمت نشان می دهیم که بین روش ACCSP و CSP ، دوگانگی وجود دارد. فرض کنید سیگنال تصادفی مختلطی با میانگین صفر مربوط به کلاس a داریم. z\left( t \right) ={ z }_{ r }\left( t \right) +j{ z }_{ i }\left( t \right) که در آن { z }_{ r }\left( t \right) و { z }_{ i }\left( t \right)به ترتیب، قسمت حقیقی و موهومی این سیگنال هستند و هر دو، سیگنال هایی تصادفی می باشند. فُرم افزوده ی این سیگنال مختلط را به صورت \hat { { Z }_{ a } } =\begin{bmatrix} z\left( t \right) \\ { z }^{ * }\left( t \right) \end{bmatrix} تعریف می کنیم. این ماتریس را می توان به صورت زیر نیز نمایش داد:
اگر ماتریس های زیر را تعریف کنیم:
با استفاده از این تعریف ها، ماتریس کواریانس { \hat { Z } }_{ a } به صورت به دست می آید:
در این رابطه، ماتریس Z مقادیری حقیقی دارد. پس { Z }^{ H }={ Z }^{ T } می باشد.
پس از این که تبدیل سفید کنندگی را به ماتریس کواریانس { \check { Z } }_{ a } (رابطه‌ی(39)) اعمال کنیم، این ماتریس، یک ماتریس همانی می‌شود. به عبارت دیگر E\left[ { \check { Z } }_{ a }{ \check { Z } }_{ a }^{ H } \right] =I می شود. با این کار رابط ی (39 ) به صورت I\quad =\quad \Phi E\left[ \check { Z } { \check { Z } }^{ T } \right] { \Phi }^{ H } تبدیل می شود. پس می توانیم بگوییم:
(40) { \Phi }^{ -1 }{ \left( { \Phi }^{ H } \right) }^{ -1 }=E\left[ \check { Z } { \check { Z } }^{ T } \right]
(41) { \Phi }^{ -1 }{ \left( { \Phi }^{ H } \right) }^{ -1 }=\begin{bmatrix} \frac { 1 }{ 2 } & \frac { 1 }{ 2 } \\ -j\frac { 1 }{ 2 } & j\frac { 1 }{ 2 } \end{bmatrix}\times \begin{bmatrix} \frac { 1 }{ 2 } & j\frac { 1 }{ 2 } \\ \frac { 1 }{ 2 } & -j\frac { 1 }{ 2 } \end{bmatrix}=\begin{bmatrix} \frac { 1 }{ 2 } & 0 \\ 0 & \frac { 1 }{ 2 } \end{bmatrix}=\frac { 1 }{ 2 } I\quad \quad \Longrightarrow \quad \quad E\left[ \check { Z } { \check { Z } }^{ T } \right] =\frac { 1 }{ 2 } I
رابطه ی (41) به ما می گوید که در روش ACCSP ، زمانی که کواریانس یک ماتریس کواریانس \hat { { Z }_{ a } } را قطری می کنیم، ماتریس E\left[ \check { Z } { \check { Z } }^{ T } \right] (با مقادیری حقیقی) نیز قطری می شود. این موضوع نشان می دهد که عملکرد ACCSP مانند CSP است. با این تفاوت که در ACCSP پس از اعمال تبدیل سفید کنندگی، قطری شدن ماتریس کواریانس مقادیر حقیقی، با ضریب 0.5 همراه است و همانی نمی شود!
۸. آزمایشها
داده ها
در این پروژه، از دو نوع داده به عنوان ورودی برای فیلتر های فضایی و تعمیم های آن استفاده شده است. در این بخش نحوه ی ایجاد داده های ساختگی و نحوه‌ی استفاده از داده های واقعی را شرح می دهیم.
الف- داده های ساختگی
برای هر یک از کلاس ها، تعدادی مجموعه داده ی ساختگی می سازیم. این داده ها مجموعی از سینوس ها هستند که دامنه ی فرکانسی در رنج دامنه ی فرکانسی امواج آلفا و بتای سیگنال EEG دارند ( 8Hz - 30Hz ). این سینوسی ها به نویز سفید تصادفی آلوده شده اند. این سیگنال ها را به صورت زیر تعریف می‌کنیم:
1) کلاس a: تعداد 100 دسته داده را به صورت زیر تعریف می کنیم:
2) کلاس b: تعداد 100 دسته داده را به صورت زیر تعریف می کنیم:
که در آن { f }_{ 3 }=9㎐، { f }_{ 2 }=19㎐، { f }_{ 1 }=10㎐ و { f }_{ 4 }=17㎐ هستند.
و در ضمن { v }_{ 1 }\left( t \right)، ... ، { v }_{ 1 }\left( t \right)، { v }_{ 1 }\left( t \right) نویز های سفید متفاوتی هستند به طوری که نسبت سیگنال به نویز (SNR) همه ی کانال ها بین 9dB- تا 15dB- تغییر می کند.
در نرم افزار متلب، فرکانس نمونه برداری را 100Hz انتخاب می کنیم و داده ها را به مدت 100 ثانیه شبیه سازی می کنیم. تابعی تعریف می کنیم که به ازای ورودی یک سیگنال و یک عدد، خروجی آن جمع آن سیگنال با نویز سفید باشد به طوری که SNR خروجی، برابر عدد وارد شده در ورودی باشد. تابع add_wgn را به صورت زیر در متلب تعریف می کنیم:
function [X] = add_wgn(original_signal,final_SNR)
Pn=10*log10(var(original_signal))-final_SNR;
X=original_signal+wgn(1,length(original_signal),Pn);
end
در این تابع، originalsignal سیگنال اصلی ما و final_SNR ، SNR سیگنال نهایی است. این تابع توان سیگنال ورودی را بر حسب دسیبل به دست آورده و طبق رابطه ی $SNR\quad =\quad { P }{ signal }-{ P }_{ noise }$ ( بر حسب دسیبل )، توان نویز را به دست می آورد. هم چنین از تابع wgn که از توابع نرم افزار متلب می باشد نیز استفاده شده است:
y=wgn(m,n,p)
خروجی این تابع، ماتریسی m\times n از نویز سفید گوسی با توان p است.
همچنین تابع دیگری با عنوان randnum تعریف می شود. هدف از این تابع ایجاد یک عدد تصادفی بین بازه ی وارد شده در ورودی آن است تا بتوان از آن یک عدد تصادفی بین بازه ی 9dB- تا 15dB- به عنوان SNR ایجاد کرد. تابع randnum را به صورت زیر تعریف می کنیم:
function [y] = randnum(a,b)
y = a + (b-a).*rand(1,1);
end
که a و b ، ورودی ها به عنوان بازه است و y عدد تصادفی در بازه ی ورودی می باشد. همچنین از تابعی از توابع متلب به نام rand نیز استفاده کرده ایم:
r = rand(m,n)
این تابع، ماتریسی m\times n از اعداد طبیعی در بازه ی [0,1 ] به عنوان خروجی ایجاد می کند.
حال برای شبیه سازی مجموعه داده های کلاس a و b ، از کد زیر استفاده می کنیم:
f1=10;f2=19;f3=10.1;f4=18.9;Fs=100;
t=0:(1/Fs):100;
for i=1:100 class_a{i}=[add_wgn(1.00*sin(2*pi*f1*t+pi/8)+1.05*sin(2*pi*f2*t),randnum(-9,-15)); add_wgn(1.11*sin(2*pi*f1*t)+1.15*sin(2*pi*f2*t),randnum(-9,-15)); add_wgn(1.95*sin(2*pi*f1*t+pi/12)+0.05*sin(2*pi*f2*t),randnum(-9,-15));
add_wgn(2.13*sin(2*pi*f1*t)+1.03*sin(2*pi*f2*t-pi/3),randnum(-9,-15))];
class_b{i}=[add_wgn(1.18*sin(2*pi*f3*t)+1.17*sin(2*pi*f4*t),randnum(-9,-15)); add_wgn(1.02*sin(2*pi*f3*t)+1.04*sin(2*pi*f4*t),randnum(-9,-15));
add_wgn(1.45*sin(2*pi*f3*t)+1.23*sin(2*pi*f4*t-pi/5),randnum(-9,-15)); add_wgn(0.98*sin(2*pi*f3*t)+1.14*sin(2*pi*f4*t),randnum(-9,-15))];
end
save('Synthetic_data.mat','class_a','class_b','t');
این m-file ، به ازای هر کلاس، یک سلول برای ما ایجاد می کند. هر کدام از این سلول ها 100 عضو دارند (که مسلماً می توانند تعداد مختلف و نابرابری داشته باشتند). هر عضو سلول ها، ماتریسی 4\times T است. یعنی هر داده 4 کانال دارد و هر کدام از این کانال ها، T نمونه دارند. با توجه به تعریف فرکانس نمونه برداری و مدت زمان شبیه سازی، T=1001 است.
ب- داده های واقعی
برای داده های واقعی از مجموعه داده های اول رقابت BCI چهارم استفاده می کنیم. داده ها از 7 نفر داوطلب سالم به وسیله ی BrainAmp MR و تقویت کننده ی EEG گرفته شده است. این ثبت 59 کاناله بوده و فرکانس نمونه برداری در آن 1000Hz بوده است. البته فرکانس نمونه برداری داده هایی که در این پروژه استفاده شده است، به 100Hz کاهش یافته است.
آزمایش انجام شده از این قرار است: از هر یک از داوطلبان خواسته می شود که به دلخواه خود، دو تا از سه تصور حرکتی را انتخاب کرده و در هر بار آزمایش یکی از آن ها را انجام دهدند: تصور حرکت دست چپ، دست راست و پاها. برای هر کدام از داوطلبان این آزمایش 200 بار انجام می شود. در هر بار آزمایش داوطلب به صورتی که مانیتور رو به روی او نشان می دهد، به مدت 4 ثانیه یکی از دو تصور را انجام می هد. انجام این آزمایشات طوری انجام می شود که هر داوطلب، هر کدام از تصورات را 100 بار انجام دهد. این در حالیست که این تصورات پشت سر هم اتفاق نمی افتد و به نحوی تصادفی انتخاب می شوند [12] .20
به دلیل مفاهیم نئوروفیزیولوژیکی ، به جای استفاده از 59 کانال، از 11 کانال برای انجام آنالیز ها استفاده می کنیم. کانال های انتخابی “FC3” ، “FC4” ، “Cz” ، “C3” ، “C4” ، “C5” ، “C6” ، “T7” ، “T8” ، “CCP3” و “CCP4” بر اساس سیستم 10-20 هستند. چرا که در [13] و [14] آمده که تصورات حرکتی در درجه ی اول با مناطق مرکزی مغز ارتباط دارند!
سیگنال ها در این آزمایشات پیوسته هستند. برای هر داوطلب یک 59 کانال سیگنال وجود دارد. متغیری به نام mrk ، زمان و تصور انجام شده هر آزمایش را مشخص می کند. داده ها برای هر داوطلب در متلب دارای متغیر های زیر هستند:
الف- cnt: ماتریسی شامل 59 کانال سیگنال EEG پیوسته T نمونه است. این داده ها با فرمت INT16 هستند.
ب- mrk: یک ساختار است که حاوی اطلاعاتی راجع به مکان اتقاق افتادن کلاس ها می باشد. این ساختار شامل رشته های زیر است:
ا.pos: برداری شامل مکان انجام هر یک از کلاس ها در سیگنال EEG داوطلب است. این مکان را در واحد نمونه به ما می دهد. یعنی در کدام نمونه از سیگنال EEG کلاس انجام شده است.
ا.y: برداری شامل برچسب هر یک از کلاس ها در سیگنال EEG داوطلب است. مقادیر آن 1- برای کلاس اول و 1 برای کلاس دوم است.
ج- nfo: یک ساختار است که اطلاعات بیشتری را به صورت رشته های زیر به ما می دهد:
ا- { f }_{ s }فرکانس نمونه برداری است.
ا- clab سلولی شامل برچسب کانال ها است.
ا- classes سلولی راجع به نام هر یک از کلاس های مربوط به تصورات حرکتی.
ا- xpos موقعیت بردار x هر الکترود در صفحه ی دو بعدی (بر حسب سیستم 20-10).
ا- ypos موقعیت بردار y الکترود در صفحه ی دو بعدی (بر حسب سیستم 20-10).
می خواهیم تابعی به نام real_data تعریف کنیم که با گرفتن داده های رقابت BCI به عنوان ورودی، دو خروجی به ما بدهد. مانند داده های ساختگی می خواهیم به ازای هر کلاس، یک سلول داشته باشیم. پس هر سلول حاوی 100 ماتریس می شود. چرا که آزمایشات به نحوی تنظیم شده اند که از هر کلاس 100 داده داشته باشیم. هر کدام از این ماتریس ها یک ماتریس N\times T است که تعداد کانال ها در این جا 11 تا و T تعداد نمونه هایی است که داوطلب حرکت مربوطه را تصور کرده است. در این تابع 11 کانال (که پیش تر گفته شدند) از 59 کانال را انتخاب می کنیم. این انتخاب بر اساس اطلاعاتی است که در nfo.clab راجع به برچسب کانال ها وجود دارد. همچنین برای تبدیل مقادیر cnt به mV از دستور زیر استفاده می کنیم:
cnt= 0.001*double(cnt);
همچنین می دانیم که تصورات حرکتی در سیگنال EEG، بازه ی فرکانسی در حدود 8Hz-30Hz دارند [1] و [2] . پس با استفاده از یک فیلتر میان گذر باترورث درجه ی 6 ، سیگنال ها را در بازه ی 8Hz-30Hz فیلتر می کنیم. برای این کار از دو تابع butter و filter در متلب استفاده می کنیم:
[b,a]=butter(n,Wn,'ftype')
y = filter(b,a,X)
در تابع butter ، داریم n به عنوان درجه ی فیلتر، { W }_{ n } به عنوان فرکانس قطع ( بازه ی فرکانس قطع) و 'ftype' به عنوان نوع فیلتر، ورودی های این تابع هستند. پس ما باید بازه ی فرکانس قطع را به صورت { W }_{ n }=\left[ \frac { 8 }{ \sfrac { { f }_{ s } }{ 2 } } ,\frac { 30 }{ \sfrac { { f }_{ s } }{ 2 } } \right] انتخاب کنیم. چرا که می خواهیم از فیلتر میان گذر استفاده کنیم، پس دو فرکانس قطع می خواهیم. همچنین برای استفاده از تابع butter ، فرکانس (های) قطع باید عددی (اعدادی) بین 0 و 1 باشند (نرمالیزه شده). از آنجا که طبق قضیه ی نایکویست می دانیم که بزرگترین فرکانس موجود دا سیگنال های گسسته زمان، نصف فرکانس نمونه برداری شده است [15] ، پس فرکانس های قطع را بر \sfrac { { f }_{ s } }{ 2 } تقسیم می‌کنیم. همچنین برای فیلتر میان گذر مقدار 'ftype' = 'bandpass' را تعیین می کنیم.
تابع filter در متلب، ضرایب فیلتر ( b به عنوان ضرایب صورت و a به عنوان ضرایب مخرج فیلتر) و سیگنال اصلی را به عنوان ورودی گرفته و سیگنال فیلتر شده را به ما تحویل می دهد.
پس تابع real_data را به صورت زیر تعریف می کنیم:
function [class_a,class_b] = real_data(x) cnt=x.cnt; mrk=x.mrk; nfo=x.nfo;
cnt=double(cnt)/1000;
cnt2(:,1)=(cnt(:,27));
cnt2(:,2)=(cnt(:,31));
cnt2(:,3)=cnt(:,11);
cnt2(:,4)=cnt(:,15);
cnt2(:,5)=cnt(:,29);
cnt2(:,6)=cnt(:,26);
cnt2(:,7)=cnt(:,32);
cnt2(:,8)=cnt(:,25);
cnt2(:,9)=cnt(:,33);
cnt2(:,10)=cnt(:,36);
cnt2(:,11)=cnt(:,39);
[b,a]=butter(3,[8/(nfo.fs/2),30/(nfo.fs/2)],'bandpass');
i_a=1;
i_b=1;
for i=1:length(mrk.y)
n=mrk.pos(i)+nfo.fs*4;
clear cnt21;
for j=1:size(cnt2,2)
end
if mrk.y(i)==-1
class_a{i_a}=cnt21';
i_a=i_a+1;
elseif mrk.y(i)==1
class_b{i_b}=cnt21';
i_b=i_b+1;
end
end
ورودی این تابع تعریف شده، ساختاری است که رشته های آن، mrk، cnt و nfo هستند. همان طور که گفته شد، mrk.yحاوی برچسب کلاس هر آزمایش است. همچنین می دانیم که هر داوطلب، در هر آزمایش 4 ثانیه تصور حرکت مورد نظر ما را انجام داده است. با یک حساب ساده ی ریاضی، در می یابیم که انجام هر آزمایش، 400 نمونه طول می کشد. چرا که در هر ثانیه، 100 نمونه برداشته می شود، پس مشخص است که در 4 ثانیه تعداد 400 نمونه برداشته می شوند. پس می دانیم که مکان هر آزمایش در سیگنال EEG ،cnt\left( mrk.pos\left( i \right) :n \right) است. mrk.pos\left( i \right)مکان آزمایش iاُم در سیگنال هاست . یعنی از زمان شروع آزمایش، تا n نمونه بعد از آن! لازم به ذکر است که cnt2 ، کانال های انتخابی از کل کانال ها بر حسب mV است.
می دانیم که mrk.y حاوی برچسب کلاس هر آزمایش است. پس، به ازای هر آزمایش، چک می کنیم که آیا آزمایش انجام شده به کلاس اول 1- و یا کلاس دوم 1 تعلق دارد. پس بازه ی انتخاب شده از سیگنال EEG در هر آزمایش را در یکی از دو سلول class_a و class_b می ریزیم. لازم به ذکر است که برای کمتر شدن حجم محاسبات، به جای فیلتر کردن کل سیگنال، فقط بازه ی انتخابی در هر بار آزمایش را فیلتر می کنیم.
پس خروجی تابع real_data به ازای ورودی مناسب، دو کلاس class_a و class_b است. در حقیقت این تابع، با گرفتن کل سیگنال و اطلاعات آن به عنوان ورودی، قسمت هایی از سیگنال که تصور مربوط به کلاس اول و یا دوم انجام شده را در سلول مربوط به آن کلاس می ریزد. چرا که هدف ما در این پروژه پردازش و کلاس بندی داده ها به صورت آنلاین نیست!
هدف، بررسی عملکرد روش های به کار گرفته شده است.
۹. شبیه سازی فیلتر ها فضایی
در این قست نحوه‌ی شبیه‌سازی و فیلترینگ، انتخاب ویژگی و کلاس‌بندی را شرح می‌دهیم.
الف- فیلترهای فضایی مشترک (CSP)
روش اصلی و اولیه‌ی فیلترهای فضایی CSP است. الگوریتم این روش به صورت زیر است [1],[2]:
1) ماتریس داده‌های میانگین صفر { Z }_{ a }\in { R }^{ N\times T } و { Z }_{ b } را تشکیل می‌دهیم. N تعداد کانال‌ها و T تعداد نمونه‌های هر سیگنال است. منظور از R، مجموعه‌ی اعداد حقیقی است.
2) از رابطه‌ی { C }_{ a }=\frac { { z }_{ a }{ z }_{ a }^{ H } }{ tr\left( { z }_{ a }{ z }_{ a }^{ H } \right) } و { C }_{ b }=\frac { { z }_{ b }{ z }_{ b }^{ H } }{ tr\left( { z }_{ b }{ z }_{ b }^{ H } \right) }، ماتریس کواریانس هر کلاس برای تمامی داده‌ها را به دست آورده و سپس میانگین آنها یعنی \overline { { C }_{ a } } و \overline { { C }_{ b } } را می‌یابیم.
3) ماتریس کواریانس مرکب را به صورت { C }_{ c }=\overline { { C }_{ a } } +\overline { { C }_{ b } } تعریف می‌کنیم.
4) ماتریس بردارهای ویژه { U }_{ c } و ماتریس مقادیر ویژه { \Lambda }_{ c } مربوط به ماتریس کواریانس مرکب را به دست می‌آوریم.
5) تبدیل سفید کنندگی را به صورت G=\sqrt { { \Lambda }_{ c }^{ -1 } } { U }_{ c }^{ H } تعریف می‌کنیم.
6) حال { S }_{ a }=G\overline { { C }_{ a } } { G }^{ H } و { S }_{ b }=G\overline { { C }_{ b } } { G }^{ H } را تعریف می‌کنیم.
7) ماتریس مقادیر ویژه { \Lambda }_{ a } و بردار های ویژه B را برای { S }_{ a } می‌یابیم.
8) و { \Lambda }_{ a } (و متناسب با آن B) را به صورت نزولی مرتب می‌کنیم. می‌دانیم که ماتریس ویژه‌ی { S }_{ b } نیز هست.
9) فیلتر های فضایی را به صورت { W }^{ H }={ B }^{ H }G تعریف می‌کنیم.
10) در { W }^{ H }، m سطر اول و m سطر آخر را انتخاب می کنیم.
11) برای هر داده ی جدید Z ، داده های فیلتر شده را به صورت V={ W }^{ H }Z می یابیم. ماتریس جدید V یک ماتریس 2m\times T است.
12) با استفاده از رابطه ی { f }_{ p }=\log { \left( \frac { var\left( { v }_{ p } \right) }{ \sum _{ i=1,...,m\quad and\quad N-m+1,...,N }^{ }{ var\left( { v }_{ i } \right) } } \right) }، به ازای هر سطر V، یک ویژگی تعریف می کنیم.{ v }_{ p } سطر pاُم ماتریس V است.
می دانیم که آزمایشات متعلق به کلاس a و b ، دو سلول هستند. چه از داده‌های ساختگی استفاده کنیم و چه از داده‌های واقعی، در قسمت قبل، طوری داده‌ها را ساختیم که هر کلاس، یک سلول داشته باشد، تعداد اعضای این سلول به تعداد آزمایشات مربوط به آن کلاس بوده و هر عضو آن، ماتریسی N\times T است که در آن N تعداد کانال ها و T تعداد نمونه هایی است که طول کشیده تا داوطلب حرکت مربوطه(کلاس مربوطه) را تصور کند.
پیش از شبیه سازی روش CSP ، می بایست داده های هر کلاس را به دو مجموعه داده ی آموزشی و آزمایشی تقسیم کنیم. چرا که در یک مسئله ی شناسایی آماری الگو، می بایست کلاسیفایر را به وسیله ی داده های آموزشی، آموزش داد و سپس با استفاده از داده های آزمایشی، عملکرد را سنجید. پس در این جا برای به دست آوردن فیلتر های فضایی (روش CSP) و آموزش کلاسیفایر از داده های آموزشی استفاده می‌کنیم. سپس داده‌های آزمایشی را به وسیله ی این فیلتر های به دست آمده فیلتر کرده و آن ها را با کلاسیفایر آموزش داده شده، ارزیابی می کنیم.
یکی از نکات مهم این است که در داده هایی که در به دست آوردن فیلتر ها (و آموزش کلاسیفایر) استفاده شده اند، نباید در ارزیابی روش، مورد استفاده قرار گیرند. چرا که بحث Overfitting مطرح می شود [16] . ارزیابی به وسیله ی داده هایی که در آموزش روش استفاده شده اند، عملکرد حداکثر را به همراه دارد. ولی اگر این فیلتر ها و کلاسیفایر روی داده های جدیدی آزمایش شود، خطا و کلاس بندی اشتباه زیادی مرتکب می شوند! پس داده های آموزشی و آزمایشی باید از هم مستقل باشند. بحث دیگری که در این کار مطرح است، این است که چگونه داده های آموزشی و آزمایشی را انتخاب کنیم؟ برای مثال اگر بخواهیم داده ها را به 80% داده ی آموزشی و 20 % داده ی آزمایشی تقسیم کنیم، چگونه این تقسیم بندی اتفاق بیافتد. تقسیم‌بندی و شکستن داده ها به طور مختلف، عملکردهای متفاوتی را به همراه است. پس در این پروژه، شکستن داده ها به داده های آموزشی و آزمایشی (با نسبت خاصی) به صورت تصادفی انجام می شود. سپس برای تایید این نتایج، آزمایش به تعداد مشخصی تکرار و هر بار داده ها به صورت تصادفی (با نسبت خاصی) به داده های آموزشی و آزمایشی تقسیم می شود.
فرض کنیم که class_a و class_b ، سلول های مربوط به داده های کلاس اول و کلاس دوم و متغیر perc نسبت تقسیم داده ها باشند. می خواهیم به صورت تصادفی و با نسبت ، داده های آموزشی و آزمایشی را برای هر کلاس انتخاب کنیم. تابعی به نام train_test تعریف می کنیم که به ازای یک بردار با طول مشخص و دو عدد به عنوان ورودی، بردار را به دو بردار با طول هایی برابر دو عدد ورودی به طور تصادفی تقسیم کند. به عنوان مثال اگر برداری با طول 10 و دو عدد 8 و 2
در ورودی تابع وارد شد، خروجی تابع (به طور تصادفی) برداری با طول 8 و بردار دیگری با طول 2 باشد. پس تابع train_test را به صورت زیر تعریف می کنیم:
function [ training,test ] = train_test( a,n_training,n_test )
n=length(a);
if ((n_training>0) && (n_test >0) && (n_training+n_test==n))
training=[];
test=[];
for i=1:n_training
j=randi(n-i+1);
training(i)=a(j);
a(j)=[];
end
test=a;
else
error('Inputs didnt enter correctly!')
end
end
در این تابع بردار مورد نظر، n_training تعداد داده ها در داده های آموزشی و n_test تعداد داده ها در داده های آزمایشی است. این تابع چک می کند که آیا هر کدام از تعداد داده های آموزشی و آزمایشی ورودی، اعدادی مثبت هستند و این که مجموع این اعداد با مجموع داده های بردار برابر هستند یا خیر. اگر این شرط ها برقرار باشند، دو بردار به نام training و test تعریف می شوند. سپس به ازای i از 1 تا تعداد داده های آزمایشی، هر بار یک عدد صحیح j بین 1 تا n-i+1 ایجاد می شود، سپس عضو jاُم بردار a را داخل بردار training می ریزیم و سپس عضو jاُم بردار a را حذف می کنیم. این کار باعث می شود که به صورت تصادفی، تعداد n_training از بردار a را انتخاب شده و این مقادیر انتخاب شده از a حذف شوند. در انتها مقادیر باقی مانده در بردار a را در بردار test می ریزیم. در نتیجه تابع train_test ، بردار را به دو بردار مستقل با طول دلخواه تبدیل می کند.
لازم به ذکر است که تابع randi در متلب، عدد صحیحی بین 1 تا ورودی تابع برای ما ایجاد می کند. یعنی:
r = randi(imax,n);
همان طور که گفته شد این تابع بردار n\times n شامل اعداد 1 تا imax به صورت تصادفی برای ما ایجاد می کند. اگر n را به این تابع ندهیم، خروجی یک عدد اسکالر می شود!
حال برای تقسیم داده ها در هر کلاس به دو دسته داده ی آموزشی و آزمایشی از کُد زیر استفاده می کنیم:
n_a=size(class_a,2);
n_b=size(class_b,2);
n_train_a=fix(n_a*perc);
n_test_a=n_a-n_train_a;
n_train_b=fix(n_b*perc);
n_test_b=n_b-n_train_b;
[train_a,test_a]=train_test(1:n_a,n_train_a,n_test_a);
[train_b,test_b]=train_test(1:n_b,n_train_b,n_test_b);
class_a_train={};
class_a_test={};
class_b_train={};
class_b_test={};
N=size(class_a{train_a(1)},1);
%class A
%training
for i=1:n_train_a
class_a_train{i}=class_a{train_a(i)};
if size (class_a_train{i},1)~=N
error('All the trials must have a equal number of channels' )
end
end
%test
for i=1:n_test_a
class_a_test{i}=class_a{test_a(i)};
if size (class_a_test{i},1)~=N
error('All the trials must have a equal number of channels' )
end
end
%class B
%training
for i=1:n_train_b
class_b_train{i}=class_b{train_b(i)};
if size (class_b_train{i},1)~=N
error('All the trials must have a equal number of channels' )
end
end
%test
for i=1:n_test_b class_b_test{i}=class_b{test_b(i)};
if size (class_b_test{i},1)~=N
error('All the trials must have a equal number of channels' )
end
end
در na و n_b تعداد آزمایشات متعلق در کلاس a و b برای داوطلب مورد نظر می باشند. چرا که هر کدام از class_a و class_b، یک سلول با ابعاد 1\times n هستند و برای هر سلول، تعداد آزمایشات (داده ها) متعلق به کلاس مربوطه است. سپس متغیر را به عنوان تعداد داده های آموزشی کلاس تعریف می کنیم. همان طور که پیش تر گفته شد perc نسبت تقسیم داده ها به داده های آموزشی و آزمایشی است و تابع fix که از توابع متلب است،اعداد را گرد می کند! سپس n_test_a=n_a-n_train_a را تعداد داده های آزمایشی کلاس a تعریف می کنیم. همین کار برای کلاس b نیز انجام می شود. سپس با استفاده از تعداد داده های آموزشی، تعداد داده های آزمایشی هر کلاس و برداری شامل اعداد 1 تا تعداد داده های کل یک کلاس و تابع train_test ، دو بردار train_a و test_a برای کلاس a و دو بردار train_b و test_b برای کلاس ایجاد می کنیم.
سپس چهار سلول class_a train،class_ test، class_b_train و class_b_test را تعریف کرده و در هر کلاس برای هر کدام از داده های آموزشی و آزمایشی، اعضای موجود در بردار های test_a، train_a،test_b و train_b را از class_a و class_b انتخاب کرده و داخل سلول های مربوطه می ریزیم. این کار باعث می شود که داده های کلاس a و کلاس b، به دو دسته داده ی آموزشی و آزمایشی برای هر کلاس به طور مستقل تقسیم شوند. یعنی در هر کلاس، هیچ یک از داده ها در داده های آموزشی و آزمایشی مشترک نیست! یادآوری می شود که هر کدام از اعضای سلول ها (برای هر کلاس) ماتریسی V است که N تعداد کانال ها و T تعداد نمونه های آن داده ی خاص است.
حال دو سلول برای داده های آموزشی و آزمایشی کلاس اول و دو سلول برای داده های آموزشی و آزمایشی کلاس دوم ایجاد کرده ایم. اکنون باید ماتریس کواریانس هر یک از داده های آموزشی (در هر کلاس) را به دست آورده و میانگین آن ها را حساب کنیم. با توجه به خلاصه ی الگوریتم CSP که در ابتدا ارائه شد، از کُد زیر استفاده می کنیم:
%Class A
Ca=0;
Z=0;
for i=1:n_train_a
Z=class_a_train{i};
mn=mean(Z')';
for j=1:size(Z,2)
Z(:,j)=Z(:,j)-mn
end Ca=Ca+(Z*Z')/trace(Z*Z');
end Ca=Ca/n_train_a;
%class B
Cb=0;
Z=0;
for i=1:n_train_b
Z=class_b_train{i};
mn=mean(Z')';
for j=1:size(Z,2)
Z(:,j)=Z(:,j)-mn;
end Cb=Cb+(Z*Z')/trace(Z*Z');
end
Cb=Cb/n_train_b;
در این کُد برای هر کلاس، ابتدا هر یک از اعضای سلول داده های آزمایشی انتخاب می شده و در Z ریخته می شوند. mnکه برداری N\times 1 است و هر کدام از سطر های آن، میانگین کل نمونه های کانال مربوطه است، حساب می شود (واضح است که N تعداد کانال ها است، پس بردار mn باید N سطر داشته باشد!). سپس برای میانگین صفر شدن، هر یک از داده ها از میانگین خود (mn) کم می شود. همان طور که در [1] و [2] آمده است، داده ها باید میانگین صفر داشته باشند تا ماتریس کواریانس آن ها با ماتریس همبستگی آن ها برابر شود. متغیر { C }_{ a } و { C }_{ b } که به عنوان ماتریس‌های کواریانس هر کلاس در نظر گرفته شده، در ابتدا صفر فرض می شوند. برای هر داده از طریق رابطه ی { \left( Z\times \acute { Z } \right) }/trace{ \left( Z\times \acute { Z } \right) } ، ماتریس کواریانس حساب شده و با مقدار قبلی جمع می شود. trace که از توابع متلب است، مجموع عناصر روی قطر اصلی ماتریس ورودی خود را حساب می کند. همچنین mean ( از توابع متلب ) ، به صورت زیر استفاده می شود:
M = mean(A)
اگر ماتریس A ، ماتریسی m\times n باشد، برداری سطری شامل میانگین هر یک از ستون های ماتریس A است. برای این که میانگین سطر های ماتریس را به صورت برداری ستونی داشته باشیم از کُد mn=mean{ \left( { Z }^{ \prime } \right) }^{ \prime } استفاده کرده ایم. لازم به ذکر است که ، عملگر ترانهاده ی مزدوج مختلط در متلب است(بدیهی است که اگر این عملگر روی داده های حقیقی اعمال شود، معادل با ترانهاده است! ).
در انتهای این کد، با تقسیم هر یک از { C }_{ a } و { C }_{ a } بر تعداد داده های آموزشی هر کلاس، ماتریس های کواریانس میانگین ها هر کلاس به دست می آیند.
ماتریس کواریانس مرکب را به صورت زیر تعریف می کنیم:
حال می خواهیم ماتریس بردار های ویژه و ماتریس مقادیر ویژه ی ماتریس کواریانس مرکب را به دست آوریم. برای این کار به صورت زیر از تابع eig (از توابع متلب) استفاده می کنیم:
[Uc,Ac]=eig(Cc);
این تابع ماتریس (مربعی) کواریانس مرکب را گرفته و در خروجی، { U }_{ c } را به عنوان ماتریس بردار های ویژه و { A }_{ c } را به عنوان ماتریس مقادیر ویژه بر می گرداند.
حال تبدیل سفید کنندگی را به صورت زیر در متلب انجام می دهیم:
G=(Ac^-0.5)*Uc';
حال با استفاده از کد زیر،{ S }_{ a } و { S }_{ b } را در متلب تعریف کرده، ماتریس بردار های ویژه و ماتریس مقادیر ویژه ی ماتریس { S }_{ a } را به دست می آوریم. سپس این ماتریس مقادیر ویژه (و متناسب با آن ماتریس بردار های ویژه) را به صورت نزولی مرتب می کنیم:
Sa=G*Ca*G';
Sb=G*Cb*G';
% Sort descending
indx=0;
[B,L1]=eig(Sa);
[L1,indx]=sort(diag(L1),'descend');
B=B(:,indx);
L1=diag(L1);
L2=eye(size(Sa,1))-L1;
همان طور که در خلاصه ی الگوریتم گفته شد، می توان ثابت کرد که { S }_{ a } و { S }_{ b } ماتریس بردار های ویژه ای را به اشتراک می گذارند! پس ماتریس B ، ماتریس بردار های ویژه ی هر دوی { S }_{ a } و { S }_{ a } است!!
در کُد بالا دیده می شود که برای مرتب کردن ماتریس مقادیر ویژه (و متناسب با آن ماتریس بردار های ویژه) ، از تابع sort (از توابع متلب) استفاده شده است. اگر ورودی تابع sort یک بردار باشد، اعضای بردار به صورت صعودی مرتب می شوند و اگر ورودی تابع یک ماتریس باشد، ستون های ماتریس به صورت صعودی مرتب می شوند. می دانیم که ماتریس مقادیر ویژه { L }_{ 1 }، یک ماتریس قطری است. پس به وسیله ی دستور diag\left( { L }_{ 1 } \right) ، عناصر قطر اصلی { L }_{ 1 } را تبدیل به یک بردار کرده و به تابع sort می دهیم. همچنین از رشته ی 'descend' در ورودی تابع استفاده کرده تا مرتب کردن به صورت نزولی باشد. خروجی این تابع { L }_{ 1 } مرتب شده (به صورت نزولی) و بردار index است.
این بردار در واقع ترتیب المان های { L }_{ 1 } در چیده شدن به صورت نزولی را در بر دارد. همچنین می دانیم که هر یک از ستون های ماتریس بردار های ویژه ی B مربوط به یکی از مقادیر ویژه هستند. پس برای مرتب کردن B متناسب با { L }_{ 1 }، از دستور B=B\left( :,index \right) استفاده می کنیم. سپس { L }_{ 1 } را به فرم ماتریسی خود بر می گردانیم. لازم به ذکر است که تابع diag ، بردار ورودی خود را به ماتریسی قطری تبدیل کرده و عناصر قطر اصلی ماتریس ورودی خود را به صورت بردار باز می گرداند.
در انتها { L }_{ 2 } را از کد زیر به دست می آوریم (Iماتریس همانی است) :
L2=eye(size(Sa,1))-L1;
فیلتر های فضایی را به صورت زیر تعریف می کنیم:
WT=B'*G;
در این کد، WT ماتریس فیلتر های فضایی مورد نظر می باشد.
همان طور که پیش تر در خلاصه ی الگوریتم CSP گفته شد، پس از تعیین فیلتر های فضایی مطلوب از روی داده های آموزشی دو کلاس، باید m سطر اول و m سطر آخر این فیلتر ها انتخاب شوند و سپس با استفاده از این سطر های انتخاب شده، داده ها فیلتر شوند. برای این انتخاب از کد زیر استفاده می کنیم:
WT_save=WT;
WT(m+1:N-m,:)=[];
برای احتیاط ماتریس فیلتر ها را در متغیری دیگر ریخته و سپس، به ازای همه ی ستون های WT ، سطر m+1 تا سطر N-m+1 آن را حذف می کنیم. می دانیم که به هر یک از سطر های WT ، یک فیلتر فضایی می گویند.
اکنون نوبت به فیلتر کردن داده ها و ایجاد بردار های ویژگی است. می دانیم که برای استفاده از کلاسیفایر ها، به دو بردار ویژگی نیاز داریم. یک بردار ویژگی به عنوان آموزش کلاسیفایر و یک بردار ویژگی برای آزمایش کلاسیفایر. پس در این جا می بایست یک بار داده های آموزشی را به وسیله ی WT فیلتر کنیم و از نتیجه ی به دست آمده بردار ویژگی بسازیم. سپس داده های آزمایشی را فیلتر کرده و بردار ویژگی آزمایشی را ایجاد کنیم. پس ابتدا برای ساختن بردار ویژگی برای داده های آموزشی از کد زیر استفاده می کنیم:
%training features
%class A
fa_training=[];
for i=1:n_train_a
summ=0;
ya=[];
Z=0;
Z=class_a_train{i};
ya=WT*Z;
%filtering
for j=1:size(ya,1)
fa_training(j,i)=var(ya(j,:));
summ=summ+var(ya(j,:));
end
fa_training(:,i)=log10((fa_training(:,i)/summ));
end
%Class B
fb_training=[];
for i=1:n_train_b
summ=0;
yb=[];
Z=0;
Z=class_b_train{i};
yb=WT*Z;
%filtering
for j=1:size(yb,1)
fb_training(j,i)=var(yb(j,:));
summ=summ+var(yb(j,:));
end
fb_training(:,i)=log10((fb_training(:,i)/summ));
end
در این کُد دو ماتریس fa_training و fb_training به عنوان ماتریس ویژگی آموزشی کلاسa و بردار ویژگی آموزشی کلاس a تعریف می شود که هر ستون آن ها، بردار ویژگی مربوط به یکی از نمونه های آموزشی کلاس مربوطه است. به طور نمونه برای کلاس a ، هر داده از داده های آموزشی در متغیر Z ریخته شده، به وسیله ی WT فیلتر می شود و ماتریس جدید فیلتر شده در ya قرار می گیرد. چون ماتریس WT پس از انتخاب انجام شده در مرحله ی قبل، یک ماتریس 2m\times N شده است و ماتریس Z ، ابعاد N\times T دارد، پس ya ماتریسی2m\times T می شود. حال طبق الگوریتم CSP ، باید به ازای هر سطر ماتریس فیلتر شده یک ویژگی بسازیم. ویژگی هایی که هر کدام واریانس یک سطر از ماتریس فیلتر شده ی ya هستند که نسبت به مجموع واریانس تمامی سطر ها نرمالیزه شده اند. متغیر summ مجموع واریانس ها را در خود نگه می دارد. به همین دلیل مقدار اولیه ی آن 0 بوده و زمانی که واریانس هر سطر به دست می آید، با مقدار قبلی summ جمع می شود. پس به ازای سطر jاُم از ماتریس فیلتر شده ی ya برای داده ی iاُم از سلول آموزشی (کلاس a) ، از کُد زیر استفاده می کنیم:
fb_training(j,i)=var(ya(j,:));
و پس از آن که واریانس تمامی سطر ها را به دست آوریم (مجموع آن ها را در summ است) ، واریانس های مربوط به آن داده را بر summ تقسیم می کنیم و از آنها لگاریتم در پایه ی نپر می گیریم. این کار به ازای تمامی اعضای سلول class_train انجام می شود.
پس از به دست آوردن بردار های ویژگی برای نمونه های آموزشی، می بایست بردار های ویژگی نمونه های آزمایشی را ایجاد کرد. روند این مرحله، کاملاً شبیه به دست آوردن بردار های ویژگی برای نمونه های آموزشی است. پس بدون توضیح اضافی از کُد زیر برای ایجاد fa_test و fb_test استفاده می کنیم.
%Test features
%class A
fa_test=[];
for i=1:n_test_a
summ=0;
Z=0;
ya=[];
Z=class_a_test{i};
ya=WT*Z;
for j=1:size(ya,1)
fa_test(j,i)=var(ya(j,:));
summ=summ+var(ya(j,:));
end
fa_test(:,i)=log10((fa_test(:,i)/summ));
end
%class B
fb_test=[];
for i=1:n_test_b
summ=0;
Z=0;
yb=[];
Z=class_b_test{i};
yb=WT*Z;
for j=1:size(yb,1)
fb_test(j,i)=var(yb(j,:));
summ=summ+var(yb(j,:));
end
fb_test(:,i)=log10((fb_test(:,i)/summ));
end
شبیه سازی روش CSP در این جا به پایان رسیده است. اما برای راحتی در آزمایشات، می خواهیم تابعی برای روش CSP تعریف کنیم که ضمن تمامی مراحل بالا، کلاس بندی داده ها به روش LDA و SVM را نیز انجام دهد.
%training features
%class A
fa_training=[];
for i=1:n_train_a summ=0;
ya=[];
Z=0;
Z=class_a_train{i};
ya=WT*Z;
%filtering
for j=1:size(ya,1)
fa_training(j,i)=var(ya(j,:));
summ=summ+var(ya(j,:));
end
fa_training(:,i)=log10((fa_training(:,i)/summ));
end
%Class B
fb_training=[];
for i=1:n_train_b
summ=0;
yb=[];
Z=0;
Z=class_b_train{i};
yb=WT*Z;
%filtering
for j=1:size(yb,1)
fb_training(j,i)=var(yb(j,:));
summ=summ+var(yb(j,:));
end
fb_training(:,i)=log10((fb_training(:,i)/summ));
end
function [ performance_LDA, performance_SVM,class_a_train] = CSP(varargin)
if nargin <2
error('Must have 2 classes for CSP!!')
end
if nargin>=4
perc=varagin{4};
m=varargin{3}
elseif nargin>=3
m=varargin{3};
perc=0.7;
else m=1;
perc=0.7;
end
%determination of class A and class B
class_a=varargin{1};
class_b=varargin{2};
%number of trials for each class
n_a=size(class_a,2);
n_b=size(class_b,2);
%mixing order of classes data and choosing training and test datasets
n_train_a=fix(n_a*perc);
n_test_a=n_a-n_train_a;
n_train_b=fix(n_b*perc);
n_test_b=n_b-n_train_b;
[train_a,test_a]=train_test(1:n_a,n_train_a,n_test_a);
[train_b,test_b]=train_test(1:n_b,n_train_b,n_test_b);
class_a_train={};
class_a_test={};
class_b_train={};
class_b_test={};
N=size(class_a{train_a(1)},1);
%class A
for i=1:n_train_a
class_a_train{i}=class_a{train_a(i)};
if size (class_a_train{i},1)~=N
error('All the trials must have a equal number of channels' )
end
end
for i=1:n_test_a
class_a_test{i}=class_a{test_a(i)};
if size (class_a_test{i},1)~=N
error('All the trials must have a equal number of channels' )
end
end
%class B
for i=1:n_train_b
class_b_train{i}=class_b{train_b(i)};
if size (class_b_train{i},1)~=N
error('All the trials must have a equal number of channels' )
end
end
for i=1:n_test_b
class_b_test{i}=class_b{test_b(i)};
if size (class_b_test{i},1)~=N
error('All the trials must have a equal number of channels' )
end
end
%covariance matrixs
Ca=0;
Z=0;
for i=1:n_train_a Z=class_a_train{i};
mn=mean(Z')';
for j=1:size(Z,2)
Z(:,j)=Z(:,j)-mn;
end
Ca=Ca+(Z*Z')/trace(Z*Z');
end
Ca=Ca/n_train_a;
Cb=0;
Z=0;
for i=1:n_train_b
Z=class_b_train{i};
mn=mean(Z')';
for j=1:size(Z,2)
Z(:,j)=Z(:,j)-mn;
end
Cb=Cb+(Z*Z')/trace(Z*Z');
end
Cb=Cb/n_train_b;
%Composite covariance matrix
Cc=Ca+Cb;
%Unitary similarity Transform ==> Cc=Uc*Ac*Uc'
[Uc,Ac]=eig(Cc);
%Whitening Transform
G=(Ac^-0.5)*Uc';
Sa=G*Ca*G';
Sb=G*Cb*G';
%So Sa and Sb share Common eigenvector matrix
% Sort descending
indx=0;
[B,L1]=eig(Sa);
[L1,indx]=sort(diag(L1),'descend');
B=B(:,indx);
L1=diag(L1);
L2=eye(size(Sa,1))-L1;
%Spatial Filter
WT=B'*G;
%Choosing 2m spatial filter
WT_save=WT;
WT(m+1:N-m,:)=[];
%Making Features
%training features
%class A
fa_training=[];
for i=1:n_train_a
summ=0;
Z=0;
ya=[];
Z=class_a_train{i};
ya=WT*Z; %filtering
for j=1:size(ya,1)
fa_training(j,i)=var(ya(j,:));
summ=summ+var(ya(j,:));
end
fa_training(:,i)=log10((fa_training(:,i)/summ))
end
%Class B
fb_training=[];
for i=1:n_train_b
summ=0;
Z=0;
yb=[];
Z=class_b_train{i};
yb=WT*Z; %filtering
for j=1:size(yb,1)
fb_training(j,i)=var(yb(j,:));
summ=summ+var(yb(j,:));
end
fb_training(:,i)=log10((fb_training(:,i)/summ));
end
%Test features
%class A
fa_test=[];
for i=1:n_test_a
summ=0;
Z=0;
ya=[];
Z=class_a_test{i};
ya=WT*Z; %filtering
for j=1:size(ya,1)
fa_test(j,i)=var(ya(j,:));
summ=summ+var(ya(j,:));
end
fa_test(:,i)=log10((fa_test(:,i)/summ));
end
%class B
fb_test=[];
for i=1:n_test_b
summ=0;
Z=0;
yb=[];
Z=class_b_test{i};
yb=WT*Z; %filtering
for j=1:size(yb,1) %tap haye har nemooneye test
fb_test(j,i)=var(yb(j,:));
summ=summ+var(yb(j,:));
end
fb_test(:,i)=log10((fb_test(:,i)/summ));
end
%classification
training_set=(horzcat(fa_training,fb_training))';
test_set=(horzcat(fa_test,fb_test))';
cls_train='';
cls_train(1:size(fa_training,2))='A';
cls_train((size(fa_training,2)+1):(size(fb_training,2)+size(fa_training,2)))='B';
cls_train=cls_train';
%LDA
cls_test_LDA=classify(test_set,training_set,cls_train);
%SVM
svm_struct=svmtrain(training_set,cls_train);
cls_test_SVM = svmclassify(svm_struct,test_set);
%Performances
n_tot=(size(fa_test,2)+size(fb_test,2));
performance_LDA=100*(length(find(cls_test_LDA(1:n_test_a)=='A'))+length(find(cls_test_LDA((n_test_a+1):n_tot)=='B')))/n_tot;
performance_SVM=100*(length(find(cls_test_SVM(1:n_test_a)=='A'))+length(find(cls_test_SVM((n_test_a+1):n_tot)=='B')))/n_tot;
هر سه مشتق فیلترهای فضایی نیز به همین منوال می باشد. لذا با تغییرات جزیی در همین کد بالا قابل اجراست به همین منظور برای جلوگیری از دوباره کاری فقط نتایج آنها اشاره شده است.
۱۰. نتایج
الف) داده های ساختگی
ابتدا نتایج کلاس بندی داده های ساختگی که در قسمت مربوطه مشخص شدند را با هم بررسی می کنیم. برای تایید نتایج، آزمایشات را 200 بار تکرار می کنیم:
ب) داده های واقعی
در این قسمت، داده های واقعی معرفی شده در قسمت مربوطه را کلاس بندی می کنیم. نتایج این آزمایشات را 100 بار تکرار کرده و میانگین و ماکزیمم هر کدام را نمایش می دهیم.
۱۱. کارهای آینده
برای بهبود این روش به نظرم قبل از استفاده از کلاسیفایرها شاید یک پری پراسسینگ با تبدیل ویولت روی داده ها می تواند نتایج را بهتر کند.
۱۲. مراجع
[1] Müller-Gerking, Johannes, Gert Pfurtscheller, and Henrik Flyvbjerg. "Designing optimal spatial filters for single-trial EEG classification in a movement task." Clinical neurophysiology 110.5 (1999): 787-798.
[2] Ramoser, Herbert, Johannes Muller-Gerking, and Gert Pfurtscheller. "Optimal spatial filtering of single trial EEG during imagined hand movement." Rehabilitation Engineering, IEEE Transactions on 8.4 (2000): 441-446.
[3] Koles, Zoltan Joseph. "The quantitative extraction and topographic mapping of the abnormal components in the clinical EEG." Electroencephalography and clinical Neurophysiology 79.6 (1991): 440-447.
[4] Ramoser, Herbert, Johannes Muller-Gerking, and Gert Pfurtscheller. "Optimal spatial filtering of single trial EEG during imagined hand movement." Rehabilitation Engineering, IEEE Transactions on 8.4 (2000): 441-446.
[5] Park, Cheolsoo, C. CHEONG TOOK, and D. Mandic. "Augmented Complex Common Spatial Patterns for Classification of Noncircular EEG from Motor Imagery Tasks." (2014): 1-1.
[6] Dornhege, Guido, et al. "Boosting bit rates in noninvasive EEG single-trial classifications by feature combination and multiclass paradigms." Biomedical Engineering, IEEE Transactions on 51.6 (2004): 993-1002.
[7] Lemm, Steven, et al. "Spatio-spectral filters for improving the classification of single trial EEG." Biomedical Engineering, IEEE Transactions on 52.9 (2005): 1541-1548.
[8] Haykin, Simon. "Adaptive filter theory." 2nd. ed., Prentice-Hall, Englewood Cliffs, NJ (1991).
[9] Huang, Norden E., et al. "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis." Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454.1971 (1998): 903-995.
[10] Yuan, Han, et al. "Cortical imaging of event-related (de) synchronization during online control of brain-computer interface using minimum-norm estimates in frequency domain." Neural Systems and Rehabilitation Engineering, IEEE Transactions on 16.5 (2008): 425-431.
[11] Pfurtscheller, Gert, and Fernando H. Lopes da Silva. "Event-related EEG/MEG synchronization and desynchronization: basic principles." Clinical neurophysiology 110.11 (1999): 1842-1857.
[12] Falzon, O., K. P. Camilleri, and J. Muscat. "Complex-valued spatial filters for task discrimination." Engineering in Medicine and Biology Society (EMBC), 2010 Annual International Conference of the IEEE. IEEE, 2010.
[13] Chebotarev, Alexander M., and Alexander E. Teretenkov. "Singular value decomposition for the Takagi factorization of symmetric matrices." Applied Mathematics and Computation 234 (2014): 380-384.
Motor Imagery
Common Spatial Patterns (CSP)
Brain-Computer Interfaces (BCI)
Electroencephalogram
Abnormality
Identity Matrix
Whitening Transform
Automatic spectral filter selection
Analytic Signal-based CSP
Augmented Complex Common Spatial Patterns
Strong Uncorrelating Transform Augmented Complex Common Spatial Patterns
Non-Circular
Correlation
Complex Statistics
Pseudocovariance
Diagonalization
‘Analytic Signal’-based CSP
Analytic
این داده ها در آدرس http://www.bbci.de/competition/iv موجود می باشند.