قرار است تصور حرکتی¹ را از سیگنال EEG تفکیک کنیم. با این روش می توانیم از سیگنال حاصله جهت استفاده در به حرکت درآوردن بازوهای مکانیکی استفاده کنیم. برای مثال می‌توان با ایجاد تصور حرکت دست‌ها و تفکیک و انتقال این سیگنال به فردی که از ناحیه دست قطع عضو گردیده یا نرون‌های عصبی آن از کار افتاده، دست مصنوعی و یا دست خود فرد را به حرکت درآوریم.

۱. مقدمه

در این گزارش، روش CSP به همراه سه تعمیم آن یعنی ACSP، ACCSP و SUTCCSP به طور کامل شرح داده می‌شود. الگو های فضایی مشترک روشی محبوب در جداسازی بین کلاس‌ها است. اعمال این روش به داده‌ها پیش از کلاس‌بندی، عملکرد را بالا تر می‌برد. این روش ماهیتی چند کاناله دارد و اغلب برای چند کانال سیگنال به کار می‌رود. یکی از مشکلات این روش عدم دخالت اطلاعات مربوط به فاز سیگنال‌ها در آن است. سیگنال‌های جهان ما اکثراً غیر‌دایره ای هستند. برای مثال در سیگنال EEG ، به دلیل اختلاف توان بین کانال‌ها و همبستگی بین آن‌ها این موضوع رویت می‌شود. سه تعمیمی که مورد بحث قرار می‌گیرند، سعی در اضافه کردن اطلاعات فاز و جدا سازی بیشتر و بهتر بین کلاس‌ها را دارند.

۲. کارهای مرتبط

روش الگو های فضایی مشترک ² روشی محبوب و پر کاربرد در جدا سازی بین داده‌های دو کلاس، مخصوصاً در تفکیک فعالیت‌های ذهنی در واسط‌های مغز کامپیوتر ³ مبتنی بر داده‌های EEG ⁴ است. عملکرد بالا در کلاس‌بندی و فیلترهای فضایی مفید به دست آمده از این روش، این الگوریتم را مشهور ساخته است [1] و [2].
این الگوریتم برای اولین بار برای تشخیص غیر طبیعی بودن ⁵ سیگنال EEG به کار گرفته شد [3] و سپس در جدا سازی بین کلاس‌های مربوط به واسط‌های مغز کامپیوتر و الگو های حرکتی ⁶ نیز مورد استفاده قرار گرفت [1] و [2]. به دلیل ماهیت چند کاناله بودن، این روش کاربرد فراوانی برای سیگنال های EEG دارد.
روش CSP، با اعمال فیلترهای فضایی به ورودی‌ها، واریانس سیگنال‌ها را در کلاس اول ماکزیمم و به طور همزمان در کلاس دیگر مینیمم می کند و سپس از سیگنال‌های فیلتر شده ویژگی‌های کلاس اول را می‌سازد. این کار به طور عکس نیز اتفاق می افتد. یک سری از فیلترهای فضایی، واریانس سیگنال را در کلاس دوم ماکزیمم و به طور همزمان در کلاس اول مینیمم می‌کنند. سپس با استفاده از این سیگنال‌های فیلتر شده، ویژگی‌های کلاس دوم ایجاد می‌شود. برای به دست آوردن این فیلترهای فضایی، CSP مجموع ماتریس‌های کواریانس دو کلاس را تبدیل به ماتریس همانی ⁷ می‌کند. این کار با اعمال تبدیلی به نام سفید کنندگی⁸ انجام می‌شود. با این کار، افزایش واریانس یک کلاس، واریانس کلاس دیگر را کاهش می‌دهد.
تعمیم‌های زیادی از این روش ارائه شده است که هر کدام سعی در بهبود کارایی و یا رفع محدودیتی در الگوریتم اصلی دارند. برای مثال در [5] حالت چند کلاسه‌ی CSP و در [6] انتخاب فیلتر طیفی به طور خود کار⁹ مطرح شده است. در این گزارش، به غیر از روش CSP ، سه تعمیم آن، یعنی ACSP ¹⁰ ، ACCSP ¹¹ و SUTCCSP¹² نیز شبیه سازی می شود.
این سه تعمیم، به غیر از اطلاعات دامنه‌ی سیگنال‌ها، اطلاعات فاز سیگنال‌ها را نیز در تفکیک کلاس‌ها دخیل می‌کنند. چرا که بسیاری از سیگنال‌های دنیای واقعی، سیگنال‌های غیر دایره ای ¹³ هستند. به عبارت دیگر، بین سیگنال‌های دو کانال، همبستگی ¹⁴ وجود دارد. روش ACSP فقط اطلاعات فاز سیگنال ها را در الگوریتم CSP دخیل می‌کند. ولی روش‌های ACCSP و SUTCCSP به غیر از این کار، از آمار درجه دوم مربوط به مجموعه اعداد مختلط¹⁵ استفاده می کنند. مباحث آماری که به ما می گوید، برای اعداد متغیر های تصادفی مختلف به غیر از کواریانس، به شبه کواریانس ¹⁶ هم نیاز داریم. شبه کواریانس داده های غیر دایره ای صفر نیستند! روش ACCSP ماتریسی حاصل از ترکیب کواریانس و شبه کواریانس را برای قطری‌سازی ¹⁷ استفاده می‌کند. در حالی که روش SUTCCSP ، به طور همزمان ماتریس های کواریانس و شبه کواریانس را قطری می سازد.

