אלגוריתם ה-FFT של קולי-טוקי

אלגוריתם קולי-טוקי, על שמם של ג'יימס קולי וג'ון טוקי, הוא אלגוריתם התמרת פורייה מהירה (FFT) הנפוץ ביותר. הוא מבטא מחדש, באופן רקורסיבי, את התמרת פורייה הבדידה (DFT) בגודל פריק שרירותי $N=N_{1}N_{2}$ במונחים של N₁ התמרות DFT קטנות יותר, כל אחת בגודל N₂, וזאת על מנת להפחית את זמן החישוב ל-O(N log N) עבור N פריק מאוד.

מכיוון שאלגוריתם קולי-טוקי מפרק את ה-DFT להתמרות DFT קטנות יותר, ניתן לשלבו עם כל אלגוריתם DFT אחר. לדוגמה, ניתן להשתמש באלגוריתמים של Rader או Bluestein לטיפול בגורמים ראשוניים גדולים שלא ניתן לפרק על ידי קולי-טוקי, או שאפשר להשתמש באלגוריתם FFT של גורם ראשוני ליעילות רבה יותר בהפרדת גורמי מספרים זרים .

האלגוריתם, כולל יישומו הרקורסיבי, הומצא בידי קרל פרידריך גאוס. 160 שנה מאוחר יותר, קולי וטוקי פיתחו אותו מחדש והפכו אותו לפופולרי .

היסטוריה

אלגוריתם זה ויישומו הרקורסיבי, הומצאו כאמור בסביבות 1805 בידי קרל פרידריך גאוס, שהשתמש בו לחישובי אינטרפולציה של מסלולי האסטרואידים פאלאס ויונו, אך עבודתו לא זכתה להכרה רחבה (היא התפרסמה רק לאחר מותו, ובנאו-לטינית ). גאוס, בכל אופן, לא ניתח את זמן החישוב האסימפטוטי. כמה צורות מוגבלות הומצאו מחדש מספר פעמים במהלך המאה ה-19 ותחילת המאה ה-20. FFT הפך לפופולרי לאחר שג'יימס קולי מ- IBM וג'ון טוקי מפרינסטון פרסמו ב-1965 מאמר הממציא מחדש את האלגוריתם ומתאר כיצד לבצע אותו בנוחות על מחשב.

על פי המדווח, טוקי הגה את הרעיון במהלך פגישה של הוועדה המייעצת למדע של הנשיא קנדי, שדנה בדרכים לאיתור ניסויים בנשק גרעיני בברית המועצות על ידי שימוש בסייסמומטרים הממוקמים מחוץ למדינה. חיישנים אלה ייצרו סדרות זמן סייסמולוגיות. הבעיה הייתה שניתוח נתונים אלה ידרוש אלגוריתמים מהירים לחישוב DFT בשל מספר החיישנים ואורך סדרות הזמן. משימה זו הייתה קריטית לצורך אשרור האיסור המוצע על ניסויים גרעיניים, כך שניתן יהיה לזהות כל הפרה ללא צורך בביקור במתקנים סובייטיים. משתתף אחר באותה פגישה היה ריצ'רד גארווין מ-IBM, שהכיר בפוטנציאל של השיטה ודאג ליצור קשר בין טוקי וקולי. גארווין וידא שקולי לא ידע את המטרה המקורית. במקום זאת, נאמר לו שזה נחוץ כדי לקבוע מחזוריות של כיווני הספין בגביש תלת-ממדי של הליום-3 . קולי וטוקי פרסמו לאחר מכן את המאמר המשותף שלהם, שאומץ במהירות ובאופן נרחב, וזאת בשל הפיתוח הבו-זמני של ממירי אנלוגי לדיגיטלי המסוגלים לדגום בקצבים של עד 300KHz.

העובדה שגאוס תיאר את אותו אלגוריתם (הגם שמבלי לנתח את עלותו האסימפטוטית) נודעה רק מספר שנים לאחר פרסום המאמר של קולי-טוקי בשנת 1965. המאמר שלהם ציין כהשראה רק את עבודתו של IJ Good שנקראת כיום אלגוריתם FFT של גורם ראשוני (PFA); למרות שהאלגוריתם של Good נחשב בתחילה כשקול לאלגוריתם קולי-טוקי, התברר במהרה ש-PFA הוא אלגוריתם שונה לחלוטין (עובד רק עבור גדלים שיש להם גורמים שהם מספרים זרים ומסתמך על משפט השאריות הסיני, בניגוד לתמיכה בכל גודל פריק ב-קולי-טוקי).

