สอบถามเรื่องการแปลงอนุกรมฟูเรีย (Fast Fourier transform)

 
คืออยากจะลองทำเครื่องวัดเสียงหน่ะครับ ว่าความถี่ไหนมันแรงมันน้อย
เจอะหลายๆที่เขาทำ ผมเลยอยากทำบ้าง แต่ติดที่เรื่อง Fast Fourier transform ครับ
คือผมไม่มีพื้นฐานเลยเรื่องนี้ จึงอยากขอข้อมูลจากพี่ๆบ้างหน่ะครับ

ความต้องการคือว่า ผมอ่านค่าอะนาล็อกโดยใช้ MCU แล้ว จากนั้นว่าจะให้ MCU มันเป็นตัวแปลง
จากโดเมนของเวลาเป็นโดเมนของความถี่(ไม่รู้เรียกถูกรึปล่าวนะ ยังไม่ได้เรียนหน่ะครับ)
เพื่อนเอามาแสดงที่ LCD ครับ ผมแค่อยากทราบว่า การคำนวน Fast Fourier transform นี่
เขาเอาข้อมูลที่เรา Sampling มาคำนวนยังไง แค่นี้หล่ะครับ นอกนั้นก็จะค่อยๆทำเอาครับ ขอบคุณมากครับ

ในโลกของอนาลอก การแปลงฟูเรียร์เป็นการแปลงฟังก์ชั่นใน time domain ( y = F(t) ) ให้ไปอยู่ใน frequency domain ( z = X(omega) ) โดย

     

แต่เมื่อใช้การ sampling ข้อมูลตามจังหวะเวลาสม่ำเสมอ จะได้ข้อมูลเป็นช่วงๆ จำนวน N ค่า การแปลงฟูเรียร์สำหรับข้อมูลลักษณะนี้จะเรียกว่า Discrete Fourier Transform (DFT) ซึ่งมีรูปแบบดังนี้

     

ซึ่งสามารถเขียนโปรแกรมเพื่อหาค่า X(omega) แต่ละค่าได้ แต่หน้าตาข้างบนจะเกิดจำนวน operation ที่แปรไปตามค่า N ยกกำลังสอง จึงมีผู้คิดหาวิธีแปลงให้เร็วกว่านี้ ที่เรียกกันทั่วๆไปว่า Fast Fourier Transform (FFT) ที่สามารถลดจำนวน operation ลงเหลือเพียง N log N เท่านั้น ซึ่งก็มีอยู่หลายวิธี แต่หลักการคล้ายๆกัน คือ ลด operation ที่ทำซ้ำๆ กัน และแยกจำนวนตัวอย่างที่ต้องการแปลงจาก N ออกเป็นกลุ่มย่อยๆ เช่น N = A x B แล้วทำการแปลงแต่ละกลุ่มก่อน แล้วจึงนำผลลัพธ์ของการแปลงแต่ละกลุ่มที่ได้มารวมกันเป็นผลลัพธ์ที่แท้จริงอีกที

ค่าของ A กับ B จะเป็นเท่าใดก็ได้ตราบเท่าที่ N = A x B ซึ่งวิธีการแบ่งกลุ่มอย่างหนึ่งที่ใช้กันอยู่ทั่วไปคือแบ่งเป็นกลุ่มเท่าๆกัน จากจำนวนตัวอย่าง N ค่า แบ่งเป็นกลุ่มละ N/2 จำนวน 2 กลุ่ม แล้วแต่ละกลุ่มแบ่งย่อยลงไปอีกเป็นกลุ่มละ N/4 ทำอย่างนี้ไปเรื่อยๆ จนเหลือกลุ่มละ 2 ค่า ซึ่งวิธีการนี้จำเป็นที่จำนวนข้อมูล N จะต้องเป็นเลขยกกำลังของ 2

ตัวอย่างโค้ดในภาษา C (จากเว็บ www.rosettacode.org)
#include <stdio.h>
#include <math.h>
#include <complex.h>

double PI;
typedef double complex cplx;

void _fft(cplx buf[], cplx out[], int n, int step)
{
  if (step < n) {
    _fft(out, buf, n, step * 2);
    _fft(out + step, buf + step, n, step * 2);

    for (int i = 0; i < n; i += 2 * step) {
      cplx t = cexp(-I * PI * i / n) * out[i + step];
      buf[i / 2]     = out[i] + t;
      buf[(i + n)/2] = out[i] - t;
    }
  }
}

void fft(cplx buf[], int n)
{
  cplx out[n];
  for (int i = 0; i < n; i++) out[i] = buf[i];

  _fft(buf, out, n, 1);
}

void show(const char * s) {
  printf("%s", s);
  for (int i = 0; i < 8; i++)
    if (!cimag(buf[i]))
      printf("%g ", creal(buf[i]));
    else
      printf("(%g, %g) ", creal(buf[i]), cimag(buf[i]));
}

int main()
{
  PI = atan2(1, 1) * 4;
  cplx buf[] = {1, 1, 1, 1, 0, 0, 0, 0};

  show("Data: ");
  fft(buf, 8);
  show("\nFFT : ");

  return 0;

}


คิดว่าต้องใช้ Fixed Point นะครับ

ขอบคุณคุณ pa_ul มากครับ
พอเป็นแนวทางได้แล้วครับ เดี่ยวผมจะพยายามศึกษาเพิ่ม หากติดขัดยังไง ก็คงต้องมารบกวนอีกที

สอบถามคุณ firmware.c หน่อยครับเจ้า Fix Point นี่เป็นแบบไหนหรอครับ
คือว่ายังไม่มีพื้นฐานเรื่องอนุกรมฟูเรียเลยครับ กำลังศึกษาอยู่ครับ ส่วนใหญ่ยังเป็นตำราหน่ะครับ
ยังไม่ได้เจอะคนที่เขาลองเอามาเขียน ที่เจอะก็เป็น Assembly ผมเลยเริ่มต้นช้าหน่อย

ขอบคุณครับ

ที่จริงสำหรับ DFT/FFT ถ้าตัดเรื่องความเร็วในการทำงานออกไปก่อน การเขียนโปรแกรมก็ไม่ได้ยุ่งยากอะไรมากมาย ค่อนข้างจะตรงไปตรงมาด้วยซ้ำไป แต่เรื่องที่ยุ่งยากซับซ้อนมากกว่านั้นก็คือการตีความและทำความเข้าใจกับผลการคำนวณที่ได้ออกมาว่ามีความหมายอย่างไร และเข้าใจถึงเหตุที่ส่งผลต่อความไม่สมบูรณ์ของผลลัพธ์ที่ได้ออกมา เช่น leakage และ aliasing ซึ่งตรงนี้ทำให้ต้องต้องศึกษาเพิ่มเติมไปถึงเรือง Nyquist sampling theory ด้วย

สำหรับตัวอย่างที่ให้ไว้ จะใช้เลขจินตภาพ (complex number) เพราะว่าการแปลงฟูเรียร์จะได้ผลลัพธ์เป็นเลขจินตภาพเสมอ เลขจินตภาพในทางคณิตศาสตร์จะสามารถมองเป็นเรื่องเฟสของความถี่ได้ (ดูเรื่อง Euler formula) ถ้าต้องการ magnitude ของสัญญาณมาใช้ในการแสดงผลก็สามารถคำนวณได้จากสมการ m = sqrt(x^2 + i^2)