۳. الگو های فضایی مشترک برای سیگنال های حقیقی و تحلیلی

الگوریتم CSP با روش CSP مختلط یا همان ACSP¹⁸ تفاوتی ندارد. ولی در روش ACSP سیگنال های ورودی پیش از اعمال به الگوریتم، به صورت تحلیلی¹⁹ در می آیند. پس ابتدا فرض کنید که سیگنال ورودی، سیگنالی حقیقی است.
ماتریس داده های { 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 }\quad =\quad I وجود دارد.
به تبدیل G\quad =\quad \sqrt { { \Lambda }_{ c }^{ -1 } } { U }_{ c }^{ H } ، تبدیل سفید کننده می گویند. چرا که:
...
یعنی این تبدیل، کواریانس مرکب را قطری ( ماتریس همانی ) می کند. پس می توان نوشت:
....
که در آن I نشان دهنده ی ماتریس همانی است.
حال اگر ... را تعریف کنیم و B یک ماتریس متعامد نرمال باشد، آنگاه داریم:
...
چرا که می دانیم...
یعنی می توان گفت که ماتریس های ... ماتریس بردار های ویژه ی مشترکی دارند. به طوری که:
...
حال با توجه به رابطه ی (10) ، اگر ماتریس مقادیر ویژه ی ... را به صورت نزولی مرتب کنیم، به طور هم زمان ... به صورت صعودی مرتب می شود. با این کار، بردار ویژه ای از ماتریس B که مربوط به بزرگ ترین مقدار ویژه ی ... می شود،
مربوط به کوچک ترین مقدار ویژه ی ... است. یعنی یک بردار ویژه از B (مربوط به بزرگترین مقدار ویژه ی ... )، به طور همزمان واریانس کلاس a را ماکزیمم و واریانس کلاس b را مینیمم می کند. پس این فیلتر فضایی را به صورت زیر تعریف می کنیم:
...
و با توجه به این تعریف، رابطه ی (9) و (11) داریم:
...
حال ماتریس داده ها را روی این فیلتر فضایی تصویر می کنیم تا داده های جدیدی به دست آوریم. در حقیقت داده های ورودی را توسط فیلتر های فضایی فیلتر می کنیم. یعنی:
...
هر کدام از ستون های ماتریس ...را یک فیلتر فضایی می نامند.
برای جدا سازی بین دو کلاس، از واریانس سیگنال های فیلتر شده در رابطه ی (12) استفاده می کنیم. می دانیم که هر سطر از ماتریس V ، یک سیگنال جدید است. به منظور کاهش بعد فضای ویژگی، m سطر اول V ( که وابسته به مقدار ویژه ی بزرگ ... است، یعنی بیشترین واریانس را در کلاس a ایجاد می کنند ) و m سطر آخر V ( که وابسته به m مقدار ویژه ی بزرگ ... است یعنی بیشترین واریانس را در کلاس b ایجاد می کنند ) را انتخاب می کنیم. به عبارت دیگر، سطر های
ماتریس ... ، اختلاف واریانس بین دو کلاس را ماکزیمم می کنند. چرا که m سطر اول بیشترین واریانس در کلاس a و کمترین واریانس در کلاس b را دارند و m سطر آخر، دارای بیشترین واریانس در کلاس b و کمترین واریانس در کلاس a هستند. ویژگی های مورد نظر ما از رابطه ی زیر به دست می آید:
...
که در این رابطه p ، ویژگی مربوط به سطر pاُم ماتریس و منظور از ... ، واریانس سطر 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)، به صورت زیر تعریف می شود:
...
در رابطه ی بالا، ... و ... به ترتیب، دامنه و فاز سیگنال مختلط z(t) هستند که به صورت زیر تعریف می شوند:
...
و ... ، تبدیل هیلبرت سیگنال ورودی است:
...
که منظور از P.V. ، مقدار اصلی کوشی ( روشی برای مقدار دادن به انتگرال های ناسره ) است.
در خصوص استفاده از تبدیل هیلبرت باید به این موضوع توجه داشت که این تبدیل فقط برای سیگنال هایی با باند فرکانسی باریک استفاده می شود و نمی تواند برای سایر سیگنال ها، اطلاعات کاملی راجع به محتوای فرکانسی آن ها ارائه کند. همچنین زمانی که تغییرات دامنه سریع تر از تغییرات فاز است (مثلاً سیگنال هایی که میانگین صفر ندارند)، تبدیل هیلبرت، فرکانس لحظه ای منفی ایجاد می کند [9] . به همین جهت، برای رفع محدودیت پهنای باند تبدیل هیلبرت، باید از تبدیل فوریه و یا EMD در کنار ACSP استفاده کرد [5] .
طبق [5] ، المان های بردار ویژگی به دست آمده از ACSP برای هر کلاس را می توان از روابط زیر به دست آورد:
...