radix-2 FFT

FFT מסוג Radix-2 הוא הגרסה הפשוטה והנפוצה ביותר של אלגוריתם קולי-טוקי. יישומי קולי-טוקי בעלי אופטימיזציה גבוהה משתמשים לרוב בגרסאות אחרות של האלגוריתם, כפי שיתואר בהמשך. בשיטה זו גודל הקלט N הוא חזקה של שתיים. בכל שלב רקורסיבי, Radix-2 מחלק DFT בגודל N לשני DFT שזורים (interleaved) בגודל N/2.

התמרת פורייה הבדידה (DFT) מוגדרת בידי הנוסחה:

$X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk},$

כאשר $k$ הוא מספר שלם הנע בין 0 ל $N-1$ .

אנו מניחים כאן ש-N הוא חזקה של שתיים. Radix-2 מחשב תחילה את ה-DFT של ערכי הקלט בעלי האינדקס הזוגי $(x_{2m}=x_{0},x_{2},\ldots ,x_{N-2})$ ושל ערכי הקלט בעלי האינדקס האי-זוגי $(x_{2m+1}=x_{1},x_{3},\ldots ,x_{N-1})$, ולאחר מכן משלב את התוצאות הללו כדי לחשב את ה-DFT של הרצף כולו. את התהליך הזה ניתן לבצע באופן רקורסיבי ובכך להוריד את זמן הריצה הכולל ל-O(NlogN). היישום יכול לרוב, הגם שלא תמיד, לבחור את מספר נקודות הדגימה N באופן חופשי (למשל בעזרת שיטות כגון שינוי קצב הדגימה, גודל החלון, או zero-padding). לכן, לעיתים קרובות, הגודל של חזקה של שתיים אינו מגבלה חשובה.

אלגוריתם radix-2 מארגן מחדש את ה-DFT של הפונקציה $x_{n}$ בשני חלקים: סכום של האינדקסים הזוגיים $n={2m}$ וסכום של האי-זוגיים $n={2m+1}$ :

${\begin{matrix}X_{k}&=&\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N}}(2m)k}+\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N}}(2m+1)k}\end{matrix}}$

בסכום השני אפשר למצוא מכפיל משותף, $e^{-{\frac {2\pi i}{N}}k}$, כמוצג במשוואה למטה. שני הסכומים במשוואה זו הם ה-DFT של המרכיבים הזוגיים $x_{2m}$ וה-DFT של המרכיבים האי-זוגיים $x_{2m+1}$ של הפונקציה $x_{n}$ .

נסמן את ה-DFT של המרכיבים הזוגיים $x_{2m}$ כ-$E_{k}$.

נסמן את ה-DFT של המרכיבים האי-זוגיים $x_{2m+1}$ כ-$O_{k}$.

ונקבל:

${\begin{matrix}X_{k}=\underbrace {\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}} _{\mathrm {DFT\;of\;even-indexed\;part\;of\;} x_{n}}{}+e^{-{\frac {2\pi i}{N}}k}\underbrace {\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}} _{\mathrm {DFT\;of\;odd-indexed\;part\;of\;} x_{n}}=E_{k}+e^{-{\frac {2\pi i}{N}}k}O_{k}\qquad {\text{ for }}k=0,\dots ,{\frac {N}{2}}-1.\end{matrix}}$

השוויון אמנם מתקיים עבור $k=0,\dots ,N-1$, אך הנקודה המכרעת היא שאת הערכים $E_{k}$ ו $O_{k}$ יש לחשב רק עבור $k=0,\dots ,{\frac {N}{2}}-1$. הודות למחזוריות של המעריך המרוכב, $X_{k+{\frac {N}{2}}}$ מחושב גם הוא מ $E_{k}$ ו $O_{k}$, כלהלן :

