Multi-sine signal with Matlab


ผมมักมีปัญหาเรื่องการใ้ช้คำสั่งของ Matlab อยู่เนือง ๆ เนื่องจากนาน ๆ จะใช้คำนั้น ๆ ทีนึง ส่วนใหญ่เปิดคู่มือดูก็รู้เรื่องเข้าใจ มีอยู่คำสั่งหนึ่งใช้ทีไรปวดหัวทุกที นั่นคือคำสั่ง idinput เพื่อใช้ในการสร้างสัญญาณ Multi-sine ซึ่งสัญญาณนี้สามารถแสดงได้ดังสมการต่อไปนี้

u(k) = A\displaystyle{\sum^d_{i=1}\sin(\omega k+ \phi_i)}

ในการสั่งงาน Matlab เราต้องมั่นใจว่าย่านความถี่ของสัญญาณที่เราเลือกนั้นมันตรงตามที่เราต้องการหรือไม่ ซึ่งทำได้โดยการแสดง Power spectrum ของสัญญาณ $$u(k)$$ ตรงนี้ใช้ทีไรก็ลืมทุกที คำสั่ง Matlab นั้นสั่งดังนี้ครับ

N = 4000;         % จำนวนแซมเปิลของสัญญาณ
u = idinput(N,'sin',[0 0.2],[-1 1],[30 10 1]);

สำหรับตัวเลือกของ idinput นั้น [0 0.2] คือช่วงความถี่ที่ต้องการ [30 10 1] 30 คือจำนวนสัญญา sine ที่ต้องการ 10 คือ จำนวนของการลอง และ 1 คือ grid skip สำหรับรายละเอียดของเรื่องพวกนี้ไปดูในคู่มือ Matlab เอาเองนะครับ

ปัญหาที่จะพูดถึงคือเรื่องช่วงความถี่ ดูคำอธิบายของ Matlab ทีไรก็ปวดหัวทุกที เอาง่าย ๆ ถ้าต้องการสัญญาณความถี่สูงสุดที่ 50 Hz ความถี่ในการชักข้อมูล (Sampling) ก็ต้องเป็นสองเท่า นั่นคือ 100 Hz หรือทุก ๆ 0.01 วินาที ช่วงความถี่ที่ต้องการคือ [0 0.2] นั้นหมายความว่าความถี่สูงสุดคือ 50*0.2 หรือ 10 Hz นั่นเอง

ในทางปฏิบัติเราจะถูกกำหนดอัตราการชักข้อมูลมาหรือ 0.01 วินาที ถ้าต้องการความถี่สูงสุดของสัญญาณเป็น 30 Hz เราก็ต้องเอา 0.01 วินาทีไปหาความถี่ในการชักข้อมูลก่อนคือ 100 Hz แล้วหารด้วยของเป็น 50 Hz แล้วกำหนดช่วงความถี่ในที่นี้คือ [0 0.6] นั่นเอง

เราจะรู้ได้อย่างไรว่าความถี่ของสัญญาณ Multi-sine อยู่ในช่วงที่ถูกต้อง (ถ้าคำนวณตามข้างบนมันก็ต้องถูกเป็นธรรมดา) แต่เพื่อความมั่นใจในทฤษฏี ปฏิบัติก็ต้องได้ด้วย วิธีการพล็อต Spectrum นั้นก็คือการหาค่ากำลังสองของค่าเฉลี่ยของค่าสมบูรณ์ของค่า FFT ของสัญญาณนั่นเอง ซึ่งทำได้ดังนี้ (ใช้ตัวอย่างความถี่สูงสุด 30 Hz)

u = idinput(N,'sine',[0 0.6],[-1 1],[50 10 1]);
ts = 0.01;                       % sampling time
T = (N-1)*ts;                  % เวลาในหนึ่งคาบ หน่วยเป็นวินาที
                                    % นับจากศูนย์
p = abs(fft(u))/(N/2);       % ค่าเฉลี่ยของค่าสมบูรณ์ของสัมประสิทธิ์ FFT
                                    % ของ u คิดแค่ครึ่งคาบบวก
p = p(1:N/2).^2;             % ค่ากำลังสอง
freq = [0:N/2-1]/T;            % คำนวณค่าความถี่
semilogy(freq,p);              % พล็อตกราฟ semilog แกน x เป็นความถี่
xlabel('Hz');

ผลที่ได้แสดงได้ตามรูปข้างล่าง

หวังว่าคงได้ประโยชน์บ้างนะครับ

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s