۵. الگو های فضایی مشترک مختلط افزوده (ACCSP)

برای یک بردار میانگین صفر مختلط مانند ، ماتریس کواریانس به صورت ... محاسبه می شود. مباحث آماری مربوط به اعداد مختلط، به طور کامل تعمیمی از مباحث آماری اعداد حقیقی نیستند! پس ماتریس کواریانس C به طور کامل آمار درجه دوم بردار z را توصیف نمی کند!! پس برای برداری مانند z که مقادیر مختلط دارد، به ماتریس شبه کواریانس یعنی ... نیز احتیاج داریم.
متغیر تصادفی ... را در نظر بگیرید. کواریانس این ماتریس برابر ... می شود ( که همیشه مثبت است مگر آن که ... باشد ). منظور از ... ، عملگر مزدوج مختلط است. این در حالیست که شبه کواریانس که به صورت ... تعریف می شود، هنگامی بی اثر است که به طور همزمان، ... با یک دیگر ناهمبسته بوده و واریانس آن ها برابر باشد( [ ] [ ]) . به سیگنال هایی که شبه کواریانس صفر دارند، دایره ا ی مرتبه دوم یا مناسب می گویند. با این حال، حتی اگر داده ای در دنیای واقعی، دایره ا ی باشد، به دلیل مشاهده ی سیگنال در پنجره های کوچک، نویز های نا همسانگر و توان متفاوت کانال های مختلف، شبه کواریانس برابر صفر نیست [5] . شکل 1، نمایشی هندسی از توضیح داده های دایره ا ی و غیر دایره ا ی را نشان می دهد. میزان دایره ا ی بودن بستگی به میزان همبستگی بین قسمت حقیقی و موهومی داده ها دارد.همبستگی کمتر، دایره ا یبودن بیشتر و همبستگی بیشتر، غیر دایره ا یبودن بیشتری را ایجاد می کند [5] .