${\begin{aligned}X_{k+{\frac {N}{2}}}&=\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}m(k+{\frac {N}{2}})}+e^{-{\frac {2\pi i}{N}}(k+{\frac {N}{2}})}\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}m(k+{\frac {N}{2}})}\\&=\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}e^{-2\pi mi}+e^{-{\frac {2\pi i}{N}}k}e^{-\pi i}\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}e^{-2\pi mi}\\&=\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}-e^{-{\frac {2\pi i}{N}}k}\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}\\&=E_{k}-e^{-{\frac {2\pi i}{N}}k}O_{k}\end{aligned}}$

אנחנו יכולים לכתוב מחדש את $X_{k}$ ו $X_{k+{\frac {N}{2}}}$ כלהלן:

${\begin{matrix}X_{k}&=&E_{k}+e^{-{\frac {2\pi i}{N}}{k}}O_{k}\\X_{k+{\frac {N}{2}}}&=&E_{k}-e^{-{\frac {2\pi i}{N}}{k}}O_{k}\end{matrix}}$

תוצאה זו, ביטויו הרקורסיבי של DFT באורך N במונחי שני DFT באורך N/2, היא הליבה של התמרת פורייה המהירה radix-2. מקור המהירות של האלגוריתם הוא השימוש החוזר בתוצאות חישובי ביניים כדי לחשב כמה תוצאות DFT. התוצאות הסופיות מתקבלות על ידי שילוב +/− של $E_{k}$ ו $O_{k}\exp(-2\pi ik/N)$, שהוא פשוט DFT בגודל 2 (הנקרא לעיתים פרפר בהקשר הזה); כאשר מכלילים לבסיסים גדולים יותר כמתואר להלן, ה-DFT בגודל 2 מוחלף בידי DFT גדול יותר (שבעצמו ניתן לחישוב בעזרת FFT).

דיאגרמת זרימת נתונים עבור N =8. אלגוריתם ה-FFT מסוג radix-2 מפרק DFT באורך N לשני DFT באורך N/2 ולאחר מכן משלב מספר DFT באורך 2 הנקראים פעולות "פרפר" (מקור השם בצורה של דיאגרמות זרימת הנתונים).

תהליך זה מדגים את הטכניקה הכללית של אלגוריתמי הפרד ומשול ; עם זאת, ביישומים רבים, נמנעת הרקורסיה המפורשת, ובמקום זאת עוברים על העץ החישובי באלגוריתם חיפוש לרוחב .

הביטוי-מחדש של DFT בגודל N כשני DFT בגודל N/2, כמפורט לעיל, נקרא לעיתים הלמה של Danielson–Lanczos, מאחר ששני המחברים כתבו על זהות זו ב-1942 (בהשפעת עבודתו של Runge משנת 1903 ). הם יישמו את הלמה שלהם כרקורסיה "לאחור", כך שגודל ה-DFT הוכפל שוב ושוב עד שספקטרום ההתמרה התכנס (הגם שנראה שהם לא הבינו שהם השיגו את הסיבוכיות האסימפטוטית O(NlogN)). עבודתם קדמה לזמינות נרחבת של מחשבים מכניים או אלקטרוניים והצריכה חישוב ידני (אפשר שעם עזרים מכניים כגון מכונת חיבור); הם דיווחו על זמן חישוב של 140 דקות עבור DFT בגודל 64 הפועל על קלט ממשי בדיוק של 3–5 ספרות. המאמר של קולי-טוקי משנת 1965 דיווח על זמן ריצה של 0.02 דקות עבור DFT מרוכב בגודל 2048 בריצה על IBM 7094 (כנראה תוך שימוש בנקודה צפה של 36 סיביות, כ-8 ספרות). בשינוי קנה המידה של הזמן לפי מספר הפעולות, דיווח זה מתאים להכפלת מהירות בערך פי 800,000. (כדי לתת פרספקטיבה לזמן החישוב הידני, 140 דקות עבור גודל 64 תואמות לממוצע של לא יותר מ-16 שניות לכל פעולת נקודה צפה, כש-20% מהן הן פעולות כפל. )

