FFT หรือ Fast Fourier transform เป็นหนึ่งในเครื่องมือสำคัญที่สุดในโลกของการประมวลผลสัญญาณ ช่วยให้เราวิเคราะห์องค์ประกอบความถี่จากสัญญาณดิบในโดเมนเวลาได้อย่างรวดเร็ว ใช้ตรวจจับความถี่ได้เร็วดีเลย แต่กระบวนการนี้มันหนักมากกก กินทรัพยากรทั้งพลังประมวลผลและหน่วยความจำ (RAM) โดยเฉพาะบนบอร์ดไมโครคอนโทรลเลอร์ทรัพยากรจำกัดอย่าง Arduino ทำให้ใช้งานยากสำหรับงานแบบเรียลไทม์ที่ต้องการผล FFT หลายครั้งต่อวินาที
ในบทสอนนี้ เราจะพาน้องไปรู้จักกับโปรแกรมที่เร็วที่สุดสำหรับทำ FFT บน Arduino พี่ตั้งชื่อมันว่า ApproxFFT เพราะมีการประมาณค่า (approximation) ไว้หลายขั้นตอนมากก แต่ด้วยการประนีประนอมความแม่นยำเพียงนิดเดียว เราก็ได้ความเร็วเพิ่มขึ้นมาประมาณ 3 เท่าเลยทีเดียว
ถ้าน้องอยากรู้วิธีใช้แบบตรงๆ ข้ามไปที่ข้อ 5 ได้เลยจ้า
โครงการนี้พี่อธิบายไว้ในวิดีโอความยาวครึ่งชั่วโมง น้องจะเลือกดูวิดีโอหรืออ่านบทสอนนี้ก็ได้ตามสะดวก
ก่อนหน้านี้พี่เคยเตรียมโค้ดสองตัวสำหรับงานเดียวกันไว้: EasyFFT: โค้ดตัวนี้ให้ผลลัพธ์ FFT ที่แม่นยำ แต่มันค่อนข้างช้าและกินหน่วยความจำสูง แนะนำให้อ่านบทสอนตัวนี้เพื่อทำความเข้าใจพื้นฐาน FFT ก่อน QuickFFT: โค้ดตัวนี้ให้ผลลัพธ์เร็วมากกก แต่ค่าอัมปลิจูด (amplitude) เบี่ยงเบนไปไกลและใช้กับงานส่วนใหญ่ไม่ได้ มันยังให้พีค (peak) หลายจุดซึ่งอาจทำให้ผลลัพธ์เข้าใจผิดได้
1: ความจำเป็นของฟังก์ชัน ApproxFFT
ถ้าน้องดูรูปที่แนบมา จะเห็นชัดว่า QuickFFT มีข้อได้เปรียบสำคัญเหนือวิธีดั้งเดิม (คำนวณความเร็วด้วยฟังก์ชัน sine/cosine ในตัวของ EasyFFT) มันเร็วเกือบ 4 เท่า แต่ขาดความแม่นยำ อย่างที่เห็นในรูปที่สอง ค่าอัมปลิจูดเบี่ยงเบนไปไกล สาเหตุหลักคือเราใช้คลื่นสี่เหลี่ยม (ประมาณเป็นคลื่น sine !!!) ซึ่งประกอบด้วยหลายความถี่ (ฮาร์มอนิก) ทำให้มีหลายคลื่นปรากฏในผลลัพธ์ สร้างความยุ่งยากสำหรับการประยุกต์ใช้งานต่างๆ
ตรงนี้เราเคยประมาณคลื่น sine/cosine เป็นคลื่นสี่เหลี่ยม ถ้าเราใช้การประมาณที่ดีกว่านี้สำหรับคลื่นเหล่านี้ เราก็จะเข้าใกล้ผลลัพธ์ที่แท้จริงมากขึ้น ในโค้ดนี้เราได้พิจารณาการประมาณที่ดีกว่าสำหรับคลื่นเหล่านี้ และเรายังมีตัวเลือกในการควบคุมความแม่นยำของคลื่นเหล่านี้ด้วย การปรับแต่งการตั้งค่านี้ทำให้เราสามารถเล่นกับความสัมพันธ์ระหว่างความแม่นยำ/ความเร็วได้หลายแบบเลย สู้งานนะน้อง!
2: ฟังก์ชันไซน์และโคไซน์แบบเร็วปึ๊ก
ถ้าน้องย้อนกลับไปดูฟังก์ชัน EasyFFT ในตอนก่อนหน้า เรามีสองตัวเลือกให้ใช้ ตัวนึงคือแบบแม่นยำ ใช้ฟังก์ชันไซน์/โคไซน์ในตัวของ Arduino ซึ่งแม่นแต่ช้าครับ อีกวิธีคือใช้ตารางค้นหา (lookup table) ซึ่งเร็วขึ้นนิดหน่อยแต่ความแม่นยำก็ลดลงไปบ้าง
ในวิธีนี้เราใช้แนวคิดที่ต่างออกไปโดยสิ้นเชิงในการคำนวณค่า ซึ่งทำให้มันเร็วขึ้นได้ประมาณ 10-12 เท่าเมื่อเทียบกับวิธีปกติครับ ฟังก์ชันพวกนี้ใช้การดำเนินการบิตชิฟต์ (bitshift operation) อย่างหนักสำหรับการคูณและการหาร (แบบประมาณ) ซึ่งการดำเนินการแบบนี้เร็วกว่าการคูณหารปกติแบบสุดๆ แต่มันก็อาจทำให้ความแม่นยำลดลง ซึ่งเราต้องจัดการให้ดี
ในฟังก์ชัน FFT ทั่วไป ส่วนที่กินเวลามากที่สุดคือการคูณเลขทศนิยม (Float) ด้วยค่าไซน์/โคไซน์ของมุมต่างๆ ในขั้นตอนที่เรียกว่า "การคำนวณผีเสื้อ (Butterfly Calculation)" ตามโค้ดด้านล่าง ทั้งสองกระบวนการนี้ช้าและใช้เวลามาก
tr=c*out_r[i10+n1]-s*out_im[i10+n1];
ti=s*out_r[i10+n1]+c*out_im[i10+n1]; // c is cosine and s is sine of some value
ใน ApproxFFT นี้ พี่ใช้วิธีที่ต่างออกไปซึ่งเร็วขึ้น 10-12 เท่า โดยใช้ การดำเนินการบิตชิฟต์ (Bitshift Operations) ร่วมกับอัลกอริทึมการค้นหาแบบทวิภาค (Binary Search) เพื่อหาค่า $A \cdot \sin(\theta)$ โดยตรง โดยไม่ต้องทำการคูณเลขทศนิยมเลย
หลักการทำงานของ Binary Search Sine:
เราเก็บตาราง isin_data ซึ่งแมปมุมกับค่าไซน์ในรูปแบบจำนวนเต็ม (0-1024 แทน 0-360 องศา) และใช้วิธีการแบ่งครึ่งช่วง (halving range) เพื่อหาค่าแอมพลิจูด ที่นี่เราใช้สิ่งที่คล้ายกับอัลกอริทึมการค้นหาแบบทวิภาค มันจะให้ผลลัพธ์ประมาณสำหรับ A*sin(θ) โดยตรง และเราไม่ต้องทำการคูณ ในกราฟแรกจะแสดงพล็อตของ ArcsinX เทียบกับ X พล็อตที่คล้ายกันนี้แหละที่ถูกเก็บเป็นตารางใน Arduino โดยที่มุมถูกแมปจาก 0-1024 ซึ่งเทียบเท่ากับ 0-360 องศา ค่า arcsine ก็ถูกพิจารณาเป็นผลคูณของ 128 ด้วย
นี่คือขั้นตอนการคำนวณ (เราจะมองข้ามค่าลบและพิจารณาแค่ 0-90 องศา เช่น 500*sin 50 (หน่วยองศา)):
สำหรับค่าใดๆ (เช่น 500) ผลคูณของไซน์สามารถเป็นได้ตั้งแต่ 0 (ที่ 0 องศา) ถึงค่านั้น (ที่ 90 องศา) (เช่น 0-500) เราจะตรวจสอบค่ามุมที่จุดกึ่งกลาง เราจะหารค่าอินพุต (500) ลงครึ่งหนึ่งและตรวจสอบมุมสำหรับ 0.5 จากนั้นเราสามารถตรวจสอบได้ว่าผลลัพธ์จะอยู่ระหว่าง 0-จุดกึ่งกลาง หรือ จุดกึ่งกลาง-ค่าสูงสุด (ในกรณีของเรา มุมสำหรับ 0.5 จะเป็น 30 องศาซึ่งน้อยกว่า 50 ดังนั้นเราสามารถสรุปได้ว่าผลลัพธ์จะอยู่ที่ไหนสักแห่งระหว่าง 250-500) ณ จุดนี้ เราได้ลดขนาดของช่วงที่เป็นไปได้ลงครึ่งหนึ่งแล้ว (จากเดิม 0-500 เหลือ 250-500)
ในช่วงที่กำหนดใหม่ เราจะทำขั้นตอนที่คล้ายกันซ้ำอีกครั้ง ในตัวอย่างของเรา จุดกึ่งกลางจะอยู่ที่ 375 ซึ่งเราจะตรวจสอบมุมสำหรับ .75 นั่นจะได้ 48 องศา ดังนั้นคำตอบของเราจะอยู่ระหว่าง 375 ถึง 500
ต้องทำขั้นตอนที่คล้ายกันนี้ซ้ำหลายครั้งเพื่อให้ได้ผลลัพธ์ที่แม่นยำ อย่างที่แสดงในภาพที่สอง (กราฟ) ยิ่งเราเพิ่มระดับ (iteration) มากเท่าไหร่ ผลลัพธ์ก็จะยิ่งใกล้เคียงกับค่าจริงมากขึ้น โดยที่ไม่ต้องพึ่งพา FPU (Floating Point Unit) ของ Arduino เลยสักนิด
ที่นี่เราสามารถดำเนินการทั้งหมดได้ด้วยแค่การบิตชิฟต์ และเราก็กำจัดขั้นตอนการคูณเลขทศนิยมออกไปได้แล้ว งานนี้จัดไปวัยรุ่น! สู้งานนะน้อง
3: Fast Root of Squared Sum (FastRSS)
ฟังก์ชันนี้ทำการประมาณค่า RSS ครับ มันไม่ได้ช่วยให้โปรแกรมเร็วขึ้นแบบสุดๆ แต่ก็ประหยัดเวลาไปได้เป็นมิลลิวินาทีเหมือนกัน ดีกว่าไม่มีเนอะ
หลังจากได้ค่า Real กับ Imaginary มาจาก FFT แล้ว ขั้นตอนต่อไปคือหาขนาด (Magnitude) ของสัญญาณด้วยสูตร $\sqrt{R^2 + I^2}$ ซึ่งถ้าใช้ฟังก์ชัน sqrt() มาตรฐานเนี่ย มันช้ามากเลยครับพี่
กราฟที่แนบมาแสดงค่าความผิดพลาด ถ้าเราสมมติว่าคำตอบของรากที่สองของผลบวกกำลังสอง เท่ากับค่าที่มากกว่าเลย ความผิดพลาดจะลดลงเมื่ออัตราส่วนระหว่างค่าทั้งสองเพิ่มขึ้น ผมเลยสร้างฟังก์ชัน fastRSS ขึ้นมาเพื่อประมาณค่านี้โดยใช้ความสัมพันธ์ของอัตราส่วนระหว่างค่ามากและค่าน้อย ถ้าค่าหนึ่งมากกว่าอีกค่ามากๆ เราก็ใช้ค่าที่มากกว่าเป็นคำตอบได้เลยโดยมีความผิดพลาดนิดเดียว แต่ถ้าค่าทั้งสองใกล้เคียงกัน เราจะใช้ตาราง RSSdata และการปรับสเกลมาช่วยชดเชยค่าให้ใกล้เคียงความจริงมากที่สุด ซึ่งช่วยประหยัดเวลาไปได้หลายมิลลิวินาทีเลย ในโค้ดนี้ ถ้าอัตราส่วนเกินค่าที่กำหนดไว้ มันก็จะตอบเป็นค่าที่มากกว่าเลยเพื่อประหยัดเวลา แต่ถ้าอัตราส่วนต่ำกว่านั้น ก็จะใช้การปรับสเกลที่กำหนดไว้ล่วงหน้าตามค่าในตาราง RSSdata นั่นแหละ
4: โฟลว์การคำนวณโดยรวม:
โฟลว์การคำนวณทั้งหมดจะเป็นตามนี้ครับ โปรแกรมจะทำงานตามลำดับขั้นตอนนี้เพื่อประสิทธิภาพสูงสุด:
- ปรับสเกลข้อมูลอินพุตให้อยู่ในช่วง +512 ถึง -512: ขั้นตอนนี้สำคัญเพื่อรักษาความแม่นยำของข้อมูล เพราะเราใช้การดำเนินการบิตชิฟต์ ถ้าเป็นเลขคี่จะมีความแม่นยำหาย ยิ่งตัวเลขใหญ่ การสูญเสียน้อยลง
- สร้างลำดับบิตกลับ (Bit reversed order)
- เรียงลำดับข้อมูลอินพุตตามลำดับบิตกลับที่สร้างไว้
- ทำการแปลงฟูเรียร์ (Frequency transform) โดยใช้ fast sine และ fast cosine ในจุดที่จำเป็น หลังจากทุกลูป จะมีการตรวจสอบว่าค่าแอมพลิจูดเกิน 15000 หรือไม่ ถ้าเกิน ทั้งอาร์เรย์ Real และ Imaginary จะถูกหารด้วย 2 ทันที (เพื่อป้องกันการล้นของ integer 16-bit) และบันทึกค่าสเกลที่ใช้หารไว้
- คำนวณรากที่สองของผลบวกกำลังสองโดยใช้ฟังก์ชัน FastRSS
- ตรวจจับค่าสูงสุดและส่งคืนเป็นผลลัพธ์
5: วิธีใช้งาน
- ต้องประกาศข้อมูลลุคอัพ (Lookup data) เป็นตัวแปรโกลบอล แค่ก็อปปี้ข้อมูลด้านล่างไปวางไว้ด้านบนสุดของโค้ดเลย
//---------------------------------lookup data------------------------------------//
byte isin_data[128]=
{0, 1, 3, 4, 5, 6, 8, 9, 10, 11, 13, 14, 15, 17, 18, 19, 20,
22, 23, 24, 26, 27, 28, 29, 31, 32, 33, 35, 36, 37, 39, 40, 41, 42,
44, 45, 46, 48, 49, 50, 52, 53, 54, 56, 57, 59, 60, 61, 63, 64, 65,
67, 68, 70, 71, 72, 74, 75, 77, 78, 80, 81, 82, 84, 85, 87, 88, 90,
91, 93, 94, 96, 97, 99, 100, 102, 104, 105, 107, 108, 110, 112, 113, 115, 117,
118, 120, 122, 124, 125, 127, 129, 131, 133, 134, 136, 138, 140, 142, 144, 146, 148,
150, 152, 155, 157, 159, 161, 164, 166, 169, 171, 174, 176, 179, 182, 185, 188, 191,
195, 198, 202, 206, 210, 215, 221, 227, 236};
unsigned int Pow2[14]={1,2,4,8,16,32,64,128,256,512,1024,2048,4096};
byte RSSdata[20]={7,6,6,5,5,5,4,4,4,4,3,3,3,3,3,3,3,2,2,2};
//---------------------------------------------------------------------------------//
ฟังก์ชัน ApproxFFT, fast_sine, fast_cosine และ fastRSS ต้องก็อปปี้ไปวางไว้ท้ายสุดของโค้ด
วิธีเรียกใช้ฟังก์ชัน:
float f=Approx_FFT(data,sample,sampling_rate);
ฟังก์ชันนี้จะคืนค่าความถี่ที่มีแอมพลิจูดสูงสุดเป็นค่าเริ่มต้น เหมือนกับฟังก์ชัน EasyFFT และ QuickFFT เลย
อาร์กิวเมนต์แรกคืออาร์เรย์ของข้อมูลที่ต้องการทำ FFT
อาร์กิวเมนต์ที่สองคือจำนวนตัวอย่าง (sample) ที่ควรเป็น 2^n เช่น 2, 4, 8, 16, 32, 64, 128,... โดยจำนวนตัวอย่างสูงสุดจะถูกจำกัดโดยหน่วยความจำที่มีอยู่
อินพุตตัวที่สามคืออัตราการสุ่มตัวอย่าง (sampling rate)
- การตั้งค่าความแม่นยำ (Accuracy setting):
ค่าตัวนี้ปรับได้ตั้งแต่ 1 ถึง 7 ยิ่งตัวเลขสูงยิ่งแม่นยำมากขึ้น การเพิ่มความแม่นยำจะทำให้คลื่นไซน์ประมาณ (approximate sine wave) ของเราพอดีกับค่าจริงมากขึ้น โดยค่าเริ่มต้น (default) จะตั้งไว้ที่ 5 ซึ่งในกรณีส่วนใหญ่ค่านี้ก็ใช้ได้ดีอยู่แล้ว ผลลัพธ์และความเร็วที่อ้างอิงทั้งหมดก็มาจากการตั้งค่านี้แหละ ความเร็วและความแม่นยำสามารถปรับแต่งได้ด้วยการเปลี่ยนค่านี้
int fast_sine(int Amp, int th)
{
int temp3,m1,m2;
byte temp1,temp2, test,quad,accuracy;
accuracy=5; // set it value from 1 to 7, where 7 being most accurate but slowest
// accuracy value of 5 recommended for typical applicaiton
while(th>1024){th=th-1024;} // here 1024 = 2*pi or 360 deg
while(th<0){th=th+1024;}
.
.
.
- โค้ดนี้สามารถดัดแปลงเพื่อแสดง (และใช้) ผลลัพธ์ดิบ (Raw output) หรือแอมพลิจูดของความถี่ต่างๆ ได้ โดยผลลัพธ์จะแสดงเป็นผลคูณของ 2^n เนื่องจากตัวแปรจำนวนเต็ม (integer) เก็บค่าได้สูงสุดแค่ +-32000 ผลลัพธ์เลยถูกแสดงในรูปแบบผลคูณแทน
6: สรุป
การทดสอบทั้งหมดที่เห็นในภาพด้านบนทำบน Arduino mega ด้วยระดับความแม่นยำที่ 5 มาดูกันว่าผลทดสอบเจ๋งแค่ไหน:
จากการทดสอบบน Arduino Mega ที่ระดับความแม่นยำ 5 พบว่า:
- ความเร็วในการประมวลผล เร็วกว่า FFT ทั่วไปมากกว่า 3 เท่า งานไวปรื๊ด!
- การใช้หน่วยความจำลดลงเกือบ 50% ประหยัด RAM ไปได้เยอะ
- ผลลัพธ์ใกล้เคียงกับค่าจริงมาก (คลาดเคลื่อนต่ำ) เหมาะจะเอาไปใช้ในโปรเจคฝังตัว (embedded) จริงๆ
หวังว่าโค้ดนี้จะมีประโยชน์กับโปรเจคของน้องๆ นะครับ ถ้ามีคำถามหรือข้อเสนอแนะอะไร ก็คอมเมนต์มาคุยกันได้เลย