The Fast Fourier Transform Chapitre 12 - Analog Devices - Revenir à l'accueil

Farnell Element 14 :

Farnell-NA555-NE555-..> 08-Sep-2014 07:33 1.5M

Farnell-AD9834-Rev-D..> 08-Sep-2014 07:32 1.2M

Farnell-MSP430F15x-M..> 08-Sep-2014 07:32 1.3M

Farnell-AD736-Rev-I-..> 08-Sep-2014 07:31 1.3M

Farnell-AD8307-Data-..> 08-Sep-2014 07:30 1.3M

Farnell-Single-Chip-..> 08-Sep-2014 07:30 1.5M

Farnell-Quadruple-2-..> 08-Sep-2014 07:29 1.5M

Farnell-ADE7758-Rev-..> 08-Sep-2014 07:28 1.7M

Farnell-MAX3221-Rev-..> 08-Sep-2014 07:28 1.8M

Farnell-USB-to-Seria..> 08-Sep-2014 07:27 2.0M

Farnell-AD8313-Analo..> 08-Sep-2014 07:26 2.0M

Farnell-SN54HC164-SN..> 08-Sep-2014 07:25 2.0M

Farnell-AD8310-Analo..> 08-Sep-2014 07:24 2.1M

Farnell-AD8361-Rev-D..> 08-Sep-2014 07:23 2.1M

Farnell-2N3906-Fairc..> 08-Sep-2014 07:22 2.1M

Farnell-AD584-Rev-C-..> 08-Sep-2014 07:20 2.2M

Farnell-ADE7753-Rev-..> 08-Sep-2014 07:20 2.3M

Farnell-TLV320AIC23B..> 08-Sep-2014 07:18 2.4M

Farnell-AD586BRZ-Ana..> 08-Sep-2014 07:17 1.6M

Farnell-STM32F405xxS..> 27-Aug-2014 18:27 1.8M

Farnell-MSP430-Hardw..> 29-Jul-2014 10:36 1.1M

Farnell-LM324-Texas-..> 29-Jul-2014 10:32 1.5M

Farnell-LM386-Low-Vo..> 29-Jul-2014 10:32 1.5M

Farnell-NE5532-Texas..> 29-Jul-2014 10:32 1.5M

Farnell-Hex-Inverter..> 29-Jul-2014 10:31 875K

Farnell-AT90USBKey-H..> 29-Jul-2014 10:31 902K

Farnell-AT89C5131-Ha..> 29-Jul-2014 10:31 1.2M

Farnell-MSP-EXP430F5..> 29-Jul-2014 10:31 1.2M

Farnell-Explorer-16-..> 29-Jul-2014 10:31 1.3M

Farnell-TMP006EVM-Us..> 29-Jul-2014 10:30 1.3M

Farnell-Gertboard-Us..> 29-Jul-2014 10:30 1.4M

Farnell-LMP91051-Use..> 29-Jul-2014 10:30 1.4M

Farnell-Thermometre-..> 29-Jul-2014 10:30 1.4M

Farnell-user-manuel-..> 29-Jul-2014 10:29 1.5M

Farnell-fx-3650P-fx-..> 29-Jul-2014 10:29 1.5M

Farnell-2-GBPS-Diffe..> 28-Jul-2014 17:42 2.7M

Farnell-LMT88-2.4V-1..> 28-Jul-2014 17:42 2.8M

Farnell-Octal-Genera..> 28-Jul-2014 17:42 2.8M

Farnell-Dual-MOSFET-..> 28-Jul-2014 17:41 2.8M

Farnell-TLV320AIC325..> 28-Jul-2014 17:41 2.9M

Farnell-SN54LV4053A-..> 28-Jul-2014 17:20 5.9M

Farnell-TAS1020B-USB..> 28-Jul-2014 17:19 6.2M

Farnell-TPS40060-Wid..> 28-Jul-2014 17:19 6.3M

Farnell-TL082-Wide-B..> 28-Jul-2014 17:16 6.3M

Farnell-RF-short-tra..> 28-Jul-2014 17:16 6.3M

Farnell-maxim-integr..> 28-Jul-2014 17:14 6.4M

Farnell-TSV6390-TSV6..> 28-Jul-2014 17:14 6.4M

Farnell-Fast-Charge-..> 28-Jul-2014 17:12 6.4M

Farnell-NVE-datashee..> 28-Jul-2014 17:12 6.5M

Farnell-Excalibur-Hi..> 28-Jul-2014 17:10 2.4M

Farnell-Excalibur-Hi..> 28-Jul-2014 17:10 2.4M

Farnell-REF102-10V-P..> 28-Jul-2014 17:09 2.4M

Farnell-TMS320F28055..> 28-Jul-2014 17:09 2.7M

Farnell-MULTICOMP-Ra..> 22-Jul-2014 12:35 5.9M

Farnell-RASPBERRY-PI..> 22-Jul-2014 12:35 5.9M

Farnell-Dremel-Exper..> 22-Jul-2014 12:34 1.6M

Farnell-STM32F103x8-..> 22-Jul-2014 12:33 1.6M

Farnell-BD6xxx-PDF.htm 22-Jul-2014 12:33 1.6M

Farnell-L78S-STMicro..> 22-Jul-2014 12:32 1.6M

Farnell-RaspiCam-Doc..> 22-Jul-2014 12:32 1.6M

Farnell-SB520-SB5100..> 22-Jul-2014 12:32 1.6M

Farnell-iServer-Micr..> 22-Jul-2014 12:32 1.6M

Farnell-LUMINARY-MIC..> 22-Jul-2014 12:31 3.6M

Farnell-TEXAS-INSTRU..> 22-Jul-2014 12:31 2.4M

Farnell-TEXAS-INSTRU..> 22-Jul-2014 12:30 4.6M

Farnell-CLASS 1-or-2..> 22-Jul-2014 12:30 4.7M

Farnell-TEXAS-INSTRU..> 22-Jul-2014 12:29 4.8M

Farnell-Evaluating-t..> 22-Jul-2014 12:28 4.9M

Farnell-LM3S6952-Mic..> 22-Jul-2014 12:27 5.9M

Farnell-Keyboard-Mou..> 22-Jul-2014 12:27 5.9M

Farnell-Full-Datashe..> 15-Jul-2014 17:08 951K

Farnell-pmbta13_pmbt..> 15-Jul-2014 17:06 959K

Farnell-EE-SPX303N-4..> 15-Jul-2014 17:06 969K

Farnell-Datasheet-NX..> 15-Jul-2014 17:06 1.0M

Farnell-Datasheet-Fa..> 15-Jul-2014 17:05 1.0M

Farnell-MIDAS-un-tra..> 15-Jul-2014 17:05 1.0M

Farnell-SERIAL-TFT-M..> 15-Jul-2014 17:05 1.0M

Farnell-MCOC1-Farnel..> 15-Jul-2014 17:05 1.0M

Farnell-TMR-2-series..> 15-Jul-2014 16:48 787K

Farnell-DC-DC-Conver..> 15-Jul-2014 16:48 781K

Farnell-Full-Datashe..> 15-Jul-2014 16:47 803K

Farnell-TMLM-Series-..> 15-Jul-2014 16:47 810K

Farnell-TEL-5-Series..> 15-Jul-2014 16:47 814K

Farnell-TXL-series-t..> 15-Jul-2014 16:47 829K

Farnell-TEP-150WI-Se..> 15-Jul-2014 16:47 837K

Farnell-AC-DC-Power-..> 15-Jul-2014 16:47 845K

Farnell-TIS-Instruct..> 15-Jul-2014 16:47 845K

Farnell-TOS-tracopow..> 15-Jul-2014 16:47 852K

Farnell-TCL-DC-traco..> 15-Jul-2014 16:46 858K

Farnell-TIS-series-t..> 15-Jul-2014 16:46 875K

Farnell-TMR-2-Series..> 15-Jul-2014 16:46 897K

Farnell-TMR-3-WI-Ser..> 15-Jul-2014 16:46 939K

Farnell-TEN-8-WI-Ser..> 15-Jul-2014 16:46 939K

Farnell-Full-Datashe..> 15-Jul-2014 16:46 947K

Farnell-HIP4081A-Int..> 07-Jul-2014 19:47 1.0M

Farnell-ISL6251-ISL6..> 07-Jul-2014 19:47 1.1M

Farnell-DG411-DG412-..> 07-Jul-2014 19:47 1.0M

Farnell-3367-ARALDIT..> 07-Jul-2014 19:46 1.2M

Farnell-ICM7228-Inte..> 07-Jul-2014 19:46 1.1M

Farnell-Data-Sheet-K..> 07-Jul-2014 19:46 1.2M

Farnell-Silica-Gel-M..> 07-Jul-2014 19:46 1.2M

Farnell-TKC2-Dusters..> 07-Jul-2014 19:46 1.2M

Farnell-CRC-HANDCLEA..> 07-Jul-2014 19:46 1.2M

Farnell-760G-French-..> 07-Jul-2014 19:45 1.2M

Farnell-Decapant-KF-..> 07-Jul-2014 19:45 1.2M

Farnell-1734-ARALDIT..> 07-Jul-2014 19:45 1.2M

Farnell-Araldite-Fus..> 07-Jul-2014 19:45 1.2M

Farnell-fiche-de-don..> 07-Jul-2014 19:44 1.4M

Farnell-safety-data-..> 07-Jul-2014 19:44 1.4M

Farnell-A-4-Hardener..> 07-Jul-2014 19:44 1.4M

Farnell-CC-Debugger-..> 07-Jul-2014 19:44 1.5M

Farnell-MSP430-Hardw..> 07-Jul-2014 19:43 1.8M

Farnell-SmartRF06-Ev..> 07-Jul-2014 19:43 1.6M

Farnell-CC2531-USB-H..> 07-Jul-2014 19:43 1.8M

Farnell-Alimentation..> 07-Jul-2014 19:43 1.8M

Farnell-BK889B-PONT-..> 07-Jul-2014 19:42 1.8M

Farnell-User-Guide-M..> 07-Jul-2014 19:41 2.0M

Farnell-T672-3000-Se..> 07-Jul-2014 19:41 2.0M

Farnell-0050375063-D..> 18-Jul-2014 17:03 2.5M

Farnell-Mini-Fit-Jr-..> 18-Jul-2014 17:03 2.5M

Farnell-43031-0002-M..> 18-Jul-2014 17:03 2.5M

Farnell-0433751001-D..> 18-Jul-2014 17:02 2.5M

Farnell-Cube-3D-Prin..> 18-Jul-2014 17:02 2.5M

Farnell-MTX-Compact-..> 18-Jul-2014 17:01 2.5M

Farnell-MTX-3250-MTX..> 18-Jul-2014 17:01 2.5M

Farnell-ATtiny26-L-A..> 18-Jul-2014 17:00 2.6M

Farnell-MCP3421-Micr..> 18-Jul-2014 17:00 1.2M

Farnell-LM19-Texas-I..> 18-Jul-2014 17:00 1.2M

Farnell-Data-Sheet-S..> 18-Jul-2014 17:00 1.2M

Farnell-LMH6518-Texa..> 18-Jul-2014 16:59 1.3M

Farnell-AD7719-Low-V..> 18-Jul-2014 16:59 1.4M

Farnell-DAC8143-Data..> 18-Jul-2014 16:59 1.5M

Farnell-BGA7124-400-..> 18-Jul-2014 16:59 1.5M

Farnell-SICK-OPTIC-E..> 18-Jul-2014 16:58 1.5M

Farnell-LT3757-Linea..> 18-Jul-2014 16:58 1.6M

Farnell-LT1961-Linea..> 18-Jul-2014 16:58 1.6M

Farnell-PIC18F2420-2..> 18-Jul-2014 16:57 2.5M

Farnell-DS3231-DS-PD..> 18-Jul-2014 16:57 2.5M

Farnell-RDS-80-PDF.htm 18-Jul-2014 16:57 1.3M

Farnell-AD8300-Data-..> 18-Jul-2014 16:56 1.3M

Farnell-LT6233-Linea..> 18-Jul-2014 16:56 1.3M

Farnell-MAX1365-MAX1..> 18-Jul-2014 16:56 1.4M

Farnell-XPSAF5130-PD..> 18-Jul-2014 16:56 1.4M

Farnell-DP83846A-DsP..> 18-Jul-2014 16:55 1.5M

Farnell-Dremel-Exper..> 18-Jul-2014 16:55 1.6M

Farnell-MCOC1-Farnel..> 16-Jul-2014 09:04 1.0M

Farnell-SL3S1203_121..> 16-Jul-2014 09:04 1.1M

Farnell-PN512-Full-N..> 16-Jul-2014 09:03 1.4M

Farnell-SL3S4011_402..> 16-Jul-2014 09:03 1.1M

Farnell-LPC408x-7x 3..> 16-Jul-2014 09:03 1.6M

Farnell-PCF8574-PCF8..> 16-Jul-2014 09:03 1.7M

Farnell-LPC81xM-32-b..> 16-Jul-2014 09:02 2.0M

Farnell-LPC1769-68-6..> 16-Jul-2014 09:02 1.9M

Farnell-Download-dat..> 16-Jul-2014 09:02 2.2M

Farnell-LPC3220-30-4..> 16-Jul-2014 09:02 2.2M

Farnell-LPC11U3x-32-..> 16-Jul-2014 09:01 2.4M

Farnell-SL3ICS1002-1..> 16-Jul-2014 09:01 2.5M

Farnell-T672-3000-Se..> 08-Jul-2014 18:59 2.0M

Farnell-tesaÂ®pack63..> 08-Jul-2014 18:56 2.0M

Farnell-Encodeur-USB..> 08-Jul-2014 18:56 2.0M

Farnell-CC2530ZDK-Us..> 08-Jul-2014 18:55 2.1M

Farnell-2020-Manuel-..> 08-Jul-2014 18:55 2.1M

Farnell-Synchronous-..> 08-Jul-2014 18:54 2.1M

Farnell-Arithmetic-L..> 08-Jul-2014 18:54 2.1M

Farnell-NA555-NE555-..> 08-Jul-2014 18:53 2.2M

Farnell-4-Bit-Magnit..> 08-Jul-2014 18:53 2.2M

Farnell-LM555-Timer-..> 08-Jul-2014 18:53 2.2M

Farnell-L293d-Texas-..> 08-Jul-2014 18:53 2.2M

Farnell-SN54HC244-SN..> 08-Jul-2014 18:52 2.3M

Farnell-MAX232-MAX23..> 08-Jul-2014 18:52 2.3M

Farnell-High-precisi..> 08-Jul-2014 18:51 2.3M

Farnell-SMU-Instrume..> 08-Jul-2014 18:51 2.3M

Farnell-900-Series-B..> 08-Jul-2014 18:50 2.3M

Farnell-BA-Series-Oh..> 08-Jul-2014 18:50 2.3M

Farnell-UTS-Series-S..> 08-Jul-2014 18:49 2.5M

Farnell-270-Series-O..> 08-Jul-2014 18:49 2.3M

Farnell-UTS-Series-S..> 08-Jul-2014 18:49 2.8M