יישומי FFT בעלי ביצועים גבוהים עורכים שינויים רבים במימוש האלגוריתם בהשוואה למימוש הבסיסי הפשוט. לדוגמה, אפשר להשתמש במקרה בסיס גדול מ- N =1 להפחתת התקורה של הרקורסיה, גורם הסיבוב (twiddle factor) $\exp[-2\pi ik/N]$ ניתן לחישוב מראש, ולעיתים קרובות נעשה שימוש בבסיסים גדולים יותר משיקולי זיכרון מטמון ; אופטימיזציות אלו ואחרות יכולות ביחד לשפר את הביצועים בסדר גודל או יותר. (במימושים של האלגוריתם בספרי לימוד רבים, הרקורסיה של האלגוריתם חיפוש לעומק נשמטת לטובת הגישה הלא-רקורסיבית של אלגוריתם חיפוש לרוחב, אם כי נטען כי לרקורסיה של חיפוש לעומק ישנה מקומיות זיכרון טובה יותר. )

הרעיון

ניתן לראות את הצעד הבסיסי של קולי-טוקי FFT כפירוש מחדש של DFT חד ממדי כדומה ל-DFT דו ממדי. מערך הקלט החד ממדי באורך N=N₁N₂ מפורש מחדש כמטריצת N₁×N₂ המאוחסנת בסדר עמודה ראשי. אנו מחשבים DFT קטנים יותר לאורך כיוון N₂ (הכיוון הלא רציף), ואז מכפילים בגורמי סיבוב (Twiddle factors), ולבסוף מחשבים DFT לאורך כיוון N₁. ניתן לבצע את שלב השחלוף (transposition) באמצע, כמוצג כאן, או בהתחלה או בסוף. התהליך הזה הוא רקורסיבי.

איור של סדר עיקרי בשורה ובעמודה

כללית, אלגוריתמים של קולי-טוקי מבטאים-מחדש באופן רקורסיבי DFT בגודל מורכב N=N₁N₂ כ:

1. חשב N₁ התמרות DFT בגודל N₂.
2. הכפל בגורמי סיבוב.
3. חשב N₂ התמרות DFT בגודל N₁.

לרוב, N1 או N2 הוא גורם קטן (לא בהכרח ראשוני), הנקרא בסיס והוא יכול להיות שונה בין שלבי הרקורסיה. גורמי הסיבוב הם שורשי יחידה מרוכבים, שהם נקודות על מעגל היחידה. דהיינו, הכפלתם במספר מרוכב אחר, משנה רק את הזווית שלו. גורמי הסיבוב במקרה של DFT בגודל 2 הם 1-/+. ה-DFT הקטן של הבסיס נקרא לעיתים פרפר, בשל צורת תרשים זרימת הנתונים במקרה של radix-2.

וריאציות

ישנן וריאציות רבות של אלגוריתם קולי-טוקי. יישומי האלגוריתם עם בסיסים מעורבים מטפלים בגדלים מורכבים שאפשר לפרק לכמה גורמים (לרוב קטנים) בנוסף ל-2, בדרך כלל (אך לא תמיד) תוך שימוש באלגוריתם הישיר, בסיבוכיות O(N ²), עבור מקרי הבסיס העיקריים של הרקורסיה ( אפשר גם להשתמש באלגוריתם O(NlogN) עבור מקרי הבסיס העיקריים, כגון האלגוריתמים של Rader או Bluestein ) . בסיס מפוצל ממזג בסיסים 2 ו-4, תוך ניצול העובדה שההתמרה הראשונה של בסיס 2 אינה דורשת גורם סיבוב, על מנת להשיג את מה שהיה זמן רב מספר הפעולות האריתמטיות הנמוך ביותר הידוע עבור גדלים שהם חזקה של 2, למרות שווריאציות שפותחו לאחרונה השיגו מספר פעולות נמוך עוד יותר. (במחשבים של ימינו, הביצועים נקבעים יותר על ידי שיקולי זיכרון מטמון וצינורות עיבוד נתונים מאשר על ידי ספירת פעולות קפדנית; יישומי FFT אופטימליים משתמשים לעיתים קרובות בבסיסים גדולים יותר ו/או התמרות בסיסיות מקודדות בגודל משמעותי. ).

