Fast Fourier Transform
October 14, 2012 Leave a comment
Fast Fourier Transform adalah suatu algoritma yang digunakan untuk merepresentasikan sinyal dalam domain waktu diskrit dan domain frekuensi. DFT merupakan metode transformasi matematis untuk sinyal waktu diskrit ke dalam domain frekuensi. Secara sederhana dapat dikatakan bahwa DFT merupakan metode transformasi matematis sinyal waktu diskrit, sementara FFT adalah algoritma yang digunakan untuk melakukan transformasi tersebut.
Dari gambar di atas, maka didapat pola angka bahwa angka urutan dari input awal adalah bit reverse-nya.
Misal dari penjelasan di atas pada fig 7.6 kita ambil contoh FFT delapan masukan maka akan terbentuk persamaan-persamaan.
Missal y adalah variable output dan x variable input. m adalah perpangkatan dan n adalah urutan data berdasarkan pola pada fig 7.6.
Dari pola perhitungan butterfy dari gambar di atas dan maka jika diperinci perhitungannya akan menjadi sebagai berikut.
Kita sudah mendapat persamaan untuk keseluruhan algoritma FFT-nya. Akan tetapi bagaimana kita sekarang mendapatkan pola bilangan n yang tepat. Sekarang kita misalkan lagi pola bilangan n sebagai berikut agar kita bisa mecari algoritma pola n ini. Agar mendapat gambaran yang lebih luas kita buat table pola n dengan m sampai 3
Dari table di atas mari kita buat suatu algoritma pembentuk himpunan n. Dalam kode matlab berikut algoritma n.
pangkat = 3; N = 2^pangkat; n = zeros(1,(N/2)); for level2=1:pangkat intervaln = 2^(level2-1); level1=1; countbil = 0; while (level1<=(N/2)) for level3=1:intervaln countbil=countbil+1; n(level1)= countbil; level1=level1+1; end for level3=1:intervaln countbil=countbil+1; end end m = level2-1; for level4=1:(N/2) n(level4) end end
sekarang sudah kita dapatkan algoritma pembentuk himpunan n. Selanjutnya langsung saja dibuat kode untuk keseluruhan proses FFT.
%--------------------------bit-reverse bilreverse = zeros(1,2^pangkat); for val = 1:(2^pangkat) outval = 0; for nmak=1:pangkat andprocc=bitand((val-1),2^(nmak-1)); if (andprocc > 0) outval = outval + (2^(pangkat-nmak)); end end bilreverse(val) = outval + 1; end %--------------------------susun yin = zeros(1,N); for s=1:N yin(s) = y(bilreverse(s)); end %--------------------------jum2 yjum = zeros(1,N); for s=2:2:N yjum(s-1) = yin(s-1) + yin(s); yjum(s) = yin(s-1) - yin(s); end %---------------------------n genenator n = zeros(1,(N/2)); yjumreal = zeros(1,N); yjumim = zeros(1,N); ybuffreal = zeros(1,N); ybuffim = zeros(1,N); yakhir = zeros(1,N); ybuffreal = yjum; for level2=2:pangkat intervaln = 2^(level2-1); level1=1; countbil = 0; while (level1<=(N/2)) for level3=1:intervaln countbil=countbil+1; n(level1)= countbil; level1=level1+1; end for level3=1:intervaln countbil=countbil+1; end end for level4=1:(N/2) %(n(level4)) n %level2 pangkat NN = 2^(level2); inin = n(level4); inim = inin+(2^(level2-1)); sudutinin = (6.2832*(inin-1))/NN; sudutinim = (6.2832*(inim-1))/NN; yjumreal(inin) = ybuffreal(inin) + ( ybuffreal(inim) * cos(sudutinin)) + (ybuffim(inim) * sin(sudutinin)); yjumreal(inim) = ybuffreal(inin) + ( ybuffreal(inim) * cos(sudutinim)) + (ybuffim(inim) * sin(sudutinim)); yjumim(inin) = ybuffim(inin) + ( ybuffim(inim) * cos(sudutinin)) + (ybuffreal(inim) * ((-1)*sin(sudutinin))); yjumim(inim) = ybuffim(inin) + ( ybuffim(inim) * cos(sudutinim)) + (ybuffreal(inim) * ((-1)*sin(sudutinim))); end for level5=1:N ybuffreal(level5) = yjumreal(level5); ybuffim(level5) = yjumim(level5); end end for level5=1:N yakhir(level5) = ((yjumreal(level5)^2)+(yjumim(level5)^2)) /N; end
untuk menguji kode tersebut maka kita buat simulasi dengan masukan yang kita buat. Dan bandingkan dengan kode fungsi matlab. Missal kita coba dengan input FFT-8
Kode fungsi matlab.
t = 1:600; y = sind(t*256)+sind(t*512); Y = fft(y,8); Pyy = Y.* conj(Y) / 8; f = 1000*(0:4)/8; stem(f,Pyy(1:5)) title('Frequency content of y') xlabel('frequency (Hz)')
sekarang mari kita uji apakah sama hasilnya dengan kode di atas.
t = 1:600; y = sind(t*256)+sind(t*512); pangkat = 3; N = 2^pangkat; buffpiN = 6.2832/N; %--------------------------bit-reverse bilreverse = zeros(1,2^pangkat); for val = 1:(2^pangkat) outval = 0; for nmak=1:pangkat andprocc=bitand((val-1),2^(nmak-1)); if (andprocc > 0) outval = outval + (2^(pangkat-nmak)); end end bilreverse(val) = outval + 1; end %--------------------------susun yin = zeros(1,N); for s=1:N yin(s) = y(bilreverse(s)); end %--------------------------jum2 yjum = zeros(1,N); for s=2:2:N yjum(s-1) = yin(s-1) + yin(s); yjum(s) = yin(s-1) - yin(s); end %---------------------------n genenator n = zeros(1,(N/2)); yjumreal = zeros(1,N); yjumim = zeros(1,N); ybuffreal = zeros(1,N); ybuffim = zeros(1,N); yakhir = zeros(1,N); ybuffreal = yjum; for level2=2:pangkat intervaln = 2^(level2-1); level1=1; countbil = 0; while (level1<=(N/2)) for level3=1:intervaln countbil=countbil+1; n(level1)= countbil; level1=level1+1; end for level3=1:intervaln countbil=countbil+1; end end for level4=1:(N/2) %(n(level4)) n %level2 pangkat NN = 2^(level2); inin = n(level4); inim = inin+(2^(level2-1)); sudutinin = (6.2832*(inin-1))/NN; sudutinim = (6.2832*(inim-1))/NN; yjumreal(inin) = ybuffreal(inin) + ( ybuffreal(inim) * cos(sudutinin)) + (ybuffim(inim) * sin(sudutinin)); yjumreal(inim) = ybuffreal(inin) + ( ybuffreal(inim) * cos(sudutinim)) + (ybuffim(inim) * sin(sudutinim)); yjumim(inin) = ybuffim(inin) + ( ybuffim(inim) * cos(sudutinin)) + (ybuffreal(inim) * ((-1)*sin(sudutinin))); yjumim(inim) = ybuffim(inin) + ( ybuffim(inim) * cos(sudutinim)) + (ybuffreal(inim) * ((-1)*sin(sudutinim))); end for level5=1:N ybuffreal(level5) = yjumreal(level5); ybuffim(level5) = yjumim(level5); end end for level5=1:N yakhir(level5) = ((yjumreal(level5)^2)+(yjumim(level5)^2)) /N; end f = 1000*(0:(N/2))/N; stem(f,yakhir(1:((N/2)+1))) title('Frequency content of y') xlabel('frequency (Hz)')
dan kita bandingkan juga dari segi waktu eksekusi dengan contoh kodingan DFT.
% DFT "FIX 3 V" t = 1:600; y = sind(t*256)+sind(t*512); N = 8; x_real = zeros(1,N); x_im = zeros(1,N); x_val = zeros(1,(N/2)+1); buffpiN = 6.2832/N; for r=1:N/2 buff_fasa_r = buffpiN*(r-1); for k=1:N buff_fasa_k = buff_fasa_r*(k-1); x_real(r) = x_real(r) + ( y(k) * cos(buff_fasa_k) ); x_im(r) = x_im(r) + ( y(k) * sin(buff_fasa_k) ); end x_val(r) = ((x_real(r)^2)+(x_im(r)^2))/N; end f = 1000*(0:(N/2))/N; stem(f,x_val(1:((N/2)+1))) title('Frequency content of y') xlabel('frequency (Hz)')
dari waktu eksekusi FFT lebih cepat daripada DFT.
Berikut kode FFT dalam bahasa C.
//=================== INIT ==================== #define pangkat 3 //pangkat input #define size_pangkat 8 //2^pangkat_val //=================== INPUT ==================== float y[size_pangkat]; //=================== BIT REVERSE ==================== unsigned int bit_reverse[size_pangkat]; void func_bit_reverse() { unsigned int val,outval,nmak,andprocc; for(val=0;val<size_pangkat;val++) { outval = 0; for(nmak=0;nmak 0) outval = outval + (unsigned int)(pow(2,pangkat-nmak-1)); } bit_reverse[val] = outval + 1; } } //=================== END BIT REVERSE ==================== //=================== SUSUN ============================== float yin[size_pangkat]; void susun() { unsigned char s; unsigned int buff_bit; for(s=0;s<size_pangkat;s++) { buff_bit = bit_reverse[s]-1; yin[s]=y[buff_bit]; } } //=================== END SUSUN ============================== //=================== JUM2 ============================== float yjum[size_pangkat]; void jum2() { unsigned char s; for(s=1;s<size_pangkat;s+=2) { yjum[s-1] = yin[s-1] + yin[s]; yjum[s] = yin[s-1] - yin[s]; } } //=================== END JUM2 ============================== //=================== N GENERATOR ============================== float yjumreal[size_pangkat], yjumim[size_pangkat], ybuffreal[size_pangkat], ybuffim[size_pangkat], yakhir[size_pangkat]; unsigned int n[size_pangkat/2]; unsigned char level2; unsigned int intervaln, level1, countbil, level3, level4, level5, inin, inim; float NN, sudutinin, sudutinim; void n_generator() { for(level5=0;level5<size_pangkat;level5++) ybuffreal[level5] = yjum[level5]; for(level2=2;level2<=pangkat;level2++) { intervaln = (unsigned int)pow(2,(level2-1)); level1 = 1; countbil = 0; while(level1<=(size_pangkat/2)) { for(level3=1;level3<=intervaln;level3++) { countbil=countbil+1; n[level1-1]= countbil; level1++; } for(level3=1;level3<=intervaln;level3++) countbil++; } for(level4=1;level4<=(size_pangkat/2);level4++) { NN = pow(2,(float)level2); inin = n[level4-1]; inim = inin + (unsigned int)pow(2,((float)level2-1)); sudutinin = (6.2832*((float)inin-1))/NN; sudutinim = (6.2832*((float)inim-1))/NN; yjumreal[inin-1] = ybuffreal[inin-1] + ( ybuffreal[inim-1] * cos(sudutinin)) + (ybuffim[inim-1] * sin(sudutinin)); yjumreal[inim-1] = ybuffreal[inin-1] + ( ybuffreal[inim-1] * cos(sudutinim)) + (ybuffim[inim-1] * sin(sudutinim)); yjumim[inin-1] = ybuffim[inin-1] + ( ybuffim[inim-1] * cos(sudutinin)) + (ybuffreal[inim-1] * ((-1)*sin(sudutinin))); yjumim[inim-1] = ybuffim[inin-1] + ( ybuffim[inim-1] * cos(sudutinim)) + (ybuffreal[inim-1] * ((-1)*sin(sudutinim))); } for(level5=0;level5<size_pangkat;level5++) { ybuffreal[level5] = yjumreal[level5]; ybuffim[level5] = yjumim[level5]; } } for(level5=0;level5<size_pangkat;level5++) yakhir[level5] = (pow(yjumreal[level5],2) + pow(yjumim[level5],2))/(float)size_pangkat; } //=================== END N GENERATOR ============================== //=================== FFT FULL ============================== void fft_full() { func_bit_reverse(); susun(); jum2(); n_generator(); } //=================== END FFT FULL ==============================
Adapun kode dalam C yang saya uji dengan program CodeVision AVR dan Proteus dijalankan pada ATMega128 dengan input awal sudah didefinisikan.
/***************************************************** Chip type : ATmega128 Program type : Application Clock frequency : 11.059200 MHz Memory model : Small External SRAM size : 0 Data Stack size : 1024 *****************************************************/ #include #include #include #define RXB8 1 #define TXB8 0 #define UPE 2 #define OVR 3 #define FE 4 #define UDRE 5 #define RXC 7 #define FRAMING_ERROR (1<<FE) #define PARITY_ERROR (1<<UPE) #define DATA_OVERRUN (1<<OVR) #define DATA_REGISTER_EMPTY (1<<UDRE) #define RX_COMPLETE (1<<RXC) // USART0 Receiver buffer #define RX_BUFFER_SIZE0 8 char rx_buffer0[RX_BUFFER_SIZE0]; #if RX_BUFFER_SIZE0<256 unsigned char rx_wr_index0,rx_rd_index0,rx_counter0; #else unsigned int rx_wr_index0,rx_rd_index0,rx_counter0; #endif // This flag is set on USART0 Receiver buffer overflow bit rx_buffer_overflow0; // USART0 Receiver interrupt service routine interrupt [USART0_RXC] void usart0_rx_isr(void) { char status,data; status=UCSR0A; data=UDR0; if ((status & (FRAMING_ERROR | PARITY_ERROR | DATA_OVERRUN))==0) { rx_buffer0[rx_wr_index0]=data; if (++rx_wr_index0 == RX_BUFFER_SIZE0) rx_wr_index0=0; if (++rx_counter0 == RX_BUFFER_SIZE0) { rx_counter0=0; rx_buffer_overflow0=1; }; }; } #ifndef _DEBUG_TERMINAL_IO_ // Get a character from the USART0 Receiver buffer #define _ALTERNATE_GETCHAR_ #pragma used+ char getchar(void) { char data; while (rx_counter0==0); data=rx_buffer0[rx_rd_index0]; if (++rx_rd_index0 == RX_BUFFER_SIZE0) rx_rd_index0=0; #asm("cli") --rx_counter0; #asm("sei") return data; } #pragma used- #endif // USART0 Transmitter buffer #define TX_BUFFER_SIZE0 8 char tx_buffer0[TX_BUFFER_SIZE0]; #if TX_BUFFER_SIZE0<256 unsigned char tx_wr_index0,tx_rd_index0,tx_counter0; #else unsigned int tx_wr_index0,tx_rd_index0,tx_counter0; #endif // USART0 Transmitter interrupt service routine interrupt [USART0_TXC] void usart0_tx_isr(void) { if (tx_counter0) { --tx_counter0; UDR0=tx_buffer0[tx_rd_index0]; if (++tx_rd_index0 == TX_BUFFER_SIZE0) tx_rd_index0=0; }; } #ifndef _DEBUG_TERMINAL_IO_ // Write a character to the USART0 Transmitter buffer #define _ALTERNATE_PUTCHAR_ #pragma used+ void putchar(char c) { while (tx_counter0 == TX_BUFFER_SIZE0); #asm("cli") if (tx_counter0 || ((UCSR0A & DATA_REGISTER_EMPTY)==0)) { tx_buffer0[tx_wr_index0]=c; if (++tx_wr_index0 == TX_BUFFER_SIZE0) tx_wr_index0=0; ++tx_counter0; } else UDR0=c; #asm("sei") } #pragma used- #endif // Standard Input/Output functions #include unsigned char kata[32],i; // Alphanumeric LCD Module functions #asm .equ __lcd_port=0x1B ;PORTA #endasm #include //=================== INIT ==================== #define pangkat 3 //pangkat input #define size_pangkat 8 //2^pangkat_val //=================== INPUT ==================== float y[size_pangkat]; //=================== BIT REVERSE ==================== unsigned int bit_reverse[size_pangkat]; void func_bit_reverse() { unsigned int val,outval,nmak,andprocc; for(val=0;val<size_pangkat;val++) { outval = 0; for(nmak=0;nmak 0) outval = outval + (unsigned int)(pow(2,pangkat-nmak-1)); } bit_reverse[val] = outval + 1; } } //=================== END BIT REVERSE ==================== //=================== SUSUN ============================== float yin[size_pangkat]; void susun() { unsigned char s; unsigned int buff_bit; for(s=0;s<size_pangkat;s++) { buff_bit = bit_reverse[s]-1; yin[s]=y[buff_bit]; } } //=================== END SUSUN ============================== //=================== JUM2 ============================== float yjum[size_pangkat]; void jum2() { unsigned char s; for(s=1;s<size_pangkat;s+=2) { yjum[s-1] = yin[s-1] + yin[s]; yjum[s] = yin[s-1] - yin[s]; } } //=================== END JUM2 ============================== //=================== N GENERATOR ============================== float yjumreal[size_pangkat], yjumim[size_pangkat], ybuffreal[size_pangkat], ybuffim[size_pangkat], yakhir[size_pangkat]; unsigned int n[size_pangkat/2]; unsigned char level2; unsigned int intervaln, level1, countbil, level3, level4, level5, inin, inim; float NN, sudutinin, sudutinim; void n_generator() { for(level5=0;level5<size_pangkat;level5++) ybuffreal[level5] = yjum[level5]; for(level2=2;level2<=pangkat;level2++) { intervaln = (unsigned int)pow(2,(level2-1)); level1 = 1; countbil = 0; while(level1<=(size_pangkat/2)) { for(level3=1;level3<=intervaln;level3++) { countbil=countbil+1; n[level1-1]= countbil; level1++; } for(level3=1;level3<=intervaln;level3++) countbil++; } for(level4=1;level4<=(size_pangkat/2);level4++) { NN = pow(2,(float)level2); inin = n[level4-1]; inim = inin + (unsigned int)pow(2,((float)level2-1)); sudutinin = (6.2832*((float)inin-1))/NN; sudutinim = (6.2832*((float)inim-1))/NN; yjumreal[inin-1] = ybuffreal[inin-1] + ( ybuffreal[inim-1] * cos(sudutinin)) + (ybuffim[inim-1] * sin(sudutinin)); yjumreal[inim-1] = ybuffreal[inin-1] + ( ybuffreal[inim-1] * cos(sudutinim)) + (ybuffim[inim-1] * sin(sudutinim)); yjumim[inin-1] = ybuffim[inin-1] + ( ybuffim[inim-1] * cos(sudutinin)) + (ybuffreal[inim-1] * ((-1)*sin(sudutinin))); yjumim[inim-1] = ybuffim[inin-1] + ( ybuffim[inim-1] * cos(sudutinim)) + (ybuffreal[inim-1] * ((-1)*sin(sudutinim))); } for(level5=0;level5<size_pangkat;level5++) { ybuffreal[level5] = yjumreal[level5]; ybuffim[level5] = yjumim[level5]; } } for(level5=0;level5<size_pangkat;level5++) yakhir[level5] = (pow(yjumreal[level5],2) + pow(yjumim[level5],2))/(float)size_pangkat; } //=================== END N GENERATOR ============================== //=================== FFT FULL ============================== void fft_full() { func_bit_reverse(); susun(); jum2(); n_generator(); } //=================== END FFT FULL ============================== void main(void) { // Declare your local variables here // Input/Output Ports initialization // Port A initialization // Func7=In Func6=In Func5=In Func4=In Func3=In Func2=In Func1=In Func0=In // State7=T State6=T State5=T State4=T State3=T State2=T State1=T State0=T PORTA=0x00; DDRA=0x00; // Port B initialization // Func7=In Func6=In Func5=In Func4=In Func3=In Func2=In Func1=In Func0=In // State7=T State6=T State5=T State4=T State3=T State2=T State1=T State0=T PORTB=0x00; DDRB=0x00; // Port C initialization // Func7=In Func6=In Func5=In Func4=In Func3=In Func2=In Func1=In Func0=In // State7=T State6=T State5=T State4=T State3=T State2=T State1=T State0=T PORTC=0x00; DDRC=0x00; // Port D initialization // Func7=In Func6=In Func5=In Func4=In Func3=In Func2=In Func1=In Func0=In // State7=T State6=T State5=T State4=T State3=T State2=T State1=T State0=T PORTD=0x00; DDRD=0x00; // Port E initialization // Func7=In Func6=In Func5=In Func4=In Func3=In Func2=In Func1=In Func0=In // State7=T State6=T State5=T State4=T State3=T State2=T State1=T State0=T PORTE=0x00; DDRE=0x00; // Port F initialization // Func7=In Func6=In Func5=In Func4=In Func3=In Func2=In Func1=In Func0=In // State7=T State6=T State5=T State4=T State3=T State2=T State1=T State0=T PORTF=0x00; DDRF=0x00; // Port G initialization // Func4=In Func3=In Func2=In Func1=In Func0=In // State4=T State3=T State2=T State1=T State0=T PORTG=0x00; DDRG=0x00; // Timer/Counter 0 initialization // Clock source: System Clock // Clock value: Timer 0 Stopped // Mode: Normal top=FFh // OC0 output: Disconnected ASSR=0x00; TCCR0=0x00; TCNT0=0x00; OCR0=0x00; // Timer/Counter 1 initialization // Clock source: System Clock // Clock value: Timer 1 Stopped // Mode: Normal top=FFFFh // OC1A output: Discon. // OC1B output: Discon. // OC1C output: Discon. // Noise Canceler: Off // Input Capture on Falling Edge // Timer 1 Overflow Interrupt: Off // Input Capture Interrupt: Off // Compare A Match Interrupt: Off // Compare B Match Interrupt: Off // Compare C Match Interrupt: Off TCCR1A=0x00; TCCR1B=0x00; TCNT1H=0x00; TCNT1L=0x00; ICR1H=0x00; ICR1L=0x00; OCR1AH=0x00; OCR1AL=0x00; OCR1BH=0x00; OCR1BL=0x00; OCR1CH=0x00; OCR1CL=0x00; // Timer/Counter 2 initialization // Clock source: System Clock // Clock value: Timer 2 Stopped // Mode: Normal top=FFh // OC2 output: Disconnected TCCR2=0x00; TCNT2=0x00; OCR2=0x00; // Timer/Counter 3 initialization // Clock source: System Clock // Clock value: Timer 3 Stopped // Mode: Normal top=FFFFh // Noise Canceler: Off // Input Capture on Falling Edge // OC3A output: Discon. // OC3B output: Discon. // OC3C output: Discon. // Timer 3 Overflow Interrupt: Off // Input Capture Interrupt: Off // Compare A Match Interrupt: Off // Compare B Match Interrupt: Off // Compare C Match Interrupt: Off TCCR3A=0x00; TCCR3B=0x00; TCNT3H=0x00; TCNT3L=0x00; ICR3H=0x00; ICR3L=0x00; OCR3AH=0x00; OCR3AL=0x00; OCR3BH=0x00; OCR3BL=0x00; OCR3CH=0x00; OCR3CL=0x00; // External Interrupt(s) initialization // INT0: Off // INT1: Off // INT2: Off // INT3: Off // INT4: Off // INT5: Off // INT6: Off // INT7: Off EICRA=0x00; EICRB=0x00; EIMSK=0x00; // Timer(s)/Counter(s) Interrupt(s) initialization TIMSK=0x00; ETIMSK=0x00; // USART0 initialization // Communication Parameters: 8 Data, 1 Stop, No Parity // USART0 Receiver: On // USART0 Transmitter: On // USART0 Mode: Asynchronous // USART0 Baud rate: 9600 UCSR0A=0x00; UCSR0B=0xD8; UCSR0C=0x06; UBRR0H=0x00; UBRR0L=0x47; // Analog Comparator initialization // Analog Comparator: Off // Analog Comparator Input Capture by Timer/Counter 1: Off ACSR=0x80; SFIOR=0x00; // Global enable interrupts #asm("sei") // LCD module initialization lcd_init(16); lcd_gotoxy(0,0); lcd_putsf("FFT"); y[0] = -0.5008; y[1] = -0.3596; y[2] = 1.7377; y[3] = -1.7562; y[4] = 0.3008; y[5] = 0.7866; y[6] = -0.4148; y[7] = -0.2325; fft_full(); for(i=0;i<8;i++) { printf("yak=%.4f \r\n",yakhir[i]); } while (1) { // Place your code here }; }
________________________________________________________
Daftar Pustaka
http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
Hayes, Monson H.. 1999. Schaum’s Outline of Theory and Problems of Digital Signal Processing. The McGraw-Hill Companies, Inc.