Farnell-Tiva-C-Serie..> 08-Jul-2014 18:49 2.6M

Farnell-UTO-Souriau-..> 08-Jul-2014 18:48 2.8M

Farnell-Clipper-Seri..> 08-Jul-2014 18:48 2.8M

Farnell-SOURIAU-Cont..> 08-Jul-2014 18:47 3.0M

Farnell-851-Series-P..> 08-Jul-2014 18:47 3.0M

Farnell-SL59830-Inte..> 06-Jul-2014 10:07 1.0M

Farnell-ALF1210-PDF.htm 06-Jul-2014 10:06 4.0M

Farnell-AD7171-16-Bi..> 06-Jul-2014 10:06 1.0M

Farnell-Low-Noise-24..> 06-Jul-2014 10:05 1.0M

Farnell-ESCON-Featur..> 06-Jul-2014 10:05 938K

Farnell-74LCX573-Fai..> 06-Jul-2014 10:05 1.9M

Farnell-1N4148WS-Fai..> 06-Jul-2014 10:04 1.9M

Farnell-FAN6756-Fair..> 06-Jul-2014 10:04 850K

Farnell-Datasheet-Fa..> 06-Jul-2014 10:04 861K

Farnell-ES1F-ES1J-fi..> 06-Jul-2014 10:04 867K

Farnell-QRE1113-Fair..> 06-Jul-2014 10:03 879K

Farnell-2N7002DW-Fai..> 06-Jul-2014 10:03 886K

Farnell-FDC2512-Fair..> 06-Jul-2014 10:03 886K

Farnell-FDV301N-Digi..> 06-Jul-2014 10:03 886K

Farnell-S1A-Fairchil..> 06-Jul-2014 10:03 896K

Farnell-BAV99-Fairch..> 06-Jul-2014 10:03 896K

Farnell-74AC00-74ACT..> 06-Jul-2014 10:03 911K

Farnell-NaPiOn-Panas..> 06-Jul-2014 10:02 911K

Farnell-LQ-RELAYS-AL..> 06-Jul-2014 10:02 924K

Farnell-ev-relays-ae..> 06-Jul-2014 10:02 926K

Farnell-ESCON-Featur..> 06-Jul-2014 10:02 931K

Farnell-Amplifier-In..> 06-Jul-2014 10:02 940K

Farnell-Serial-File-..> 06-Jul-2014 10:02 941K

Farnell-Both-the-Del..> 06-Jul-2014 10:01 948K

Farnell-Videk-PDF.htm 06-Jul-2014 10:01 948K

Farnell-EPCOS-173438..> 04-Jul-2014 10:43 3.3M

Farnell-Sensorless-C..> 04-Jul-2014 10:42 3.3M

Farnell-197.31-KB-Te..> 04-Jul-2014 10:42 3.3M

Farnell-PIC12F609-61..> 04-Jul-2014 10:41 3.7M

Farnell-PADO-semi-au..> 04-Jul-2014 10:41 3.7M

Farnell-03-iec-runds..> 04-Jul-2014 10:40 3.7M

Farnell-ACC-Silicone..> 04-Jul-2014 10:40 3.7M

Farnell-Series-TDS10..> 04-Jul-2014 10:39 4.0M

Farnell-03-iec-runds..> 04-Jul-2014 10:40 3.7M

Farnell-0430300011-D..> 14-Jun-2014 18:13 2.0M

Farnell-06-6544-8-PD..> 26-Mar-2014 17:56 2.7M

Farnell-3M-Polyimide..> 21-Mar-2014 08:09 3.9M

Farnell-3M-VolitionT..> 25-Mar-2014 08:18 3.3M

Farnell-10BQ060-PDF.htm 14-Jun-2014 09:50 2.4M

Farnell-10TPB47M-End..> 14-Jun-2014 18:16 3.4M

Farnell-12mm-Size-In..> 14-Jun-2014 09:50 2.4M

Farnell-24AA024-24LC..> 23-Jun-2014 10:26 3.1M

Farnell-50A-High-Pow..> 20-Mar-2014 17:31 2.9M

Farnell-197.31-KB-Te..> 04-Jul-2014 10:42 3.3M

Farnell-1907-2006-PD..> 26-Mar-2014 17:56 2.7M

Farnell-5910-PDF.htm 25-Mar-2014 08:15 3.0M

Farnell-6517b-Electr..> 29-Mar-2014 11:12 3.3M

Farnell-A-True-Syste..> 29-Mar-2014 11:13 3.3M

Farnell-ACC-Silicone..> 04-Jul-2014 10:40 3.7M

Farnell-AD524-PDF.htm 20-Mar-2014 17:33 2.8M

Farnell-ADL6507-PDF.htm 14-Jun-2014 18:19 3.4M

Farnell-ADSP-21362-A..> 20-Mar-2014 17:34 2.8M

Farnell-ALF1210-PDF.htm 04-Jul-2014 10:39 4.0M

Farnell-ALF1225-12-V..> 01-Apr-2014 07:40 3.4M

Farnell-ALF2412-24-V..> 01-Apr-2014 07:39 3.4M

Farnell-AN10361-Phil..> 23-Jun-2014 10:29 2.1M

Farnell-ARADUR-HY-13..> 26-Mar-2014 17:55 2.8M

Farnell-ARALDITE-201..> 21-Mar-2014 08:12 3.7M

Farnell-ARALDITE-CW-..> 26-Mar-2014 17:56 2.7M

Farnell-ATMEL-8-bit-..> 19-Mar-2014 18:04 2.1M

Farnell-ATMEL-8-bit-..> 11-Mar-2014 07:55 2.1M

Farnell-ATmega640-VA..> 14-Jun-2014 09:49 2.5M

Farnell-ATtiny20-PDF..> 25-Mar-2014 08:19 3.6M

Farnell-ATtiny26-L-A..> 13-Jun-2014 18:40 1.8M

Farnell-Alimentation..> 14-Jun-2014 18:24 2.5M

Farnell-Alimentation..> 01-Apr-2014 07:42 3.4M

Farnell-Amplificateu..> 29-Mar-2014 11:11 3.3M

Farnell-An-Improved-..> 14-Jun-2014 09:49 2.5M

Farnell-Atmel-ATmega..> 19-Mar-2014 18:03 2.2M

Farnell-Avvertenze-e..> 14-Jun-2014 18:20 3.3M

Farnell-BC846DS-NXP-..> 13-Jun-2014 18:42 1.6M

Farnell-BC847DS-NXP-..> 23-Jun-2014 10:24 3.3M

Farnell-BF545A-BF545..> 23-Jun-2014 10:28 2.1M

Farnell-BK2650A-BK26..> 29-Mar-2014 11:10 3.3M

Farnell-BT151-650R-N..> 13-Jun-2014 18:40 1.7M

Farnell-BTA204-800C-..> 13-Jun-2014 18:42 1.6M

Farnell-BUJD203AX-NX..> 13-Jun-2014 18:41 1.7M

Farnell-BYV29F-600-N..> 13-Jun-2014 18:42 1.6M

Farnell-BYV79E-serie..> 10-Mar-2014 16:19 1.6M

Farnell-BZX384-serie..> 23-Jun-2014 10:29 2.1M

Farnell-Battery-GBA-..> 14-Jun-2014 18:13 2.0M

Farnell-C.A-6150-C.A..> 14-Jun-2014 18:24 2.5M

Farnell-C.A 8332B-C...> 01-Apr-2014 07:40 3.4M

Farnell-CC2560-Bluet..> 29-Mar-2014 11:14 2.8M

Farnell-CD4536B-Type..> 14-Jun-2014 18:13 2.0M

Farnell-CIRRUS-LOGIC..> 10-Mar-2014 17:20 2.1M

Farnell-CS5532-34-BS..> 01-Apr-2014 07:39 3.5M

Farnell-Cannon-ZD-PD..> 11-Mar-2014 08:13 2.8M

Farnell-Ceramic-tran..> 14-Jun-2014 18:19 3.4M

Farnell-Circuit-Note..> 26-Mar-2014 18:00 2.8M

Farnell-Circuit-Note..> 26-Mar-2014 18:00 2.8M

Farnell-Cles-electro..> 21-Mar-2014 08:13 3.9M

Farnell-Conception-d..> 11-Mar-2014 07:49 2.4M

Farnell-Connectors-N..> 14-Jun-2014 18:12 2.1M

Farnell-Construction..> 14-Jun-2014 18:25 2.5M

Farnell-Controle-de-..> 11-Mar-2014 08:16 2.8M

Farnell-Cordless-dri..> 14-Jun-2014 18:13 2.0M

Farnell-Current-Tran..> 26-Mar-2014 17:58 2.7M

Farnell-Current-Tran..> 26-Mar-2014 17:58 2.7M

Farnell-Current-Tran..> 26-Mar-2014 17:59 2.7M

Farnell-Current-Tran..> 26-Mar-2014 17:59 2.7M

Farnell-DC-Fan-type-..> 14-Jun-2014 09:48 2.5M

Farnell-DC-Fan-type-..> 14-Jun-2014 09:51 1.8M

Farnell-Davum-TMC-PD..> 14-Jun-2014 18:27 2.4M

Farnell-De-la-puissa..> 29-Mar-2014 11:10 3.3M

Farnell-Directive-re..> 25-Mar-2014 08:16 3.0M

Farnell-Documentatio..> 14-Jun-2014 18:26 2.5M

Farnell-Download-dat..> 13-Jun-2014 18:40 1.8M

Farnell-ECO-Series-T..> 20-Mar-2014 08:14 2.5M

Farnell-ELMA-PDF.htm 29-Mar-2014 11:13 3.3M

Farnell-EMC1182-PDF.htm 25-Mar-2014 08:17 3.0M

Farnell-EPCOS-173438..> 04-Jul-2014 10:43 3.3M

Farnell-EPCOS-Sample..> 11-Mar-2014 07:53 2.2M

Farnell-ES2333-PDF.htm 11-Mar-2014 08:14 2.8M

Farnell-Ed.081002-DA..> 19-Mar-2014 18:02 2.5M

Farnell-F28069-Picco..> 14-Jun-2014 18:14 2.0M

Farnell-F42202-PDF.htm 19-Mar-2014 18:00 2.5M

Farnell-FDS-ITW-Spra..> 14-Jun-2014 18:22 3.3M

Farnell-FICHE-DE-DON..> 10-Mar-2014 16:17 1.6M

Farnell-Fastrack-Sup..> 23-Jun-2014 10:25 3.3M

Farnell-Ferric-Chlor..> 29-Mar-2014 11:14 2.8M

Farnell-Fiche-de-don..> 14-Jun-2014 09:47 2.5M

Farnell-Fiche-de-don..> 14-Jun-2014 18:26 2.5M

Farnell-Fluke-1730-E..> 14-Jun-2014 18:23 2.5M

Farnell-GALVA-A-FROI..> 26-Mar-2014 17:56 2.7M

Farnell-GALVA-MAT-Re..> 26-Mar-2014 17:57 2.7M

Farnell-GN-RELAYS-AG..> 20-Mar-2014 08:11 2.6M

Farnell-HC49-4H-Crys..> 14-Jun-2014 18:20 3.3M

Farnell-HFE1600-Data..> 14-Jun-2014 18:22 3.3M

Farnell-HI-70300-Sol..> 14-Jun-2014 18:27 2.4M

Farnell-HUNTSMAN-Adv..> 10-Mar-2014 16:17 1.7M

Farnell-Haute-vitess..> 11-Mar-2014 08:17 2.4M

Farnell-IP4252CZ16-8..> 13-Jun-2014 18:41 1.7M

Farnell-Instructions..> 19-Mar-2014 18:01 2.5M

Farnell-KSZ8851SNL-S..> 23-Jun-2014 10:28 2.1M

Farnell-L-efficacite..> 11-Mar-2014 07:52 2.3M

Farnell-LCW-CQ7P.CC-..> 25-Mar-2014 08:19 3.2M

Farnell-LME49725-Pow..> 14-Jun-2014 09:49 2.5M

Farnell-LOCTITE-542-..> 25-Mar-2014 08:15 3.0M

Farnell-LOCTITE-3463..> 25-Mar-2014 08:19 3.0M

Farnell-LUXEON-Guide..> 11-Mar-2014 07:52 2.3M

Farnell-Leaded-Trans..> 23-Jun-2014 10:26 3.2M

Farnell-Les-derniers..> 11-Mar-2014 07:50 2.3M

Farnell-Loctite3455-..> 25-Mar-2014 08:16 3.0M

Farnell-Low-cost-Enc..> 13-Jun-2014 18:42 1.7M

Farnell-Lubrifiant-a..> 26-Mar-2014 18:00 2.7M

Farnell-MC3510-PDF.htm 25-Mar-2014 08:17 3.0M

Farnell-MC21605-PDF.htm 11-Mar-2014 08:14 2.8M

Farnell-MCF532x-7x-E..> 29-Mar-2014 11:14 2.8M

Farnell-MICREL-KSZ88..> 11-Mar-2014 07:54 2.2M

Farnell-MICROCHIP-PI..> 19-Mar-2014 18:02 2.5M

Farnell-MOLEX-39-00-..> 10-Mar-2014 17:19 1.9M

Farnell-MOLEX-43020-..> 10-Mar-2014 17:21 1.9M

Farnell-MOLEX-43160-..> 10-Mar-2014 17:21 1.9M

Farnell-MOLEX-87439-..> 10-Mar-2014 17:21 1.9M

Farnell-MPXV7002-Rev..> 20-Mar-2014 17:33 2.8M

Farnell-MX670-MX675-..> 14-Jun-2014 09:46 2.5M

Farnell-Microchip-MC..> 13-Jun-2014 18:27 1.8M

Farnell-Microship-PI..> 11-Mar-2014 07:53 2.2M

Farnell-Midas-Active..> 14-Jun-2014 18:17 3.4M

Farnell-Midas-MCCOG4..> 14-Jun-2014 18:11 2.1M

Farnell-Miniature-Ci..> 26-Mar-2014 17:55 2.8M

Farnell-Mistral-PDF.htm 14-Jun-2014 18:12 2.1M

Farnell-Molex-83421-..> 14-Jun-2014 18:17 3.4M

Farnell-Molex-COMMER..> 14-Jun-2014 18:16 3.4M

Farnell-Molex-Crimp-..> 10-Mar-2014 16:27 1.7M

Farnell-Multi-Functi..> 20-Mar-2014 17:38 3.0M

Farnell-NTE_SEMICOND..> 11-Mar-2014 07:52 2.3M

Farnell-NXP-74VHC126..> 10-Mar-2014 16:17 1.6M

Farnell-NXP-BT136-60..> 11-Mar-2014 07:52 2.3M

Farnell-NXP-PBSS9110..> 10-Mar-2014 17:21 1.9M

Farnell-NXP-PCA9555 ..> 11-Mar-2014 07:54 2.2M

Farnell-NXP-PMBFJ620..> 10-Mar-2014 16:16 1.7M

Farnell-NXP-PSMN1R7-..> 10-Mar-2014 16:17 1.6M

Farnell-NXP-PSMN7R0-..> 10-Mar-2014 17:19 2.1M