אפשר גם לראות את האלגוריתם של קולי-טוקי כביטוי מחדש של DFT חד-ממדי בגודל N כ-DFT דו ממדי בגודל N₁ על N₂ (בתוספת סיבובים), כאשר מטריצת הפלט עוברת שחלוף. עבור אלגוריתם radix-2, התוצאה נטו של השחלופים תואמת להיפוך סיביות של האינדקס של הקלט או הפלט. אם, במקום להשתמש בבסיס קטן, נשתמש בבסיס של בערך השורש הריבועי של N, ובשחלופים מפורשים של מטריצת קלט/פלט, נקבל אלגוריתם שנקרא אלגוריתם בן ארבעה צעדים (או שישה צעדים, תלוי במספר השחלופים). האלגוריתם הוצע בתחילה כדי לשפר את מקומיוּת הזיכרון, למשל עבור אופטימיזציה של זיכרון מטמון או פעולות מחוץ לזיכרון, ומאוחר יותר הוכח כי הוא אופטימלי כ-cache-oblivious algorithm .

הפירוק הכללי של קולי-טוקי משכתב את האינדקסים k ו-n כ$k=N_{2}k_{1}+k_{2}$ ו $n=N_{1}n_{2}+n_{1}$, בהתאמה, כאשר האינדקסים k_a ו-n_a הם בתחום שבין 0 ו-N_a-1 (כאשר a הוא 1 או 2). כלומר, הפירוק מארגן מחדש את הקלט (n) ואת הפלט (k) כמערכים דו-ממדיים של N₁ על N₂ בסדר עמודה ושורה, בהתאמה; ההבדל בין האינדקסים הללו הוא שחלוף, כמוזכר לעיל. כאשר האינדקסים החדשים האלה מוחלפים בנוסחת DFT עבור nk, הביטוי הצולב $N_{1}n_{2}N_{2}k_{1}$ נעלם (המעריך שלו הוא 1), ושאר הביטויים נותנים

$X_{N_{2}k_{1}+k_{2}}=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{1}N_{2}}}\cdot (N_{1}n_{2}+n_{1})\cdot (N_{2}k_{1}+k_{2})}$

$=\sum _{n_{1}=0}^{N_{1}-1}\left[e^{-{\frac {2\pi i}{N_{1}N_{2}}}n_{1}k_{2}}\right]\left(\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}$

$=\sum _{n_{1}=0}^{N_{1}-1}\left(\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N_{1}N_{2}}}n_{1}(N_{2}k_{1}+k_{2})}$ .

כאשר כל סכום פנימי הוא DFT בגודל N₂, כל סכום חיצוני הוא DFT בגודל N₁, והביטוי [...] בסוגריים הוא גורם הסיבוב.

אפשר להשתמש בבסיס r שרירותי (כמו גם בסיסים מעורבים), כפי שהראו הן קולי וטוקי והן גאוס (שנתן דוגמאות של שלבים של בסיס-3 ושל בסיס-6). קולי וטוקי הניחו במקור שפרפר הבסיס עבור r מחייב סיבוכיות חישובית של O(r²) ומכאן חישבו שהסיבוכיות עבור בסיס r היא O(r² N/r log_rN) = O(N log₂(N) r/log₂r); מחישוב ערכי r/log₂r עבור ערכי r בין 2 ל-12, נמצא שהבסיס האופטימלי הוא 3 (המספר השלם הקרוב ביותר ל-e, אשר ממזער את הביטוי r/log₂r ). אולם ניתוח זה היה שגוי: פרפר הבסיס הוא גם DFT וניתן לחישוב באמצעות אלגוריתם FFT בסיבוכיות O(rlogr), ומכאן שבסיס מבטל למעשה את הסיבוכיות O(r log(r) N/r log_rN), וה- r האופטימלי נקבע על סמך שיקולים מורכבים יותר. בפועל, r גדול למדי (32 או 64) חשוב על מנת לנצל למשל בצורה יעילה את המספר הגדול של האוגרים במעבדים מודרניים, ואפילו בסיס לא מוגבל r = √ N משיג גם הוא סיבוכיות של O(N log N) ויש לו יתרונות תאורטיים ומעשיים עבור N גדול.

הערות שוליים

הערך באדיבות ויקיפדיה העברית, קרדיט,
רשימת התורמים
רישיון cc-by-sa 3.0

35138013אלגוריתם ה-FFT של קולי-טוקי