از این رو، برای داده های مختلط، باید خواص آماری مشترک ... مورد بررسی قرار گیرد. پس فُرم افزوده ای از متغیر مختلط z به صورت ̂ [ ] ،تعریف می شود. ماتریس کواریانس این متغیر افزوده ( ماتریس افزوده ) به صورت زیر است:
...
این ماتریس کواریانس افزوده، خصوصیات آماری مرتبه ی دوم کامل کواریانس و شبه کواریانس را در بر دارد.
روش ACSP ( و به طور قطع روش CSP )، غیر دایره ای بودن داده ها را در نظر نمی گیرد. به همین دلیل روش الگوهای فضایی مشترک مختلط افزوده (ACCSP) مطرح می شود [5] . در این روش به جای ماتریس کواریانس، فیلتر های فضایی از طریق ماتریس کواریانس افزوده به دست می آید. ماتریس داده های تبدیل شده به فرم تحلیلی ... را داریم ( منظور از C ، مجموعه اعداد مختلط است ). این ماتریس ها را به فرم افزوده ی خود ...تبدیل کرده و کواریانس آنها را
به دست می آوریم:
...
که ̂ [ ] و ̂ [ ] هستند. ادامه ی روند، مانند الگوریتم اولیه ی CSP است. ماتریس کواریانس مرکب افزوده تعریف می شود. تبدیل سفید کننده روی آن اعمال می شود و ماتریس بردار های ویژه ی مشترک بین کلاس a و b به دست می آید. ماتریس فیلتر فضایی افزوده،...، از رابطه ی (11) به دست می آید. حال برای هر داده ی ورودی Z ، فرم افزوده ی ورودی،...، را به دست می آوریم و سیگنال ها را با استفاده از رابطه ی زیر، فیلتر می کنیم:
...
ماتریس...مقادیر مختلط دارد و مانند قبل، 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 }} را به صورت ... تجزیه کرد. قابل ذکر است که ماتریس ... می تواند مربعی نباشد، ... ماتریس های واحد و ... هم ماتریسی با ابعاد ... است که عناصر آن جذر مقادیر تکین ...هستند. اگر ... ماتریسی مربعی باشد ( که هست! )، ... ماتریسی قطری می شود.
همان طور که بیان شد، داریم:
̅ ̅ ( )
ابتدای دو طرف رابطه را ... در ضرب می کنیم :
...
در رابطه ی بالا از خاصیت واحد بودن ماتریس ... استفاده کردیم. حال انتهای دو طرف رابطه را در ... ضرب می کنیم :
( )
حال با تعریف ... ، به رابطه ی زیر می رسیم:
( )
... یک ماتریس واحد است. چرا که از ضرب دو ماتریس واحد حاصل شده است. پس رابطه ی ... بر قرار است.
همچنین می توان ثابت کرد که اگر رابطه ی ... وجود داشته باشد، رابطه ی ( ) نیز موجود است!
چرا که می توان نوشت:
( ) ( ) ( ) ( ) ( ) ( )
به همین طریق می توان به ( ) نیز رسید!
حال به رابطه ی (28) باز می گردیم:
̅ ( ) ( ) ( )
اگر ... را تعریف کنیم، به عبارت زیر می رسیم:
̅ ( )
به رابطه ی (30) تبدیل تاکاگی یا SSVD می گویند. باز هم مشخص است که Y یک ماتریس واحد است. چرا که از ضرب دو ماتریس واحد به دست می آید.
همان طور که پیش تر گفته شد، به دنبال ماتریس Q می گردیم که در رابطه ی (23) صدق کند. اگر ... را تعریف کنیم، آنگاه:
( ) ̅ ( )
پس ملاحظه می شود که با تعریف ... ، ماتریس مورد نظر در رابطه ی (23) را یافتیم. برای کواریانس می توان نوشت:
( )
اگر ... را تعریف کنیم، مانند الگوریتم اصلی CSP ، می توان گفت که ماتریس های ... ماتریس بردار های ویژه ی مشترکی دارند. به طوری که:
, ( ) ( )
به همین ترتیب، می خواهیم برای شبه کواریانس دو کلاس هم، ماتریس بردار های ویژه ی مشترکی تخمین بزنیم. پس:
...
ابتدا و انتهای عبارت بالا را در ... ضرب می کنیم:
( )
در رابطه ی بالا، ... یک ماتریس قطری با عناصر حقیقی است. پس رابطه ی ... برقرار است. برای یافتن ماتریس بردار های ویژه ی مشترک شبه کواریانس ها، به صورت زیر عمل می کنیم:
̂ ( )
اگر ̂ ̂ ̂ و ̂ ̂ ̂ را تعریف کنیم، رابطه ی (33) به صورت ̂ ̂ در می آید. پس مانند مباحث قبل، می توان برای ماتریس های ̂ و ̂ ، ماتریس بردار های ویژه ی مشترکی تخمین زد. یعنی:
̂ ̂ ̂ ̂ , ̂ ̂ ̂ ̂ ( ̂ ̂ ) ( )
در حقیقت با اعمال تبدیل SUT ، به طور هم زمان ماتریس های کواریانس و شبه کواریانس را به ماتریس همانی تبدیل کردیم ( رابطه ی (32 )و (34)). با این کار، اختلاف واریانس بین دو کلاس ماکزیمم خواهد شد. این بار به جای یک دسته فیلتر فضایی، برای هر کدام از کواریانس و شبه کواریانس، یک دسته فیلتر فضایی تعریف می کنیم:
̂ ̂ ̂ ( )
و به همین ترتیب، سیگنال ها به صورت زیر فیلتر می شوند:
̂ ̂ ( )
تا بتوانیم المان های بردار های ویژگی را طبق (17) و (18) یا (13) به دست آوریم.
زمانی که ماتریس کواریانس مرکب را با اعمال تبدیل سفید کنندگی، به ماتریس همانی تبدیل می کنیم، مجموع ماتریس بردار های ویژه کلاس و همانی می شود. این کار باعث می شود که سطر اول ماتریس فیلتر شده ی ، بیشترین واریانس را در کلاس و کمترین واریانس را در کلاس داشته باشد و همین طور سطر آخر این ماتریس فیلتر شده بیشترین واریانس را در کلاس و کمترین واریانس را در کلاس دارد. حال در این جا با اعمال تبدیل SUT ، به غیر از ماتریس کواریانس مرکب، شبه کواریانس مرکب نیز به ماتریسی همانی تبدیل می شود. پس طبیعی است که سیگنال های فیلتر شده به وسیله ی فیلتر های فضایی شبه کواریانس ) ̂ ( نیز، خاصیت های را در مورد اختلاف واریانس ها در کلاس های مختلف داشته باشند.
می دانیم که تبدیل تاکاگی منجر به رابطه ی (23) شد. نگاهی به این رابطه می اندازیم:
[ ] [( )( ) ]
می دانیم که مقادیری مختلط دارد. پس مقادیر ماتریس به صورت زیر بیان می شوند:
λ [ ] [ ] [ ] [ ]
این در حالیست که ماتریس ای که از تبدیل تاکاگی به دست می آید، مقادیری حقیقی دارد. این موضوع نشان می دهد که مقادیر ماتریس λ که به صورت [ ] [ ] هستند، اطلاعات اضافی راجع به اختلاف توان بین قسمت های حقیقی و موهومی منبع ورودی به ما می دهند. پس روش 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] .²⁰
به دلیل مفاهیم نئوروفیزیولوژیکی ، به جای استفاده از 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 است. در حقیقت این تابع، با گرفتن کل سیگنال و اطلاعات آن به عنوان ورودی، قسمت هایی از سیگنال که تصور مربوط به کلاس اول و یا دوم انجام شده را در سلول مربوط به آن کلاس می ریزد. چرا که هدف ما در این پروژه پردازش و کلاس بندی داده ها به صورت آنلاین نیست!
هدف، بررسی عملکرد روش های به کار گرفته شده است.

۹. کارهای آینده

شبیه سازی فیلتر ها فضایی
در این قست نحوه ی شبیه سازی و فیلترینگ، انتخاب ویژگی و کلاس بندی به وسیله ی هر روش را شرح می دهیم.

۱۰. مراجع

[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.