Farnell-NXP-TEA1703T..> 11-Mar-2014 08:15 2.8M

Farnell-Nilï¬-sk-E-..> 14-Jun-2014 09:47 2.5M

Farnell-Novembre-201..> 20-Mar-2014 17:38 3.3M

Farnell-OMRON-Master..> 10-Mar-2014 16:26 1.8M

Farnell-OSLON-SSL-Ce..> 19-Mar-2014 18:03 2.1M

Farnell-OXPCIE958-FB..> 13-Jun-2014 18:40 1.8M

Farnell-PADO-semi-au..> 04-Jul-2014 10:41 3.7M

Farnell-PBSS5160T-60..> 19-Mar-2014 18:03 2.1M

Farnell-PDTA143X-ser..> 20-Mar-2014 08:12 2.6M

Farnell-PDTB123TT-NX..> 13-Jun-2014 18:43 1.5M

Farnell-PESD5V0F1BL-..> 13-Jun-2014 18:43 1.5M

Farnell-PESD9X5.0L-P..> 13-Jun-2014 18:43 1.6M

Farnell-PIC12F609-61..> 04-Jul-2014 10:41 3.7M

Farnell-PIC18F2455-2..> 23-Jun-2014 10:27 3.1M

Farnell-PIC24FJ256GB..> 14-Jun-2014 09:51 2.4M

Farnell-PMBT3906-PNP..> 13-Jun-2014 18:44 1.5M

Farnell-PMBT4403-PNP..> 23-Jun-2014 10:27 3.1M

Farnell-PMEG4002EL-N..> 14-Jun-2014 18:18 3.4M

Farnell-PMEG4010CEH-..> 13-Jun-2014 18:43 1.6M

Farnell-Panasonic-15..> 23-Jun-2014 10:29 2.1M

Farnell-Panasonic-EC..> 20-Mar-2014 17:36 2.6M

Farnell-Panasonic-EZ..> 20-Mar-2014 08:10 2.6M

Farnell-Panasonic-Id..> 20-Mar-2014 17:35 2.6M

Farnell-Panasonic-Ne..> 20-Mar-2014 17:36 2.6M

Farnell-Panasonic-Ra..> 20-Mar-2014 17:37 2.6M

Farnell-Panasonic-TS..> 20-Mar-2014 08:12 2.6M

Farnell-Panasonic-Y3..> 20-Mar-2014 08:11 2.6M

Farnell-Pico-Spox-Wi..> 10-Mar-2014 16:16 1.7M

Farnell-Pompes-Charg..> 24-Apr-2014 20:23 3.3M

Farnell-Ponts-RLC-po..> 14-Jun-2014 18:23 3.3M

Farnell-Portable-Ana..> 29-Mar-2014 11:16 2.8M

Farnell-Premier-Farn..> 21-Mar-2014 08:11 3.8M

Farnell-Produit-3430..> 14-Jun-2014 09:48 2.5M

Farnell-Proskit-SS-3..> 10-Mar-2014 16:26 1.8M

Farnell-Puissance-ut..> 11-Mar-2014 07:49 2.4M

Farnell-Q48-PDF.htm 23-Jun-2014 10:29 2.1M

Farnell-Radial-Lead-..> 20-Mar-2014 08:12 2.6M

Farnell-Realiser-un-..> 11-Mar-2014 07:51 2.3M

Farnell-Reglement-RE..> 21-Mar-2014 08:08 3.9M

Farnell-Repartiteurs..> 14-Jun-2014 18:26 2.5M

Farnell-S-TRI-SWT860..> 21-Mar-2014 08:11 3.8M

Farnell-SB175-Connec..> 11-Mar-2014 08:14 2.8M

Farnell-SMBJ-Transil..> 29-Mar-2014 11:12 3.3M

Farnell-SOT-23-Multi..> 11-Mar-2014 07:51 2.3M

Farnell-SPLC780A1-16..> 14-Jun-2014 18:25 2.5M

Farnell-SSC7102-Micr..> 23-Jun-2014 10:25 3.2M

Farnell-SVPE-series-..> 14-Jun-2014 18:15 2.0M

Farnell-Sensorless-C..> 04-Jul-2014 10:42 3.3M

Farnell-Septembre-20..> 20-Mar-2014 17:46 3.7M

Farnell-Serie-PicoSc..> 19-Mar-2014 18:01 2.5M

Farnell-Serie-Standa..> 14-Jun-2014 18:23 3.3M

Farnell-Series-2600B..> 20-Mar-2014 17:30 3.0M

Farnell-Series-TDS10..> 04-Jul-2014 10:39 4.0M

Farnell-Signal-PCB-R..> 14-Jun-2014 18:11 2.1M

Farnell-Strangkuhlko..> 21-Mar-2014 08:09 3.9M

Farnell-Supercapacit..> 26-Mar-2014 17:57 2.7M

Farnell-TDK-Lambda-H..> 14-Jun-2014 18:21 3.3M

Farnell-TEKTRONIX-DP..> 10-Mar-2014 17:20 2.0M

Farnell-Tektronix-AC..> 13-Jun-2014 18:44 1.5M

Farnell-Telemetres-l..> 20-Mar-2014 17:46 3.7M

Farnell-Termometros-..> 14-Jun-2014 18:14 2.0M

Farnell-The-essentia..> 10-Mar-2014 16:27 1.7M

Farnell-U2270B-PDF.htm 14-Jun-2014 18:15 3.4M

Farnell-USB-Buccanee..> 14-Jun-2014 09:48 2.5M

Farnell-USB1T11A-PDF..> 19-Mar-2014 18:03 2.1M

Farnell-V4N-PDF.htm 14-Jun-2014 18:11 2.1M

Farnell-WetTantalum-..> 11-Mar-2014 08:14 2.8M

Farnell-XPS-AC-Octop..> 14-Jun-2014 18:11 2.1M

Farnell-XPS-MC16-XPS..> 11-Mar-2014 08:15 2.8M

Farnell-YAGEO-DATA-S..> 11-Mar-2014 08:13 2.8M

Farnell-ZigBee-ou-le..> 11-Mar-2014 07:50 2.4M

Farnell-celpac-SUL84..> 21-Mar-2014 08:11 3.8M

Farnell-china_rohs_o..> 21-Mar-2014 10:04 3.9M

Farnell-cree-Xlamp-X..> 20-Mar-2014 17:34 2.8M

Farnell-cree-Xlamp-X..> 20-Mar-2014 17:35 2.7M

Farnell-cree-Xlamp-X..> 20-Mar-2014 17:31 2.9M

Farnell-cree-Xlamp-m..> 20-Mar-2014 17:32 2.9M

Farnell-cree-Xlamp-m..> 20-Mar-2014 17:32 2.9M

Farnell-ir1150s_fr.p..> 29-Mar-2014 11:11 3.3M

Farnell-manual-bus-p..> 10-Mar-2014 16:29 1.9M

Farnell-propose-plus..> 11-Mar-2014 08:19 2.8M

Farnell-techfirst_se..> 21-Mar-2014 08:08 3.9M

Farnell-testo-205-20..> 20-Mar-2014 17:37 3.0M

Farnell-testo-470-Fo..> 20-Mar-2014 17:38 3.0M

Farnell-uC-OS-III-Br..> 10-Mar-2014 17:20 2.0M

Sefram-7866HD.pdf-PD..> 29-Mar-2014 11:46 472K

Sefram-CAT_ENREGISTR..> 29-Mar-2014 11:46 461K

Sefram-CAT_MESUREURS..> 29-Mar-2014 11:46 435K

Sefram-GUIDE_SIMPLIF..> 29-Mar-2014 11:46 481K

Sefram-GUIDE_SIMPLIF..> 29-Mar-2014 11:46 442K

Sefram-GUIDE_SIMPLIF..> 29-Mar-2014 11:46 422K

Sefram-SP270.pdf-PDF..> 29-Mar-2014 11:46 464KCHAPTER 12 The Fast Fourier Transform There are several ways to calculate the Discrete Fourier Transform (DFT), such as solving simultaneous linear equations or the correlation method described in Chapter 8. The Fast Fourier Transform (FFT) is another method for calculating the DFT. While it produces the same result as the other approaches, it is incredibly more efficient, often reducing the computation time by hundreds. This is the same improvement as flying in a jet aircraft versus walking! If the FFT were not available, many of the techniques described in this book would not be practical. While the FFT only requires a few dozen lines of code, it is one of the most complicated algorithms in DSP. But don't despair! You can easily use published FFT routines without fully understanding the internal workings. Real DFT Using the Complex DFT J.W. Cooley and J.W. Tukey are given credit for bringing the FFT to the world in their paper: "An algorithm for the machine calculation of complex Fourier Series," Mathematics Computation, Vol. 19, 1965, pp 297-301. In retrospect, others had discovered the technique many years before. For instance, the great German mathematician Karl Friedrich Gauss (1777-1855) had used the method more than a century earlier. This early work was largely forgotten because it lacked the tool to make it practical: the digital computer. Cooley and Tukey are honored because they discovered the FFT at the right time, the beginning of the computer revolution. The FFT is based on the complex DFT, a more sophisticated version of the real DFT discussed in the last four chapters. These transforms are named for the way each represents data, that is, using complex numbers or using real numbers. The term complex does not mean that this representation is difficult or complicated, but that a specific type of mathematics is used. Complex mathematics often is difficult and complicated, but that isn't where the name comes from. Chapter 29 discusses the complex DFT and provides the background needed to understand the details of the FFT algorithm. The 226 The Scientist and Engineer's Guide to Digital Signal Processing FIGURE 12-1 Comparing the real and complex DFTs. The real DFT takes an N point time domain signal and creates two N/2% 1 point frequency domain signals. The complex DFT takes two N point time domain signals and creates two N point frequency domain signals. The crosshatched regions shows the values common to the two transforms. Real DFT Complex DFT Time Domain Time Domain Frequency Domain Frequency Domain 0 N-1 0 N-1 0 N-1 0 N/2 0 N/2 0 0 N-1 N-1 N/2 N/2 Real Part Imaginary Part Real Part Imaginary Part Real Part Imaginary Part Time Domain Signal topic of this chapter is simpler: how to use the FFT to calculate the real DFT, without drowning in a mire of advanced mathematics. Since the FFT is an algorithm for calculating the complex DFT, it is important to understand how to transfer real DFT data into and out of the complex DFT format. Figure 12-1 compares how the real DFT and the complex DFT store data. The real DFT transforms an N point time domain signal into two N/2 % 1 point frequency domain signals. The time domain signal is called just that: the time domain signal. The two signals in the frequency domain are called the real part and the imaginary part, holding the amplitudes of the cosine waves and sine waves, respectively. This should be very familiar from past chapters. In comparison, the complex DFT transforms two N point time domain signals into two N point frequency domain signals. The two time domain signals are called the real part and the imaginary part, just as are the frequency domain signals. In spite of their names, all of the values in these arrays are just ordinary numbers. (If you are familiar with complex numbers: the j's are not included in the array values; they are a part of the mathematics. Recall that the operator, Im( ), returns a real number). Chapter 12- The Fast Fourier Transform 227 6000 'NEGATIVE FREQUENCY GENERATION 6010 'This subroutine creates the complex frequency domain from the real frequency domain. 6020 'Upon entry to this subroutine, N% contains the number of points in the signals, and 6030 'REX[ ] and IMX[ ] contain the real frequency domain in samples 0 to N%/2. 6040 'On return, REX[ ] and IMX[ ] contain the complex frequency domain in samples 0 to N%-1. 6050 ' 6060 FOR K% = (N%/2+1) TO (N%-1) 6070 REX[K%] = REX[N%-K%] 6080 IMX[K%] = -IMX[N%-K%] 6090 NEXT K% 6100 ' 6110 RETURN TABLE 12-1 Suppose you have an N point signal, and need to calculate the real DFT by means of the Complex DFT (such as by using the FFT algorithm). First, move the N point signal into the real part of the complex DFT's time domain, and then set all of the samples in the imaginary part to zero. Calculation of the complex DFT results in a real and an imaginary signal in the frequency domain, each composed of N points. Samples 0 through N/2 of these signals correspond to the real DFT's spectrum. As discussed in Chapter 10, the DFT's frequency domain is periodic when the negative frequencies are included (see Fig. 10-9). The choice of a single period is arbitrary; it can be chosen between -1.0 and 0, -0.5 and 0.5, 0 and 1.0, or any other one unit interval referenced to the sampling rate. The complex DFT's frequency spectrum includes the negative frequencies in the 0 to 1.0 arrangement. In other words, one full period stretches from sample 0 to sample N&1 , corresponding with 0 to 1.0 times the sampling rate. The positive frequencies sit between sample 0 and N/2 , corresponding with 0 to 0.5. The other samples, between N/2% 1 and N&1 , contain the negative frequency values (which are usually ignored). Calculating a real Inverse DFT using a complex Inverse DFT is slightly harder. This is because you need to insure that the negative frequencies are loaded in the proper format. Remember, points 0 through N/2 in the complex DFT are the same as in the real DFT, for both the real and the imaginary parts. For the real part, point N/2% 1 is the same as point N/2& 1 , point N/2% 2 is the same as point N/2& 2 , etc. This continues to point N&1 being the same as point 1. The same basic pattern is used for the imaginary part, except the sign is changed. That is, point N/2% 1 is the negative of point N/2& 1 , point N/2% 2 is the negative of point N/2& 2 , etc. Notice that samples 0 and N/2 do not have a matching point in this duplication scheme. Use Fig. 10-9 as a guide to understanding this symmetry. In practice, you load the real DFT's frequency spectrum into samples 0 to N/2 of the complex DFT's arrays, and then use a subroutine to generate the negative frequencies between samples N/2% 1 and N&1 . Table 12-1 shows such a program. To check that the proper symmetry is present, after taking the inverse FFT, look at the imaginary part of the time domain. It will contain all zeros if everything is correct (except for a few parts-permillion of noise, using single precision calculations). 228 The Scientist and Engineer's Guide to Digital Signal Processing FIGURE 12-2 The FFT decomposition. An N point signal is decomposed into N signals each containing a single point. Each stage uses an interlace decomposition, separating the even and odd numbered samples. 1 signal of 16 points 2 signals of 8 points 4 signals of 4 points 8 signals of 2 points 16 signals of 1 point 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0 2 4 6 8 10 12 14 1 3 5 7 9 11 13 15 0 4 8 12 2 6 10 14 1 5 9 13 3 7 11 15 0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15 0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15 How the FFT works The FFT is a complicated algorithm, and its details are usually left to those that specialize in such things. This section describes the general operation of the FFT, but skirts a key issue: the use of complex numbers. If you have a background in complex mathematics, you can read between the lines to understand the true nature of the algorithm. Don't worry if the details elude you; few scientists and engineers that use the FFT could write the program from scratch. In complex notation, the time and frequency domains each contain one signal made up of N complex points. Each of these complex points is composed of two numbers, the real part and the imaginary part. For example, when we talk about complex sample X[42] , it refers to the combination of ReX[42] and ImX[42]. In other words, each complex variable holds two numbers. When two complex variables are multiplied, the four individual components must be combined to form the two components of the product (such as in Eq. 9-1). The following discussion on "How the FFT works" uses this jargon of complex notation. That is, the singular terms: signal, point, sample, and value, refer to the combination of the real part and the imaginary part. The FFT operates by decomposing an N point time domain signal into N time domain signals each composed of a single point. The second step is to calculate the N frequency spectra corresponding to these N time domain signals. Lastly, the N spectra are synthesized into a single frequency spectrum. Figure 12-2 shows an example of the time domain decomposition used in the FFT. In this example, a 16 point signal is decomposed through four Chapter 12- The Fast Fourier Transform 229 Sample numbers Sample numbers in normal order after bit reversal Decimal Binary Decimal Binary 0 0000 0 0000 1 0001 8 1000 2 0010 4 0100 3 0011 12 1100 4 0100 2 0010 5 0101 10 1010 6 0110 6 0100 7 0111 14 1110 8 1000 1 0001 9 1001 9 1001 10 1010 5 0101 11 1011 13 1101 12 1100 3 0011 13 1101 11 1011 14 1110 7 0111 15 1111 15 1111 FIGURE 12-3 The FFT bit reversal sorting. The FFT time domain decomposition can be implemented by sorting the samples according to bit reversed order. separate stages. The first stage breaks the 16 point signal into two signals each consisting of 8 points. The second stage decomposes the data into four signals of 4 points. This pattern continues until there are N signals composed of a single point. An interlaced decomposition is used each time a signal is broken in two, that is, the signal is separated into its even and odd numbered samples. The best way to understand this is by inspecting Fig. 12-2 until you grasp the pattern. There are Log stages required in this decomposition, i.e., 2N a 16 point signal (24) requires 4 stages, a 512 point signal (27) requires 7 stages, a 4096 point signal (212) requires 12 stages, etc. Remember this value, Log ; it will be referenced many times in this chapter. 2N Now that you understand the structure of the decomposition, it can be greatly simplified. The decomposition is nothing more than a reordering of the samples in the signal. Figure 12-3 shows the rearrangement pattern required. On the left, the sample numbers of the original signal are listed along with their binary equivalents. On the right, the rearranged sample numbers are listed, also along with their binary equivalents. The important idea is that the binary numbers are the reversals of each other. For example, sample 3 (0011) is exchanged with sample number 12 (1100). Likewise, sample number 14 (1110) is swapped with sample number 7 (0111), and so forth. The FFT time domain decomposition is usually carried out by a bit reversal sorting algorithm. This involves rearranging the order of the N time domain samples by counting in binary with the bits flipped left-for-right (such as in the far right column in Fig. 12-3). 230 The Scientist and Engineer's Guide to Digital Signal Processing a b c d a 0 b 0 c 0 d 0 A B C D A B C D A B C D e f g h 0 e 0 f 0 g 0 h E F G H F G H E F G H × sinusoid Time Domain Frequency Domain E FIGURE 12-4 The FFT synthesis. When a time domain signal is diluted with zeros, the frequency domain is duplicated. If the time domain signal is also shifted by one sample during the dilution, the spectrum will additionally be multiplied by a sinusoid. The next step in the FFT algorithm is to find the frequency spectra of the 1 point time domain signals. Nothing could be easier; the frequency spectrum of a 1 point signal is equal to itself. This means that nothing is required to do this step. Although there is no work involved, don't forget that each of the 1 point signals is now a frequency spectrum, and not a time domain signal. The last step in the FFT is to combine the N frequency spectra in the exact reverse order that the time domain decomposition took place. This is where the algorithm gets messy. Unfortunately, the bit reversal shortcut is not applicable, and we must go back one stage at a time. In the first stage, 16 frequency spectra (1 point each) are synthesized into 8 frequency spectra (2 points each). In the second stage, the 8 frequency spectra (2 points each) are synthesized into 4 frequency spectra (4 points each), and so on. The last stage results in the output of the FFT, a 16 point frequency spectrum. Figure 12-4 shows how two frequency spectra, each composed of 4 points, are combined into a single frequency spectrum of 8 points. This synthesis must undo the interlaced decomposition done in the time domain. In other words, the frequency domain operation must correspond to the time domain procedure of combining two 4 point signals by interlacing. Consider two time domain signals, abcd and efgh. An 8 point time domain signal can be formed by two steps: dilute each 4 point signal with zeros to make it an Chapter 12- The Fast Fourier Transform 231 + + + + + + + + Eight Point Frequency Spectrum Odd- Four Point Frequency Spectrum Even- Four Point Frequency Spectrum xS xS xS xS FIGURE 12-5 FFT synthesis flow diagram. This shows the method of combining two 4 point frequency spectra into a single 8 point frequency spectrum. The ×S operation means that the signal is multiplied by a sinusoid with an appropriately selected frequency. 2 point input 2 point output xS FIGURE 12-6 The FFT butterfly. This is the basic calculation element in the FFT, taking two complex points and converting them into two other complex points. 8 point signal, and then add the signals together. That is, abcd becomes a0b0c0d0, and efgh becomes 0e0f0g0h. Adding these two 8 point signals produces aebfcgdh. As shown in Fig. 12-4, diluting the time domain with zeros corresponds to a duplication of the frequency spectrum. Therefore, the frequency spectra are combined in the FFT by duplicating them, and then adding the duplicated spectra together. In order to match up when added, the two time domain signals are diluted with zeros in a slightly different way. In one signal, the odd points are zero, while in the other signal, the even points are zero. In other words, one of the time domain signals (0e0f0g0h in Fig. 12-4) is shifted to the right by one sample. This time domain shift corresponds to multiplying the spectrum by a sinusoid. To see this, recall that a shift in the time domain is equivalent to convolving the signal with a shifted delta function. This multiplies the signal's spectrum with the spectrum of the shifted delta function. The spectrum of a shifted delta function is a sinusoid (see Fig 11-2). Figure 12-5 shows a flow diagram for combining two 4 point spectra into a single 8 point spectrum. To reduce the situation even more, notice that Fig. 12- 5 is formed from the basic pattern in Fig 12-6 repeated over and over. 232 The Scientist and Engineer's Guide to Digital Signal Processing Time Domain Data Frequency Domain Data Bit Reversal Data Sorting Overhead Overhead Calculation Decomposition Synthesis Time Domain Frequency Domain Butterfly FIGURE 12-7 Flow diagram of the FFT. This is based on three steps: (1) decompose an N point time domain signal into N signals each containing a single point, (2) find the spectrum of each of the N point signals (nothing required), and (3) synthesize the N frequency spectra into a single frequency spectrum. Loop for each Butterfly Loop for Leach sub-DFT Loop for Log2N stages This simple flow diagram is called a butterfly due to its winged appearance. The butterfly is the basic computational element of the FFT, transforming two complex points into two other complex points. Figure 12-7 shows the structure of the entire FFT. The time domain decomposition is accomplished with a bit reversal sorting algorithm. Transforming the decomposed data into the frequency domain involves nothing and therefore does not appear in the figure. The frequency domain synthesis requires three loops. The outer loop runs through the Log stages (i.e., each level in Fig. 12-2, starting from the bottom 2N and moving to the top). The middle loop moves through each of the individual frequency spectra in the stage being worked on (i.e., each of the boxes on any one level in Fig. 12-2). The innermost loop uses the butterfly to calculate the points in each frequency spectra (i.e., looping through the samples inside any one box in Fig. 12-2). The overhead boxes in Fig. 12-7 determine the beginning and ending indexes for the loops, as well as calculating the sinusoids needed in the butterflies. Now we come to the heart of this chapter, the actual FFT programs. Chapter 12- The Fast Fourier Transform 233 5000 'COMPLEX DFT BY CORRELATION 5010 'Upon entry, N% contains the number of points in the DFT, and 5020 'XR[ ] and XI[ ] contain the real and imaginary parts of the time domain. 5030 'Upon return, REX[ ] and IMX[ ] contain the frequency domain data. 5040 'All signals run from 0 to N%-1. 5050 ' 5060 PI = 3.14159265 'Set constants 5070 ' 5080 FOR K% = 0 TO N%-1 'Zero REX[ ] and IMX[ ], so they can be used 5090 REX[K%] = 0 'as accumulators during the correlation 5100 IMX[K%] = 0 5110 NEXT K% 5120 ' 5130 FOR K% = 0 TO N%-1 'Loop for each value in frequency domain 5140 FOR I% = 0 TO N%-1 'Correlate with the complex sinusoid, SR & SI 5150 ' 5160 SR = COS(2*PI*K%*I%/N%) 'Calculate complex sinusoid 5170 SI = -SIN(2*PI*K%*I%/N%) 5180 REX[K%] = REX[K%] + XR[I%]*SR - XI[I%]*SI 5190 IMX[K%] = IMX[K%] + XR[I%]*SI + XI[I%]*SR 5200 ' 5210 NEXT I% 5220 NEXT K% 5230 ' 5240 RETURN TABLE 12-2 FFT Programs As discussed in Chapter 8, the real DFT can be calculated by correlating the time domain signal with sine and cosine waves (see Table 8-2). Table 12-2 shows a program to calculate the complex DFT by the same method. In an apples-to-apples comparison, this is the program that the FFT improves upon. Tables 12-3 and 12-4 show two different FFT programs, one in FORTRAN and one in BASIC. First we will look at the BASIC routine in Table 12-4. This subroutine produces exactly the same output as the correlation technique in Table 12-2, except it does it much faster. The block diagram in Fig. 12-7 can be used to identify the different sections of this program. Data are passed to this FFT subroutine in the arrays: REX[ ] and IMX[ ], each running from sample 0 to N&1 . Upon return from the subroutine, REX[ ] and IMX[ ] are overwritten with the frequency domain data. This is another way that the FFT is highly optimized; the same arrays are used for the input, intermediate storage, and output. This efficient use of memory is important for designing fast hardware to calculate the FFT. The term in-place computation is used to describe this memory usage. While all FFT programs produce the same numerical result, there are subtle variations in programming that you need to look out for. Several of these 234 The Scientist and Engineer's Guide to Digital Signal Processing TABLE 12-3 The Fast Fourier Transform in FORTRAN. Data are passed to this subroutine in the variables X( ) and M. The integer, M, is the base two logarithm of the length of the DFT, i.e., M = 8 for a 256 point DFT, M = 12 for a 4096 point DFT, etc. The complex array, X( ), holds the time domain data upon entering the DFT. Upon return from this subroutine, X( ) is overwritten with the frequency domain data. Take note: this subroutine requires that the input and output signals run from X(1) through X(N), rather than the customary X(0) through X(N-1). SUBROUTINE FFT(X,M) COMPLEX X(4096),U,S,T PI=3.14159265 N=2**M DO 20 L=1,M LE=2**(M+1-L) LE2=LE/2 U=(1.0,0.0) S=CMPLX(COS(PI/FLOAT(LE2)),-SIN(PI/FLOAT(LE2))) DO 20 J=1,LE2 DO 10 I=J,N,LE IP=I+LE2 T=X(I)+X(IP) X(IP)=(X(I)-X(IP))*U 10 X(I)=T 20 U=U*S ND2=N/2 NM1=N-1 J=1 DO 50 I=1,NM1 IF(I.GE.J) GO TO 30 T=X(J) X(J)=X(I) X(I)=T 30 K=ND2 40 IF(K.GE.J) GO TO 50 J=J-K K=K/2 GO TO 40 50 J=J+K RETURN END of these differences are illustrated by the FORTRAN program listed in Table 12-3. This program uses an algorithm called decimation in frequency, while the previously described algorithm is called decimation in time. In a decimation in frequency algorithm, the bit reversal sorting is done after the three nested loops. There are also FFT routines that completely eliminate the bit reversal sorting. None of these variations significantly improve the performance of the FFT, and you shouldn't worry about which one you are using. The important differences between FFT algorithms concern how data are passed to and from the subroutines. In the BASIC program, data enter and leave the subroutine in the arrays REX[ ] and IMX[ ], with the samples running from index 0 to N&1 . In the FORTRAN program, data are passed in the complex array X( ), with the samples running from 1 to N. Since this is an array of complex variables, each sample in X( ) consists of two numbers, a real part and an imaginary part. The length of the DFT must also be passed to these subroutines. In the BASIC program, the variable N% is used for this purpose. In comparison, the FORTRAN program uses the variable M, which is defined to equal Log . For instance, M will be 2N Chapter 12- The Fast Fourier Transform 235 TABLE 12-4 The Fast Fourier Transform in BASIC. 1000 'THE FAST FOURIER TRANSFORM 1010 'Upon entry, N% contains the number of points in the DFT, REX[ ] and 1020 'IMX[ ] contain the real and imaginary parts of the input. Upon return, 1030 'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1. 1040 ' 1050 PI = 3.14159265 'Set constants 1060 NM1% = N%-1 1070 ND2% = N%/2 1080 M% = CINT(LOG(N%)/LOG(2)) 1090 J% = ND2% 1100 ' 1110 FOR I% = 1 TO N%-2 'Bit reversal sorting 1120 IF I% >= J% THEN GOTO 1190 1130 TR = REX[J%] 1140 TI = IMX[J%] 1150 REX[J%] = REX[I%] 1160 IMX[J%] = IMX[I%] 1170 REX[I%] = TR 1180 IMX[I%] = TI 1190 K% = ND2% 1200 IF K% > J% THEN GOTO 1240 1210 J% = J%-K% 1220 K% = K%/2 1230 GOTO 1200 1240 J% = J%+K% 1250 NEXT I% 1260 ' 1270 FOR L% = 1 TO M% 'Loop for each stage 1280 LE% = CINT(2^L%) 1290 LE2% = LE%/2 1300 UR = 1 1310 UI = 0 1320 SR = COS(PI/LE2%) 'Calculate sine & cosine values 1330 SI = -SIN(PI/LE2%) 1340 FOR J% = 1 TO LE2% 'Loop for each sub DFT 1350 JM1% = J%-1 1360 FOR I% = JM1% TO NM1% STEP LE% 'Loop for each butterfly 1370 IP% = I%+LE2% 1380 TR = REX[IP%]*UR - IMX[IP%]*UI 'Butterfly calculation 1390 TI = REX[IP%]*UI + IMX[IP%]*UR 1400 REX[IP%] = REX[I%]-TR 1410 IMX[IP%] = IMX[I%]-TI 1420 REX[I%] = REX[I%]+TR 1430 IMX[I%] = IMX[I%]+TI 1440 NEXT I% 1450 TR = UR 1460 UR = TR*SR - UI*SI 1470 UI = TR*SI + UI*SR 1480 NEXT J% 1490 NEXT L% 1500 ' 1510 RETURN 236 The Scientist and Engineer's Guide to Digital Signal Processing 2000 'INVERSE FAST FOURIER TRANSFORM SUBROUTINE 2010 'Upon entry, N% contains the number of points in the IDFT, REX[ ] and 2020 'IMX[ ] contain the real and imaginary parts of the complex frequency domain. 2030 'Upon return, REX[ ] and IMX[ ] contain the complex time domain. 2040 'All signals run from 0 to N%-1. 2050 ' 2060 FOR K% = 0 TO N%-1 'Change the sign of IMX[ ] 2070 IMX[K%] = -IMX[K%] 2080 NEXT K% 2090 ' 2100 GOSUB 1000 'Calculate forward FFT (Table 12-3) 2110 ' 2120 FOR I% = 0 TO N%-1 'Divide the time domain by N% and 2130 REX[I%] = REX[I%]/N% 'change the sign of IMX[ ] 2140 IMX[I%] = -IMX[I%]/N% 2150 NEXT I% 2160 ' 2170 RETURN TABLE 12-5 8 for a 256 point DFT, 12 for a 4096 point DFT, etc. The point is, the programmer who writes an FFT subroutine has many options for interfacing with the host program. Arrays that run from 1 to N, such as in the FORTRAN program, are especially aggravating. Most of the DSP literature (including this book) explains algorithms assuming the arrays run from sample 0 to N&1 . For instance, if the arrays run from 1 to N, the symmetry in the frequency domain is around points 1 and N/2% 1 , rather than points 0 and N/2 , Using the complex DFT to calculate the real DFT has another interesting advantage. The complex DFT is more symmetrical between the time and frequency domains than the real DFT. That is, the duality is stronger. Among other things, this means that the Inverse DFT is nearly identical to the Forward DFT. In fact, the easiest way to calculate an Inverse FFT is to calculate a Forward FFT, and then adjust the data. Table 12-5 shows a subroutine for calculating the Inverse FFT in this manner. Suppose you copy one of these FFT algorithms into your computer program and start it running. How do you know if it is operating properly? Two tricks are commonly used for debugging. First, start with some arbitrary time domain signal, such as from a random number generator, and run it through the FFT. Next, run the resultant frequency spectrum through the Inverse FFT and compare the result with the original signal. They should be identical, except round-off noise (a few parts-per-million for single precision). The second test of proper operation is that the signals have the correct symmetry. When the imaginary part of the time domain signal is composed of all zeros (the normal case), the frequency domain of the complex DFT will be symmetrical around samples 0 and N/2 , as previously described. Chapter 12- The Fast Fourier Transform 237 EQUATION 12-1 DFT execution time. The time required to calculate a DFT by correlation is proportional to the length of the DFT squared. ExecutionTime ’ kDFT N2 EQUATION 12-2 FFT execution time. The time required to calculate a DFT using the FFT is proportional to N multiplied by the logarithm of N. ExecutionTime ’ kFFT N log2N Likewise, when this correct symmetry is present in the frequency domain, the Inverse DFT will produce a time domain that has an imaginary part composes of all zeros (plus round-off noise). These debugging techniques are essential for using the FFT; become familiar with them. Speed and Precision Comparisons When the DFT is calculated by correlation (as in Table 12-2), the program uses two nested loops, each running through N points. This means that the total number of operations is proportional to N times N. The time to complete the program is thus given by: where N is the number of points in the DFT and kDFT is a constant of proportionality. If the sine and cosine values are calculated within the nested loops, kDFT is equal to about 25 microseconds on a Pentium at 100 MHz. If you precalculate the sine and cosine values and store them in a look-up-table, kDFT drops to about 7 microseconds. For example, a 1024 point DFT will require about 25 seconds, or nearly 25 milliseconds per point. That's slow! Using this same strategy we can derive the execution time for the FFT. The time required for the bit reversal is negligible. In each of the Log stages 2N there are N/2 butterfly computations. This means the execution time for the program is approximated by: The value of kFFT is about 10 microseconds on a 100 MHz Pentium system. A 1024 point FFT requires about 70 milliseconds to execute, or 70 microseconds per point. This is more than 300 times faster than the DFT calculated by correlation! Not only is NLog less than , it increases much more slowly as N 2N N 2 becomes larger. For example, a 32 point FFT is about ten times faster than the correlation method. However, a 4096 point FFT is one-thousand times faster. For small values of N (say, 32 to 128), the FFT is important. For large values of N (1024 and above), the FFT is absolutely critical. Figure 12-8 compares the execution times of the two algorithms in a graphical form. 238 The Scientist and Engineer's Guide to Digital Signal Processing Number points in DFT 8 16 32 64 128 256 512 1024 2048 4096 0.001 0.01 0.1 1 10 100 1000 FFT correlation correlation w/LUT FIGURE 12-8 Execution times for calculating the DFT. The correlation method refers to the algorithm described in Table 12-2. This method can be made faster by precalculating the sine and cosine values and storing them in a look-up table (LUT). The FFT (Table 12-3) is the fastest algorithm when the DFT is greater than 16 points long. The times shown are for a Pentium processor at 100 MHz. Execution time (seconds) Number of points in DFT 16 32 64 128 256 512 1024 0 10 20 30 40 50 60 70 FFT correlation FIGURE 12-9 DFT precision. Since the FFT calculates the DFT faster than the correlation method, it also calculates it with less round-off error. Error (parts per million) The FFT has another advantage besides raw speed. The FFT is calculated more precisely because the fewer number of calculations results in less round-off error. This can be demonstrated by taking the FFT of an arbitrary signal, and then running the frequency spectrum through an Inverse FFT. This reconstructs the original time domain signal, except for the addition of roundoff noise from the calculations. A single number characterizing this noise can be obtained by calculating the standard deviation of the difference between the two signals. For comparison, this same procedure can be repeated using a DFT calculated by correlation, and a corresponding Inverse DFT. How does the round-off noise of the FFT compare to the DFT by correlation? See for yourself in Fig. 12-9. Further Speed Increases There are several techniques for making the FFT even faster; however, the improvements are only about 20-40%. In one of these methods, the time Chapter 12- The Fast Fourier Transform 239 4000 'INVERSE FFT FOR REAL SIGNALS 4010 'Upon entry, N% contains the number of points in the IDFT, REX[ ] and 4020 'IMX[ ] contain the real and imaginary parts of the frequency domain running from 4030 'index 0 to N%/2. The remaining samples in REX[ ] and IMX[ ] are ignored. 4040 'Upon return, REX[ ] contains the real time domain, IMX[ ] contains zeros. 4050 ' 4060 ' 4070 FOR K% = (N%/2+1) TO (N%-1) 'Make frequency domain symmetrical 4080 REX[K%] = REX[N%-K%] '(as in Table 12-1) 4090 IMX[K%] = -IMX[N%-K%] 4100 NEXT K% 4110 ' 4120 FOR K% = 0 TO N%-1 'Add real and imaginary parts together 4130 REX[K%] = REX[K%]+IMX[K%] 4140 NEXT K% 4150 ' 4160 GOSUB 3000 'Calculate forward real DFT (TABLE 12-6) 4170 ' 4180 FOR I% = 0 TO N%-1 'Add real and imaginary parts together 4190 REX[I%] = (REX[I%]+IMX[I%])/N% 'and divide the time domain by N% 4200 IMX[I%] = 0 4210 NEXT I% 4220 ' 4230 RETURN TABLE 12-6 domain decomposition is stopped two stages early, when each signals is composed of only four points. Instead of calculating the last two stages, highly optimized code is used to jump directly into the frequency domain, using the simplicity of four point sine and cosine waves. Another popular algorithm eliminates the wasted calculations associated with the imaginary part of the time domain being zero, and the frequency spectrum being symmetrical. In other words, the FFT is modified to calculate the real DFT, instead of the complex DFT. These algorithms are called the real FFT and the real Inverse FFT (or similar names). Expect them to be about 30% faster than the conventional FFT routines. Tables 12-6 and 12-7 show programs for these algorithms. There are two small disadvantages in using the real FFT. First, the code is about twice as long. While your computer doesn't care, you must take the time to convert someone else's program to run on your computer. Second, debugging these programs is slightly harder because you cannot use symmetry as a check for proper operation. These algorithms force the imaginary part of the time domain to be zero, and the frequency domain to have left-right symmetry. For debugging, check that these programs produce the same output as the conventional FFT algorithms. Figures 12-10 and 12-11 illustrate how the real FFT works. In Fig. 12-10, (a) and (b) show a time domain signal that consists of a pulse in the real part, and all zeros in the imaginary part. Figures (c) and (d) show the corresponding frequency spectrum. As previously described, the frequency domain's real part has an even symmetry around sample 0 and sample N/2 , while the imaginary part has an odd symmetry around these same points. 240 The Scientist and Engineer's Guide to Digital Signal Processing Sample number 0 16 32 48 64 -1 0 1 2 63 a. Real part Freqeuncy 0 16 32 48 -8 -4 0 4 8 c. Real part (even symmetry) 63 Frequency 0 16 32 48 -8 -4 0 4 8 d. Imaginary part (odd symmetry) 63 Time Domain Frequency Domain Sample number 0 16 32 48 64 -1 0 1 2 63 b. Imaginary part FIGURE 12-10 Real part symmetry of the DFT. Amplitude Amplitude Amplitude Amplitude Now consider Fig. 12-11, where the pulse is in the imaginary part of the time domain, and the real part is all zeros. The symmetry in the frequency domain is reversed; the real part is odd, while the imaginary part is even. This situation will be discussed in Chapter 29. For now, take it for granted that this is how the complex DFT behaves. What if there is a signal in both parts of the time domain? By additivity, the frequency domain will be the sum of the two frequency spectra. Now the key element: a frequency spectrum composed of these two types of symmetry can be perfectly separated into the two component signals. This is achieved by the even/odd decomposition discussed in Chapter 6. In other words, two real DFT's can be calculated for the price of single FFT. One of the signals is placed in the real part of the time domain, and the other signal is placed in the imaginary part. After calculating the complex DFT (via the FFT, of course), the spectra are separated using the even/odd decomposition. When two or more signals need to be passed through the FFT, this technique reduces the execution time by about 40%. The improvement isn't a full factor of two because of the calculation time required for the even/odd decomposition. This is a relatively simple technique with few pitfalls, nothing like writing an FFT routine from scratch. Chapter 12- The Fast Fourier Transform 241 Sample number 0 16 32 48 64 -1 0 1 2 63 a. Real part Frequency 0 16 32 48 -8 -4 0 4 8 c. Real part (odd symmetry) 63 Frequency 0 16 32 48 -8 -4 0 4 8 d. Imaginary part (even symmetry) 63 Time Domain Frequency Domain Sample number 0 16 32 48 64 -1 0 1 2 63 b. Imaginary part FIGURE 12-11 Imaginary part symmetry of the DFT. Amplitude Amplitude Amplitude Amplitude The next step is to modify the algorithm to calculate a single DFT faster. It's ugly, but here is how it is done. The input signal is broken in half by using an interlaced decomposition. The N/2 even points are placed into the real part of the time domain signal, while the N/2 odd points go into the imaginary part. An N/2 point FFT is then calculated, requiring about one-half the time as an N point FFT. The resulting frequency domain is then separated by the even/odd decomposition, resulting in the frequency spectra of the two interlaced time domain signals. These two frequency spectra are then combined into a single spectrum, just as in the last synthesis stage of the FFT. To close this chapter, consider that the FFT is to Digital Signal Processing what the transistor is to electronics. It is a foundation of the technology; everyone in the field knows its characteristics and how to use it. However, only a small number of specialists really understand the details of the internal workings. 242 The Scientist and Engineer's Guide to Digital Signal Processing 3000 'FFT FOR REAL SIGNALS 3010 'Upon entry, N% contains the number of points in the DFT, REX[ ] contains 3020 'the real input signal, while values in IMX[ ] are ignored. Upon return, 3030 'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1. 3040 ' 3050 NH% = N%/2-1 'Separate even and odd points 3060 FOR I% = 0 TO NH% 3070 REX(I%) = REX(2*I%) 3080 IMX(I%) = REX(2*I%+1) 3090 NEXT I% 3100 ' 3110 N% = N%/2 'Calculate N%/2 point FFT 3120 GOSUB 1000 '(GOSUB 1000 is the FFT in Table 12-3) 3130 N% = N%*2 3140 ' 3150 NM1% = N%-1 'Even/odd frequency domain decomposition 3160 ND2% = N%/2 3170 N4% = N%/4-1 3180 FOR I% = 1 TO N4% 3190 IM% = ND2%-I% 3200 IP2% = I%+ND2% 3210 IPM% = IM%+ND2% 3220 REX(IP2%) = (IMX(I%) + IMX(IM%))/2 3230 REX(IPM%) = REX(IP2%) 3240 IMX(IP2%) = -(REX(I%) - REX(IM%))/2 3250 IMX(IPM%) = -IMX(IP2%) 3260 REX(I%) = (REX(I%) + REX(IM%))/2 3270 REX(IM%) = REX(I%) 3280 IMX(I%) = (IMX(I%) - IMX(IM%))/2 3290 IMX(IM%) = -IMX(I%) 3300 NEXT I% 3310 REX(N%*3/4) = IMX(N%/4) 3320 REX(ND2%) = IMX(0) 3330 IMX(N%*3/4) = 0 3340 IMX(ND2%) = 0 3350 IMX(N%/4) = 0 3360 IMX(0) = 0 3370 ' 3380 PI = 3.14159265 'Complete the last FFT stage 3390 L% = CINT(LOG(N%)/LOG(2)) 3400 LE% = CINT(2^L%) 3410 LE2% = LE%/2 3420 UR = 1 3430 UI = 0 3440 SR = COS(PI/LE2%) 3450 SI = -SIN(PI/LE2%) 3460 FOR J% = 1 TO LE2% 3470 JM1% = J%-1 3480 FOR I% = JM1% TO NM1% STEP LE% 3490 IP% = I%+LE2% 3500 TR = REX[IP%]*UR - IMX[IP%]*UI 3510 TI = REX[IP%]*UI + IMX[IP%]*UR 3520 REX[IP%] = REX[I%]-TR 3530 IMX[IP%] = IMX[I%]-TI 3540 REX[I%] = REX[I%]+TR 3550 IMX[I%] = IMX[I%]+TI 3560 NEXT I% 3570 TR = UR 3580 UR = TR*SR - UI*SI 3590 UI = TR*SI + UI*SR 3600 NEXT J% 3610 RETURN TABLE 12-7 CHAPTER 15 EQUATION 15-1 Equation of the moving average filter. In this equation, x[ ] is the input signal, y[ ] is the output signal, and M is the number of points used in the moving average. This equation only uses points on one side of the output sample being calculated. y[i ] ’ 1 M j M&1 j’ 0 x [ i %j ] y [80] ’ x [80] % x [81] % x [82] % x [83] % x [84] 5 Moving Average Filters The moving average is the most common filter in DSP, mainly because it is the easiest digital filter to understand and use. In spite of its simplicity, the moving average filter is optimal for a common task: reducing random noise while retaining a sharp step response. This makes it the premier filter for time domain encoded signals. However, the moving average is the worst filter for frequency domain encoded signals, with little ability to separate one band of frequencies from another. Relatives of the moving average filter include the Gaussian, Blackman, and multiplepass moving average. These have slightly better performance in the frequency domain, at the expense of increased computation time. Implementation by Convolution As the name implies, the moving average filter operates by averaging a number of points from the input signal to produce each point in the output signal. In equation form, this is written: Where x [ ] is the input signal, y [ ] is the output signal, and M is the number of points in the average. For example, in a 5 point moving average filter, point 80 in the output signal is given by: 278 The Scientist and Engineer's Guide to Digital Signal Processing y [80] ’ x [78] % x [79] % x [80] % x [81] % x [82] 5 100 'MOVING AVERAGE FILTER 110 'This program filters 5000 samples with a 101 point moving 120 'average filter, resulting in 4900 samples of filtered data. 130 ' 140 DIM X[4999] 'X[ ] holds the input signal 150 DIM Y[4999] 'Y[ ] holds the output signal 160 ' 170 GOSUB XXXX 'Mythical subroutine to load X[ ] 180 ' 190 FOR I% = 50 TO 4949 'Loop for each point in the output signal 200 Y[I%] = 0 'Zero, so it can be used as an accumulator 210 FOR J% = -50 TO 50 'Calculate the summation 220 Y[I%] = Y[I%] + X(I%+J%] 230 NEXT J% 240 Y[I%] = Y[I%]/101 'Complete the average by dividing 250 NEXT I% 260 ' 270 END TABLE 15-1 As an alternative, the group of points from the input signal can be chosen symmetrically around the output point: This corresponds to changing the summation in Eq. 15-1 from: j ’ 0 to M&1 , to: j ’ &(M&1) /2 to (M&1) /2 . For instance, in an 11 point moving average filter, the index, j, can run from 0 to 11 (one side averaging) or -5 to 5 (symmetrical averaging). Symmetrical averaging requires that M be an odd number. Programming is slightly easier with the points on only one side; however, this produces a relative shift between the input and output signals. You should recognize that the moving average filter is a convolution using a very simple filter kernel. For example, a 5 point filter has the filter kernel: þ 0, 0, 1/5, 1/5, 1/5, 1/5, 1/5, 0, 0 þ . That is, the moving average filter is a convolution of the input signal with a rectangular pulse having an area of one. Table 15-1 shows a program to implement the moving average filter. Noise Reduction vs. Step Response Many scientists and engineers feel guilty about using the moving average filter. Because it is so very simple, the moving average filter is often the first thing tried when faced with a problem. Even if the problem is completely solved, there is still the feeling that something more should be done. This situation is truly ironic. Not only is the moving average filter very good for many applications, it is optimal for a common problem, reducing random white noise while keeping the sharpest step response. Chapter 15- Moving Average Filters 279 Sample number 0 100 200 300 400 500 -1 0 1 2 a. Original signal Sample number 0 100 200 300 400 500 -1 0 1 2 b. 11 point moving average FIGURE 15-1 Example of a moving average filter. In (a), a rectangular pulse is buried in random noise. In (b) and (c), this signal is filtered with 11 and 51 point moving average filters, respectively. As the number of points in the filter increases, the noise becomes lower; however, the edges becoming less sharp. The moving average filter is the optimal solution for this problem, providing the lowest noise possible for a given edge sharpness. Sample number 0 100 200 300 400 500 -1 0 1 2 c. 51 point moving average Amplitude Amplitude Amplitude Figure 15-1 shows an example of how this works. The signal in (a) is a pulse buried in random noise. In (b) and (c), the smoothing action of the moving average filter decreases the amplitude of the random noise (good), but also reduces the sharpness of the edges (bad). Of all the possible linear filters that could be used, the moving average produces the lowest noise for a given edge sharpness. The amount of noise reduction is equal to the square-root of the number of points in the average. For example, a 100 point moving average filter reduces the noise by a factor of 10. To understand why the moving average if the best solution, imagine we want to design a filter with a fixed edge sharpness. For example, let's assume we fix the edge sharpness by specifying that there are eleven points in the rise of the step response. This requires that the filter kernel have eleven points. The optimization question is: how do we choose the eleven values in the filter kernel to minimize the noise on the output signal? Since the noise we are trying to reduce is random, none of the input points is special; each is just as noisy as its neighbor. Therefore, it is useless to give preferential treatment to any one of the input points by assigning it a larger coefficient in the filter kernel. The lowest noise is obtained when all the input samples are treated equally, i.e., the moving average filter. (Later in this chapter we show that other filters are essentially as good. The point is, no filter is better than the simple moving average). 280 The Scientist and Engineer's Guide to Digital Signal Processing EQUATION 15-2 Frequency response of an M point moving average filter. The frequency, f, runs between 0 and 0.5. For f ’ 0, use: H[ f ] ’ 1 H [ f ] ’ sin(Bf M ) M sin(Bf ) Frequency 0 0.1 0.2 0.3 0.4 0.5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 3 point 11 point 31 point FIGURE 15-2 Frequency response of the moving average filter. The moving average is a very poor low-pass filter, due to its slow roll-off and poor stopband attenuation. These curves are generated by Eq. 15-2. Amplitude Frequency Response Figure 15-2 shows the frequency response of the moving average filter. It is mathematically described by the Fourier transform of the rectangular pulse, as discussed in Chapter 11: The roll-off is very slow and the stopband attenuation is ghastly. Clearly, the moving average filter cannot separate one band of frequencies from another. Remember, good performance in the time domain results in poor performance in the frequency domain, and vice versa. In short, the moving average is an exceptionally good smoothing filter (the action in the time domain), but an exceptionally bad low-pass filter (the action in the frequency domain). Relatives of the Moving Average Filter In a perfect world, filter designers would only have to deal with time domain or frequency domain encoded information, but never a mixture of the two in the same signal. Unfortunately, there are some applications where both domains are simultaneously important. For instance, television signals fall into this nasty category. Video information is encoded in the time domain, that is, the shape of the waveform corresponds to the patterns of brightness in the image. However, during transmission the video signal is treated according to its frequency composition, such as its total bandwidth, how the carrier waves for sound & color are added, elimination & restoration of the DC component, etc. As another example, electromagnetic interference is best understood in the frequency domain, even if Chapter 15- Moving Average Filters 281 Sample number 0 6 12 18 24 0.0 0.1 0.2 2 pass 4 pass 1 pass a. Filter kernel Sample number 0 6 12 18 24 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1 pass 4 pass 2 pass b. Step response Frequency 0 0.1 0.2 0.3 0.4 0.5 -120 -100 -80 -60 -40 -20 0 20 40 1 pass 2 pass 4 pass d. Frequency response (dB) FIGURE 15-3 Characteristics of multiple-pass moving average filters. Figure (a) shows the filter kernels resulting from passing a seven point moving average filter over the data once, twice and four times. Figure (b) shows the corresponding step responses, while (c) and (d) show the corresponding frequency responses. FFT Integrate 20 Log( ) Amplitude Amplitude Frequency 0 0.1 0.2 0.3 0.4 0.5 0.00 0.25 0.50 0.75 1.00 1.25 1 pass 2 pass 4 pass c. Frequency response Amplitude (dB) Amplitude the signal's information is encoded in the time domain. For instance, the temperature monitor in a scientific experiment might be contaminated with 60 hertz from the power lines, 30 kHz from a switching power supply, or 1320 kHz from a local AM radio station. Relatives of the moving average filter have better frequency domain performance, and can be useful in these mixed domain applications. Multiple-pass moving average filters involve passing the input signal through a moving average filter two or more times. Figure 15-3a shows the overall filter kernel resulting from one, two and four passes. Two passes are equivalent to using a triangular filter kernel (a rectangular filter kernel convolved with itself). After four or more passes, the equivalent filter kernel looks like a Gaussian (recall the Central Limit Theorem). As shown in (b), multiple passes produce an "s" shaped step response, as compared to the straight line of the single pass. The frequency responses in (c) and (d) are given by Eq. 15-2 multiplied by itself for each pass. That is, each time domain convolution results in a multiplication of the frequency spectra. 282 The Scientist and Engineer's Guide to Digital Signal Processing Figure 15-4 shows the frequency response of two other relatives of the moving average filter. When a pure Gaussian is used as a filter kernel, the frequency response is also a Gaussian, as discussed in Chapter 11. The Gaussian is important because it is the impulse response of many natural and manmade systems. For example, a brief pulse of light entering a long fiber optic transmission line will exit as a Gaussian pulse, due to the different paths taken by the photons within the fiber. The Gaussian filter kernel is also used extensively in image processing because it has unique properties that allow fast two-dimensional convolutions (see Chapter 24). The second frequency response in Fig. 15-4 corresponds to using a Blackman window as a filter kernel. (The term window has no meaning here; it is simply part of the accepted name of this curve). The exact shape of the Blackman window is given in Chapter 16 (Eq. 16-2, Fig. 16-2); however, it looks much like a Gaussian. How are these relatives of the moving average filter better than the moving average filter itself? Three ways: First, and most important, these filters have better stopband attenuation than the moving average filter. Second, the filter kernels taper to a smaller amplitude near the ends. Recall that each point in the output signal is a weighted sum of a group of samples from the input. If the filter kernel tapers, samples in the input signal that are farther away are given less weight than those close by. Third, the step responses are smooth curves, rather than the abrupt straight line of the moving average. These last two are usually of limited benefit, although you might find applications where they are genuine advantages. The moving average filter and its relatives are all about the same at reducing random noise while maintaining a sharp step response. The ambiguity lies in how the risetime of the step response is measured. If the risetime is measured from 0% to 100% of the step, the moving average filter is the best you can do, as previously shown. In comparison, measuring the risetime from 10% to 90% makes the Blackman window better than the moving average filter. The point is, this is just theoretical squabbling; consider these filters equal in this parameter. The biggest difference in these filters is execution speed. Using a recursive algorithm (described next), the moving average filter will run like lightning in your computer. In fact, it is the fastest digital filter available. Multiple passes of the moving average will be correspondingly slower, but still very quick. In comparison, the Gaussian and Blackman filters are excruciatingly slow, because they must use convolution. Think a factor of ten times the number of points in the filter kernel (based on multiplication being about 10 times slower than addition). For example, expect a 100 point Gaussian to be 1000 times slower than a moving average using recursion. Recursive Implementation A tremendous advantage of the moving average filter is that it can be implemented with an algorithm that is very fast. To understand this Chapter 15- Moving Average Filters 283 FIGURE 15-4 Frequency response of the Blackman window and Gaussian filter kernels. Both these filters provide better stopband attenuation than the moving average filter. This has no advantage in removing random noise from time domain encoded signals, but it can be useful in mixed domain problems. The disadvantage of these filters is that they must use convolution, a terribly slow algorithm. Frequency 0 0.1 0.2 0.3 0.4 0.5 -140 -120 -100 -80 -60 -40 -20 0 20 Gaussian Blackman Amplitude (dB) y [50] ’ x [47] % x [48] % x [49] % x [50] % x [51] % x [52] % x [53] y [51] ’ x [48] % x [49] % x [50] % x [51] % x [52] % x [53] % x [54] y [51] ’ y [50] % x [54] & x [47] EQUATION 15-3 Recursive implementation of the moving average filter. In this equation, x[ ] is the input signal, y[ ] is the output signal, M is the number of points in the moving average (an odd number). Before this equation can be used, the first point in the signal must be calculated using a standard summation. y [i ] ’ y [i &1] % x [i %p] & x [i &q] q ’ p % 1 where: p ’ (M&1) /2 algorithm, imagine passing an input signal, x [ ], through a seven point moving average filter to form an output signal, y [ ]. Now look at how two adjacent output points, y [50] and y [51], are calculated: These are nearly the same calculation; points x [48] through x [53] must be added for y [50], and again for y [51]. If y [50] has already been calculated, the most efficient way to calculate y [51] is: Once y [51] has been found using y [50], then y [52] can be calculated from sample y [51], and so on. After the first point is calculated in y [ ], all of the other points can be found with only a single addition and subtraction per point. This can be expressed in the equation: Notice that this equation use two sources of data to calculate each point in the output: points from the input and previously calculated points from the output. This is called a recursive equation, meaning that the result of one calculation 284 The Scientist and Engineer's Guide to Digital Signal Processing 100 'MOVING AVERAGE FILTER IMPLEMENTED BY RECURSION 110 'This program filters 5000 samples with a 101 point moving 120 'average filter, resulting in 4900 samples of filtered data. 130 'A double precision accumulator is used to prevent round-off drift. 140 ' 150 DIM X[4999] 'X[ ] holds the input signal 160 DIM Y[4999] 'Y[ ] holds the output signal 170 DEFDBL ACC 'Define the variable ACC to be double precision 180 ' 190 GOSUB XXXX 'Mythical subroutine to load X[ ] 200 ' 210 ACC = 0 'Find Y[50] by averaging points X[0] to X[100] 220 FOR I% = 0 TO 100 230 ACC = ACC + X[I%] 240 NEXT I% 250 Y[[50] = ACC/101 260 ' 'Recursive moving average filter (Eq. 15-3) 270 FOR I% = 51 TO 4949 280 ACC = ACC + X[I%+50] - X[I%-51] 290 Y[I%] = ACC 300 NEXT I% 310 ' 320 END TABLE 15-2 CHAPTER 6 Convolution Convolution is a mathematical way of combining two signals to form a third signal. It is the single most important technique in Digital Signal Processing. Using the strategy of impulse decomposition, systems are described by a signal called the impulse response. Convolution is important because it relates the three signals of interest: the input signal, the output signal, and the impulse response. This chapter presents convolution from two different viewpoints, called the input side algorithm and the output side algorithm. Convolution provides the mathematical framework for DSP; there is nothing more important in this book. The Delta Function and Impulse Response The previous chapter describes how a signal can be decomposed into a group of components called impulses. An impulse is a signal composed of all zeros, except a single nonzero point. In effect, impulse decomposition provides a way to analyze signals one sample at a time. The previous chapter also presented the fundamental concept of DSP: the input signal is decomposed into simple additive components, each of these components is passed through a linear system, and the resulting output components are synthesized (added). The signal resulting from this divide-and-conquer procedure is identical to that obtained by directly passing the original signal through the system. While many different decompositions are possible, two form the backbone of signal processing: impulse decomposition and Fourier decomposition. When impulse decomposition is used, the procedure can be described by a mathematical operation called convolution. In this chapter (and most of the following ones) we will only be dealing with discrete signals. Convolution also applies to continuous signals, but the mathematics is more complicated. We will look at how continious signals are processed in Chapter 13. Figure 6-1 defines two important terms used in DSP. The first is the delta function, symbolized by the Greek letter delta, *[n]. The delta function is a normalized impulse, that is, sample number zero has a value of one, while 108 The Scientist and Engineer's Guide to Digital Signal Processing all other samples have a value of zero. For this reason, the delta function is frequently called the unit impulse. The second term defined in Fig. 6-1 is the impulse response. As the name suggests, the impulse response is the signal that exits a system when a delta function (unit impulse) is the input. If two systems are different in any way, they will have different impulse responses. Just as the input and output signals are often called x[n] and y[n] , the impulse response is usually given the symbol, h[n]. Of course, this can be changed if a more descriptive name is available, for instance, f [n] might be used to identify the impulse response of a filter. Any impulse can be represented as a shifted and scaled delta function. Consider a signal, a[n] , composed of all zeros except sample number 8, which has a value of -3. This is the same as a delta function shifted to the right by 8 samples, and multiplied by -3. In equation form: a[n] ’ &3*[n&8]. Make sure you understand this notation, it is used in nearly all DSP equations. If the input to a system is an impulse, such as &3*[n&8] , what is the system's output? This is where the properties of homogeneity and shift invariance are used. Scaling and shifting the input results in an identical scaling and shifting of the output. If *[n] results in h[n] , it follows that &3*[n&8] results in &3h[n&8] . In words, the output is a version of the impulse response that has been shifted and scaled by the same amount as the delta function on the input. If you know a system's impulse response, you immediately know how it will react to any impulse. Convolution Let's summarize this way of understanding how a system changes an input signal into an output signal. First, the input signal can be decomposed into a set of impulses, each of which can be viewed as a scaled and shifted delta function. Second, the output resulting from each impulse is a scaled and shifted version of the impulse response. Third, the overall output signal can be found by adding these scaled and shifted impulse responses. In other words, if we know a system's impulse response, then we can calculate what the output will be for any possible input signal. This means we know everything about the system. There is nothing more that can be learned about a linear system's characteristics. (However, in later chapters we will show that this information can be represented in different forms). The impulse response goes by a different name in some applications. If the system being considered is a filter, the impulse response is called the filter kernel, the convolution kernel, or simply, the kernel. In image processing, the impulse response is called the point spread function. While these terms are used in slightly different ways, they all mean the same thing, the signal produced by a system when the input is a delta function. Chapter 6- Convolution 109 System -2 -1 0 1 2 3 4 5 6 -1 0 1 2 -2 -1 0 1 2 3 4 5 6 -1 0 1 2 *[n] h[n] Delta Impulse Response Linear Function FIGURE 6-1 Definition of delta function and impulse response. The delta function is a normalized impulse. All of its samples have a value of zero, except for sample number zero, which has a value of one. The Greek letter delta, *[n] , is used to identify the delta function. The impulse response of a linear system, usually denoted by h[n] , is the output of the system when the input is a delta function. x[n] h[n] = y[n] x[n] y[n] Linear System h[n] FIGURE 6-2 How convolution is used in DSP. The output signal from a linear system is equal to the input signal convolved with the system's impulse response. Convolution is denoted by a star when writing equations. Convolution is a formal mathematical operation, just as multiplication, addition, and integration. Addition takes two numbers and produces a third number, while convolution takes two signals and produces a third signal. Convolution is used in the mathematics of many fields, such as probability and statistics. In linear systems, convolution is used to describe the relationship between three signals of interest: the input signal, the impulse response, and the output signal. Figure 6-2 shows the notation when convolution is used with linear systems. An input signal, x[n] , enters a linear system with an impulse response, h[n] , resulting in an output signal, y[n] . In equation form: x[n] t h[n] ’ y[n] . Expressed in words, the input signal convolved with the impulse response is equal to the output signal. Just as addition is represented by the plus, +, and multiplication by the cross, ×, convolution is represented by the star, t. It is unfortunate that most programming languages also use the star to indicate multiplication. A star in a computer program means multiplication, while a star in an equation means convolution. 110 The Scientist and Engineer's Guide to Digital Signal Processing Sample number 0 10 20 30 40 50 60 70 80 90 100 110 -2 -1 0 1 2 3 4 S 0 10 20 30 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 S 0 10 20 30 -0.02 0.00 0.02 0.04 0.06 0.08 a. Low-pass Filter b. High-pass Filter Sample number 0 10 20 30 40 50 60 70 80 -2 -1 0 1 2 3 4 Sample number 0 10 20 30 40 50 60 70 80 90 100 110 -2 -1 0 1 2 3 4 Sample number 0 10 20 30 40 50 60 70 80 -2 -1 0 1 2 3 4 Sample number Sample number Input Signal Impulse Response Output Signal Amplitude Amplitude Amplitude Amplitude Amplitude Amplitude FIGURE 6-3 Examples of low-pass and high-pass filtering using convolution. In this example, the input signal is a few cycles of a sine wave plus a slowly rising ramp. These two components are separated by using properly selected impulse responses. Figure 6-3 shows convolution being used for low-pass and high-pass filtering. The example input signal is the sum of two components: three cycles of a sine wave (representing a high frequency), plus a slowly rising ramp (composed of low frequencies). In (a), the impulse response for the low-pass filter is a smooth arch, resulting in only the slowly changing ramp waveform being passed to the output. Similarly, the high-pass filter, (b), allows only the more rapidly changing sinusoid to pass. Figure 6-4 illustrates two additional examples of how convolution is used to process signals. The inverting attenuator, (a), flips the signal top-for-bottom, and reduces its amplitude. The discrete derivative (also called the first difference), shown in (b), results in an output signal related to the slope of the input signal. Notice the lengths of the signals in Figs. 6-3 and 6-4. The input signals are 81 samples long, while each impulse response is composed of 31 samples. In most DSP applications, the input signal is hundreds, thousands, or even millions of samples in length. The impulse response is usually much shorter, say, a few points to a few hundred points. The mathematics behind convolution doesn't restrict how long these signals are. It does, however, specify the length of the output signal. The length of the output signal is Chapter 6- Convolution 111 S 0 10 20 30 -2.00 -1.00 0.00 1.00 2.00 S 0 10 20 30 -2.00 -1.00 0.00 1.00 2.00 a. Inverting Attenuator b. Discrete Derivative Sample number 0 10 20 30 40 50 60 70 80 90 100 110 -2 -1 0 1 2 3 4 Sample number 0 10 20 30 40 50 60 70 80 90 100 110 -2 -1 0 1 2 3 4 Sample number 0 10 20 30 40 50 60 70 80 -2 -1 0 1 2 3 4 Sample number 0 10 20 30 40 50 60 70 80 -2 -1 0 1 2 3 4 Input Signal Impulse Response Output Signal Sample number Sample number Amplitude Amplitude Amplitude Amplitude Amplitude Amplitude FIGURE 6-4 Examples of signals being processed using convolution. Many signal processing tasks use very simple impulse responses. As shown in these examples, dramatic changes can be achieved with only a few nonzero points. equal to the length of the input signal, plus the length of the impulse response, minus one. For the signals in Figs. 6-3 and 6-4, each output signal is: 81% 31& 1 ’ 111 samples long. The input signal runs from sample 0 to 80, the impulse response from sample 0 to 30, and the output signal from sample 0 to 110. Now we come to the detailed mathematics of convolution. As used in Digital Signal Processing, convolution can be understood in two separate ways. The first looks at convolution from the viewpoint of the input signal. This involves analyzing how each sample in the input signal contributes to many points in the output signal. The second way looks at convolution from the viewpoint of the output signal. This examines how each sample in the output signal has received information from many points in the input signal. Keep in mind that these two perspectives are different ways of thinking about the same mathematical operation. The first viewpoint is important because it provides a conceptual understanding of how convolution pertains to DSP. The second viewpoint describes the mathematics of convolution. This typifies one of the most difficult tasks you will encounter in DSP: making your conceptual understanding fit with the jumble of mathematics used to communicate the ideas. 112 The Scientist and Engineer's Guide to Digital Signal Processing 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 0 1 2 3 4 5 6 7 8 -3 -2 -1 0 1 2 3 0 1 2 3 -3 -2 -1 0 1 2 3 x[n] h[n] y[n] FIGURE 6-5 Example convolution problem. A nine point input signal, convolved with a four point impulse response, results in a twelve point output signal. Each point in the input signal contributes a scaled and shifted impulse response to the output signal. These nine scaled and shifted impulse responses are shown in Fig. 6-6. Now examine sample x[8] , the last point in the input signal. This sample is at index number eight, and has a value of -0.5. As shown in the lower-right graph of Fig. 6-6, x[8] results in an impulse response that has been shifted to the right by eight points and multiplied by -0.5. Place holding zeros have been added at points 0-7. Lastly, examine the effect of points x[0] and x[7] . Both these samples have a value of zero, and therefore produce output components consisting of all zeros. The Input Side Algorithm Figure 6-5 shows a simple convolution problem: a 9 point input signal, x[n] , is passed through a system with a 4 point impulse response, h[n] , resulting in a 9% 4& 1 ’ 12 point output signal, y[n] . In mathematical terms, x[n] is convolved with h[n] to produce y[n] . This first viewpoint of convolution is based on the fundamental concept of DSP: decompose the input, pass the components through the system, and synthesize the output. In this example, each of the nine samples in the input signal will contribute a scaled and shifted version of the impulse response to the output signal. These nine signals are shown in Fig. 6-6. Adding these nine signals produces the output signal, y[n] . Let's look at several of these nine signals in detail. We will start with sample number four in the input signal, i.e., x[4] . This sample is at index number four, and has a value of 1.4. When the signal is decomposed, this turns into an impulse represented as: 1.4*[n&4]. After passing through the system, the resulting output component will be: 1.4 h[n&4]. This signal is shown in the center box of the nine signals in Fig. 6-6. Notice that this is the impulse response, h[n] , multiplied by 1.4, and shifted four samples to the right. Zeros have been added at samples 0-3 and at samples 8-11 to serve as place holders. To make this more clear, Fig. 6-6 uses squares to represent the data points that come from the shifted and scaled impulse response, and diamonds for the added zeros. Chapter 6- Convolution 113 FIGURE 6-6 Output signal components for the convolution in Fig. 6-5. In these signals, each point that results from a scaled and shifted impulse response is represented by a square marker. The remaining data points, represented by diamonds, are zeros that have been added as place holders. 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 In this example, x[n] is a nine point signal and h[n] is a four point signal. In our next example, shown in Fig. 6-7, we will reverse the situation by making x[n] a four point signal, and h[n] a nine point signal. The same two waveforms are used, they are just swapped. As shown by the output signal components, the four samples in x[n] result in four shifted and scaled versions of the nine point impulse response. Just as before, leading and trailing zeros are added as place holders. But wait just one moment! The output signal in Fig. 6-7 is identical to the output signal in Fig. 6-5. This isn't a mistake, but an important property. Convolution is commutative: a[n]tb[n] ’ b[n]ta[n] . The mathematics does not care which is the input signal and which is the impulse response, only that two signals are convolved with each other. Although the mathematics may allow it, exchanging the two signals has no physical meaning in system theory. The input signal and impulse response are two totally different things and exchanging them doesn't make sense. What the commutative property provides is a mathematical tool for manipulating equations to achieve various results. 114 The Scientist and Engineer's Guide to Digital Signal Processing TABLE 6-1 100 'CONVOLUTION USING THE INPUT SIDE ALGORITHM 110 ' 120 DIM X[80] 'The input signal, 81 points 130 DIM H[30] 'The impulse response, 31 points 140 DIM Y[110] 'The output signal, 111 points 150 ' 160 GOSUB XXXX 'Mythical subroutine to load X[ ] and H[ ] 170 ' 180 FOR I% = 0 TO 110 'Zero the output array 190 Y(I%) = 0 200 NEXT I% 210 ' 220 FOR I% = 0 TO 80 'Loop for each point in X[ ] 230 FOR J% = 0 TO 30 'Loop for each point in H[ ] 240 Y[I%+J%] = Y[I%+J%] + X[I%]tH[J%] 250 NEXT J% 260 NEXT I% '(remember, t is multiplication in programs!) 270 ' 280 GOSUB XXXX 'Mythical subroutine to store Y[ ] 290 ' 300 END A program for calculating convolutions using the input side algorithm is shown in Table 6-1. Remember, the programs in this book are meant to convey algorithms in the simplest form, even at the expense of good programming style. For instance, all of the input and output is handled in mythical subroutines (lines 160 and 280), meaning we do not define how these operations are conducted. Do not skip over these programs; they are a key part of the material and you need to understand them in detail. The program convolves an 81 point input signal, held in array X[ ], with a 31 point impulse response, held in array H[ ], resulting in a 111 point output signal, held in array Y[ ]. These are the same lengths shown in Figs. 6-3 and 6-4. Notice that the names of these arrays use upper case letters. This is a violation of the naming conventions previously discussed, because upper case letters are reserved for frequency domain signals. Unfortunately, the simple BASIC used in this book does not allow lower case variable names. Also notice that line 240 uses a star for multiplication. Remember, a star in a program means multiplication, while a star in an equation means convolution. A star in text (such as documentation or program comments) can mean either. The mythical subroutine in line 160 places the input signal into X[ ] and the impulse response into H[ ]. Lines 180-200 set all of the values in Y[ ] to zero. This is necessary because Y[ ] is used as an accumulator to sum the output components as they are calculated. Lines 220 to 260 are the heart of the program. The FOR statement in line 220 controls a loop that steps through each point in the input signal, X[ ]. For each sample in the input signal, an inner loop (lines 230-250) calculates a scaled and shifted version of the impulse response, and adds it to the array accumulating the output signal, Y[ ]. This nested loop structure (one loop within another loop) is a key characteristic of convolution programs; become familiar with it. Chapter 6- Convolution 115 FIGURE 6-7 A second example of convolution. The waveforms for the input signal and impulse response are exchanged from the example of Fig. 6-5. Since convolution is commutative, the output signals for the two examples are identical. 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 0 1 2 3 4 5 6 7 8 -3 -2 -1 0 1 2 3 0 1 2 3 -3 -2 -1 0 1 2 3 x[n] h[n] y[n] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 contribution from x[ ] h[n- ] 0 0 1 1 2 2 3 3 Output signal components Keeping the indexing straight in line 240 can drive you crazy! Let's say we are halfway through the execution of this program, so that we have just begun action on sample X[40], i.e., I% = 40. The inner loop runs through each point in the impulse response doing three things. First, the impulse response is scaled by multiplying it by the value of the input sample. If this were the only action taken by the inner loop, line 240 could be written, Y[J%] = X[40]tH[J%]. Second, the scaled impulse is shifted 40 samples to the right by adding this number to the index used in the output signal. This second action would change line 240 to: Y[40+J%] = X[40]tH[J%]. Third, Y[ ] must accumulate (synthesize) all the signals resulting from each sample in the input signal. Therefore, the new information must be added to the information that is already in the array. This results in the final command: Y[40+J%] = Y[40+J%] + X[40]tH[J%]. Study this carefully; it is very confusing, but very important. 116 The Scientist and Engineer's Guide to Digital Signal Processing The Output Side Algorithm The first viewpoint of convolution analyzes how each sample in the input signal affects many samples in the output signal. In this second viewpoint, we reverse this by looking at individual samples in the output signal, and finding the contributing points from the input. This is important from both mathematical and practical standpoints. Suppose that we are given some input signal and impulse response, and want to find the convolution of the two. The most straightforward method would be to write a program that loops through the output signal, calculating one sample on each loop cycle. Likewise, equations are written in the form: y[n] ’ some combination of other variables. That is, sample n in the output signal is equal to some combination of the many values in the input signal and impulse response. This requires a knowledge of how each sample in the output signal can be calculated independently of all other samples in the output signal. The output side algorithm provides this information. Let's look at an example of how a single point in the output signal is influenced by several points from the input. The example point we will use is y[6] in Fig. 6-5. This point is equal to the sum of all the sixth points in the nine output components, shown in Fig. 6-6. Now, look closely at these nine output components and identify which can affect y[6] . That is, find which of these nine signals contains a nonzero sample at the sixth position. Five of the output components only have added zeros (the diamond markers) at the sixth sample, and can therefore be ignored. Only four of the output components are capable of having a nonzero value in the sixth position. These are the output components generated from the input samples: x[3], x[4], x[5], and x[6] . By adding the sixth sample from each of these output components, y[6] is determined as: y[6] ’ x[3]h[3] % x[4]h[2] % x[5]h[1] % x[6]h[0] . That is, four samples from the input signal are multiplied by the four samples in the impulse response, and the products added. Figure 6-8 illustrates the output side algorithm as a convolution machine, a flow diagram of how convolution occurs. Think of the input signal, x[n] , and the output signal, y[n] , as fixed on the page. The convolution machine, everything inside the dashed box, is free to move left and right as needed. The convolution machine is positioned so that its output is aligned with the output sample being calculated. Four samples from the input signal fall into the inputs of the convolution machine. These values are multiplied by the indicated samples in the impulse response, and the products are added. This produces the value for the output signal, which drops into its proper place. For example, y[6] i s s h own b e i n g c a l c u l a t e d f r om t h e f o u r i n p u t s amp l e s : x[3], x[4], x[5], and x[6] . To calculate y[7] , the convolution machine moves one sample to the right. This results in another four samples entering the machine, x[4] through x[7] , and the value for y[7] dropping into the proper place. This process is repeated for all points in the output signal needing to be calculated. Chapter 6- Convolution 117 0 1 2 3 4 5 6 7 8 -3 -2 -1 0 1 2 3 -3 -2 -1 0 3.0 2.0 1.0 0.0 1.0 2.0 3.0 x[n] y[n] h[n] 0 -1 -2 -3 2 3 1 (flipped) FIGURE 6-8 The convolution machine. This is a flow diagram showing how each sample in the output signal is influenced by the input signal and impulse response. See the text for details. 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 The arrangement of the impulse response inside the convolution machine is very important. The impulse response is flipped left-for-right. This places sample number zero on the right, and increasingly positive sample numbers running to the left. Compare this to the normal impulse response in Fig. 6-5 to understand the geometry of this flip. Why is this flip needed? It simply falls out of the mathematics. The impulse response describes how each point in the input signal affects the output signal. This results in each point in the output signal being affected by points in the input signal weighted by a flipped impulse response. 118 The Scientist and Engineer's Guide to Digital Signal Processing FIGURE 6-9 The convolution machine in action. Figures (a) through (d) show the convolution machine set to calculate four different output signal samples, y[0], y[3], y[8], and y[11]. 0 1 2 3 4 5 6 7 8 -3 -2 -1 0 1 2 3 -3 -2 -1 0 3.0 2.0 1.0 0.0 1.0 2.0 3.0 x[n] y[n] h[n] 0 -1 -2 -3 2 3 1 (flipped) a. Set to calculate y[0] 0 1 2 3 4 5 6 7 8 -3 -2 -1 0 1 2 3 -3 -2 -1 0 3.0 2.0 1.0 0.0 1.0 2.0 3.0 x[n] y[n] h[n] 0 -1 -2 -3 2 3 1 (flipped) b. Set to calculate y[3] 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 Figure 6-9 shows the convolution machine being used to calculate several samples in the output signal. This diagram also illustrates a real nuisance in convolution. In (a), the convolution machine is located fully to the left with its output aimed at y[0] . In this position, it is trying to receive input from samples: x[&3], x[&2], x[&1], and x[0] . The problem is, three of these samples: x[&3], x[&2], and x[&1] , do not exist! This same dilemma arises in (d), where the convolution machine tries to accept samples to the right of the defined input signal, points x[9], x[10], and x[11] . One way to handle this problem is by inventing the nonexistent samples. This involves adding samples to the ends of the input signal, with each of the added samples having a value of zero. This is called padding the signal with zeros. Instead of trying to access a nonexistent value, the convolution machine receives a sample that has a value of zero. Since this zero is eliminated during the multiplication, the result is mathematically the same as ignoring the nonexistent inputs. Chapter 6- Convolution 119 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 0 1 2 3 4 5 6 7 8 9 10 11 -3 -2 -1 0 1 2 3 0 1 2 3 4 5 6 7 8 -3 -2 -1 0 1 2 3 -3 -2 -1 0 3.0 2.0 1.0 0.0 1.0 2.0 3.0 x[n] y[n] h[n] 0 -1 -2 -3 2 3 1 (flipped) c. Set to calculate y[8] 0 1 2 3 4 5 6 7 8 -3 -2 -1 0 1 2 3 -3 -2 -1 0 3.0 2.0 1.0 0.0 1.0 2.0 3.0 x[n] y[n] h[n] 0 -1 -2 -3 2 3 1 (flipped) d. Set to calculate y[11] Figure 6-9 (continued) The important part is that the far left and far right samples in the output signal are based on incomplete information. In DSP jargon, the impulse response is not fully immersed in the input signal. If the impulse response is M points in length, the first and last M&1 samples in the output signal are based on less information than the samples between. This is analogous to an electronic circuit requiring a certain amount of time to stabilize after the power is applied. The difference is that this transient is easy to ignore in electronics, but very prominent in DSP. Figure 6-10 shows an example of the trouble these end effects can cause. The input signal is a sine wave plus a DC component. The desire is to remove the DC part of the signal, while leaving the sine wave intact. This calls for a highpass filter, such as the impulse response shown in the figure. The problem is, the first and last 30 points are a mess! The shape of these end regions can be understood by imagining the input signal padded with 30 zeros on the left side, samples x[&1] through x[&30] , and 30 zeros on the right, samples x[81] through x[110] . The output signal can then be viewed as a filtered version of this longer waveform. These "end effect" problems are widespread in 120 The Scientist and Engineer's Guide to Digital Signal Processing EQUATION 6-1 The convolution summation. This is the formal definition of convolution, written in the shorthand: y [n] ’ x [n] t h[n]. In this equation, h[n] is an M point signal with indexes running from 0 to M-1. y [i ] ’ jM&1 j ’0 h[ j ] x [i&j ] DSP. As a general rule, expect that the beginning and ending samples in processed signals will be quite useless. Now the math. Using the convolution machine as a guideline, we can write the standard equation for convolution. If x[n] is an N point signal running from 0 to N-1, and h[n] is an M point signal running from 0 to M-1, the convolution of the two: y[n] ’ x[n] t h[n], is an N+M-1 point signal running from 0 to N+M-2, given by: This equation is called the convolution sum. It allows each point in the output signal to be calculated independently of all other points in the output signal. The index, i, determines which sample in the output signal is being calculated, and therefore corresponds to the left-right position of the convolution machine. In computer programs performing convolution, a loop makes this index run through each sample in the output signal. To calculate one of the output samples, the index, j, is used inside of the convolution machine. As j runs through 0 to M-1, each sample in the impulse response, h[ j], is multiplied by the proper sample from the input signal, x[i& j ]. All these products are added to produce the output sample being calculated. Study Eq. 6-1 until you fully understand how it is implemented by the convolution machine. Much of DSP is based on this equation. (Don't be confused by the n in y[n] ’ x[n] t h[n]. This is merely a place holder to indicate that some variable is the index into the array. Sometimes the equations are written: y[ ] ’ x[ ] t h[ ], just to avoid having to bring in a meaningless symbol). Table 6-2 shows a program for performing convolutions using the output side algorithm, a direct use of Eq. 6-1. This program produces the same output signal as the program for the input side algorithm, shown previously in Table 6-1. Notice the main difference between these two programs: the input side algorithm loops through each sample in the input signal (line 220 of Table 6- 1), while the output side algorithm loops through each sample in the output signal (line 180 of Table 6-2). Here is a detailed operation of this program. The FOR-NEXT loop in lines 180 to 250 steps through each sample in the output signal, using I% as the index. For each of these values, an inner loop, composed of lines 200 to 230, calculates the value of the output sample, Y[I%]. The value of Y[I%] is set to zero in line 190, allowing it to accumulate the products inside of the convolution machine. The FOR-NEXT loop in lines 200 to 240 provide a direct implementation of Eq. 6-1. The index, J%, steps through each Chapter 6- Convolution 121 sample in the impulse response. Line 230 provides the multiplication of each sample in the impulse response, H[J%], with the appropriate sample from the input signal, X[I%-J%], and adds the result to the accumulator. In line 230, the sample taken from the input signal is: X[I%-J%]. Lines 210 and 220 prevent this from being outside the defined array, X[0] to X[80]. In other words, this program handles undefined samples in the input signal by ignoring them. Another alternative would be to define the input signal's array from X[-30] to X[110], allowing 30 zeros to be padded on each side of the true data. As a third alternative, the FOR-NEXT loop in line 180 could be changed to run from 30 to 80, rather than 0 to 110. That is, the program would only calculate the samples in the output signal where the impulse response is fully immersed in the input signal. The important thing is that you must use one of these three techniques. If you don't, the program will crash when it tries to read the out-of-bounds data. S 0 10 20 30 -0.5 0.0 0.5 1.0 1.5 Sample number 0 10 20 30 40 50 60 70 80 -4 -2 0 2 4 Sample number 0 10 20 30 40 50 60 70 80 90 100 110 -4 -2 0 2 4 Input signal Impulse response Output signal unusable usable unusable Sample number Amplitude Amplitude Amplitude FIGURE 6-10 End effects in convolution. When an input signal is convolved with an M point impulse response, the first and last M-1 points in the output signal may not be usable. In this example, the impulse response is a high-pass filter used to remove the DC component from the input signal. 100 'CONVOLUTION USING THE OUTPUT SIDE ALGORITHM 110 ' 120 DIM X[80] 'The input signal, 81 points 130 DIM H[30] 'The impulse response, 31 points 140 DIM Y[110] 'The output signal, 111 points 150 ' 160 GOSUB XXXX 'Mythical subroutine to load X[ ] and H[ ] 170 ' 180 FOR I% = 0 TO 110 'Loop for each point in Y[ ] 190 Y[I%] = 0 'Zero the sample in the output array 200 FOR J% = 0 TO 30 'Loop for each point in H[ ] 210 IF (I%-J% < 0) THEN GOTO 240 220 IF (I%-J% > 80) THEN GOTO 240 230 Y(I%) = Y(I%) + H(J%) t X(I%-J%) 240 NEXT J% 250 NEXT I% 260 ' 270 GOSUB XXXX 'Mythical subroutine to store Y[ ] 280 ' 290 END TABLE 6-2 122 The Scientist and Engineer's Guide to Digital Signal Processing The Sum of Weighted Inputs The characteristics of a linear system are completely described by its impulse response. This is the basis of the input side algorithm: each point in the input signal contributes a scaled and shifted version of the impulse response to the output signal. The mathematical consequences of this lead to the output side algorithm: each point in the output signal receives a contribution from many points in the input signal, multiplied by a flipped impulse response. While this is all true, it doesn't provide the full story on why convolution is important in signal processing. Look back at the convolution machine in Fig. 6-8, and ignore that the signal inside the dotted box is an impulse response. Think of it as a set of weighing coefficients that happen to be embedded in the flow diagram. In this view, each sample in the output signal is equal to a sum of weighted inputs. Each sample in the output is influenced by a region of samples in the input signal, as determined by what the weighing coefficients are chosen to be. For example, imagine there are ten weighing coefficients, each with a value of onetenth. This makes each sample in the output signal the average of ten samples from the input. Taking this further, the weighing coefficients do not need to be restricted to the left side of the output sample being calculated. For instance, Fig. 6-8 shows y[6] being calculated from: x[3], x[4], x[5], and x[6] . Viewing the convolution machine as a sum of weighted inputs, the weighing coefficients could be chosen symmetrically around the output sample. For example, y[6] might receive contributions from: x[4], x[5], x[6], x[7], and x[8] . Using the same indexing notation as in Fig. 6-8, the weighing coefficients for these five inputs would be held in: h[2], h[1], h[0], h[&1], and h[&2] . In other words, the impulse response that corresponds to our selection of symmetrical weighing coefficients requires the use of negative indexes. We will return to this in the next chapter. Mathematically, there is only one concept here: convolution as defined by Eq. 6-1. However, science and engineering problems approach this single concept from two distinct directions. Sometimes you will want to think of a system in terms of what its impulse response looks like. Other times you will understand the system as a set of weighing coefficients. You need to become familiar with both views, and how to toggle between them. Digital Signal Processors Digital Signal Processing is carried out by mathematical operations. In comparison, word processing and similar programs merely rearrange stored data. This means that computers designed for business and other general applications are not optimized for algorithms such as digital filtering and Fourier analysis. Digital Signal Processors are microprocessors specifically designed to handle Digital Signal Processing tasks. These devices have seen tremendous growth in the last decade, finding use in everything from cellular telephones to advanced scientific instruments. In fact, hardware engineers use "DSP" to mean Digital Signal Processor, just as algorithm developers use "DSP" to mean Digital Signal Processing. This chapter looks at how DSPs are different from other types of microprocessors, how to decide if a DSP is right for your application, and how to get started in this exciting new field. In the next chapter we will take a more detailed look at one of these sophisticated products: the Analog Devices SHARC® family. How DSPs are Different from Other Microprocessors In the 1960s it was predicted that artificial intelligence would revolutionize the way humans interact with computers and other machines. It was believed that by the end of the century we would have robots cleaning our houses, computers driving our cars, and voice interfaces controlling the storage and retrieval of information. This hasn't happened; these abstract tasks are far more complicated than expected, and very difficult to carry out with the step-by-step logic provided by digital computers. However, the last forty years have shown that computers are extremely capable in two broad areas, (1) data manipulation, such as word processing and database management, and (2) mathematical calculation, used in science, engineering, and Digital Signal Processing. All microprocessors can perform both tasks; however, it is difficult (expensive) to make a device that is optimized for both. There are technical tradeoffs in the hardware design, such as the size of the instruction set and how interrupts are handled. Even 504 The Scientist and Engineer's Guide to Digital Signal Processing Data Manipulation Math Calculation Word processing, database management, spread sheets, operating sytems, etc. Digital Signal Processing, motion control, scientific and engineering simulations, etc. data movement (A º B) value testing (If A=B then ...) addition (A+B=C ) multiplication (A×B=C ) Typical Applications Main Operations FIGURE 28-1 Data manipulation versus mathematical calculation. Digital computers are useful for two general tasks: data manipulation and mathematical calculation. Data manipulation is based on moving data and testing inequalities, while mathematical calculation uses multiplication and addition. more important, there are marketing issues involved: development and manufacturing cost, competitive position, product lifetime, and so on. As a broad generalization, these factors have made traditional microprocessors, such as the Pentium®, primarily directed at data manipulation. Similarly, DSPs are designed to perform the mathematical calculations needed in Digital Signal Processing. Figure 28-1 lists the most important differences between these two categories. Data manipulation involves storing and sorting information. For instance, consider a word processing program. The basic task is to store the information (typed in by the operator), organize the information (cut and paste, spell checking, page layout, etc.), and then retrieve the information (such as saving the document on a floppy disk or printing it with a laser printer). These tasks are accomplished by moving data from one location to another, and testing for inequalities (A=B, A