用单片机做FFT
FFT在单片机上的实现
FFT在单片机上的实现摘要音频信号分析仪是一种可广泛见于各种音响、调音和录音设备上的，能实时地采样及分析输入的音频信号的频谱，并将其显示在显示屏上的设备，使人在聆听音乐时能对音乐的高低频能有直观的了解。
本文所介绍的即是这样一个音频信号分析系统。
系统的硬件由信号调理、控制处理器、显示模块三部分组成。
信号调理电路使信号可输入300mV~3V的交流音频信号。
这里只对单路信号处理：当电压较低时使用LM324运放获得增益，对超过12800Hz的信号进行滤波处理。
另外设置输出音频接口以便监听。
控制处理器采用51内核1T单片机STC12C60A5S2，晶振频率为32.768MHz。
该单片机自带8路10位高速ADC，这里只用1路ADC的高8位。
对信号连续采32个点进行浮点型FFT运算。
一次完整采样的时间为1.25ms，最高采样频率为25600Hz，分辨频率为800Hz~12800Hz，分16级。
显示部分主体为1602液晶显示屏，其具有2行×16列的8×5点显示点阵。
16分频谱将分别以柱高形式显示在显示屏上。
程序中设置了频率下落效果以使观感更好。
另设置了对比度调节电阻，使屏幕对比度可调。
关键词：FFT 单片机音频频谱THE REALIZATION OF FFT IN THEMICROCONTROALERABSTRACTAudio signal analyzer is a kind of device which can be widely found in various of audio, mixing and recording devices, and can sampling and analysis of the spectrum of the input have an intuitive audio signal and displays it on the display in real-time, people can are listening to Music for music when high frequency understanding. What presented in this article is just such an audio signal analysis system. The hardware of the system are formed with three parts: the signal conditioner, the control processor and the display module.Signal conditioning circuit makes the signal of 300mV ~ 3V AC audio signal available for inputting. In this system,we only process with single-channel signal: When the voltage is lower the system uses LM324 op amp to gain voltage, and as to signals more than 12800Hz it filters them. In addition the system sets an output audio interface for monitoring.The control processor of the system is the 51 cores 1T MCU STC12C60A5S2, with 32.768MHz crystal frequency. The device comes with 8-channel &10-bit high-speed ADC, where only one channel ADC high 8. The signal collected 32 points in consecutive floating-point FFT operation. A complete sampling time is 1.25ms, the maximum sampling frequency is 25600Hz, and the distinguish frequency is 800Hz ~ 12800Hz, with 16 levels.The main display section is 1602 LCD screen, which has 2 rows ×16 columns - 8 × 5 dot display matrix. 16 points to the column height spectrum will be displayed on the display. The process of setting of the frequency drop in the perception of better effect. There is also contrast adjustment resistor, which makes the screen contrast adjustable.KEY WORDS：FFT,MCU,AUDIO SPECTRUM目录第1章绪论 (1)§1.1研究的背景及意义 (1)§1.1.1课题研究背景 (1)§1.1.2课题研究意义 (1)§1.2课题发展的状况 (1)§1.3设计任务 (2)第2章系统方案设计 (3)§2.1 系统方案设计 (3)§2.2系统硬件的选择 (3)§2.2.1处理器的比较与选择 (3)§2.2.2采样模块的确定 (4)§2.2.3显示器件的比较和选择 (4)第3章系统硬件设计 (5)§3.1 单片机STC12C5A60S2 (5)§3.1.1单片机STC12C5A60S2功能简介 (5)§3.1.2 单片机STC12C5A60S2引脚图 (6)§3.1.2 单片机的最小系统 (6)§3.2 显示屏LCD1602 (7)§3.2.1 LCD1602简介 (7)§3.2.2 LCD1602的硬件连接 (8)§3.3.1 LM324电压增益与偏移电路 (9)§3.3.2 滤波电路 (10)第4章系统软件设计 (11)§4.1系统软件总体设计 (11)§4.2 系统软件详细设计 (12)§4.2.1 系统的准备和初始化 (12)§4.2.2 AD采样子程序 (13)§4.2.3 蝶形运算的FFT算法 (15)§4.2.4 显示子程序 (17)第5章系统调试 (20)§5.1 信号电压调试 (20)§5.2 单频率信号测试 (20)§5.2.1 实际频率分度测试 (21)§5.2.1 频率混叠和滤波效果 (22)§5.3 实际使用效果 (22)结论 (24)参考文献 (25)致谢 (26)附录 (27)一、主程序代码 (27)二、原理图 (35)第1章绪论§1.1研究的背景及意义§1.1.1课题研究背景在家庭影院、卡拉OK等音响系统中，实时显示音乐信号的频谱将为音响系统增不少色彩。
单片机实现傅里叶变换
单片机实现傅里叶变换单片机是一种集成了微处理器、内存和输入输出设备等功能于一体的微型计算机系统。
在工程实践中，单片机广泛应用于各种控制系统中，包括自动化控制、仪器仪表控制、电力电子控制等。
傅里叶变换是一种重要的数学工具，可以将时域信号转换为频域信号，对信号的频谱特性进行分析。
本文将探讨如何利用单片机实现傅里叶变换。
傅里叶变换的基本原理是将一个周期性函数分解为一系列正弦函数的叠加。
在数字信号处理中，傅里叶变换可以通过离散傅里叶变换（DFT）来实现。
DFT是一种将离散信号转换为频域信号的方法，可以将时域上的数字序列转换为在频域上的能量谱密度。
要实现傅里叶变换，首先需要将输入信号进行采样。
采样是将连续信号离散化的过程，可以通过模数转换器（ADC）将模拟信号转换为数字信号。
在单片机中，可以通过ADC模块来实现信号的采样。
接下来，需要对采样信号进行存储和处理。
在单片机中，可以利用内存来存储采样信号，并利用处理器对信号进行处理。
通过计算，可以将离散信号转换为频域信号。
在单片机中，可以利用快速傅里叶变换（FFT）算法来实现傅里叶变换。
FFT算法是一种高效的计算DFT的方法，可以大大减少计算复杂度。
通过使用FFT算法，可以在较短的时间内完成傅里叶变换。
在单片机中，可以通过软件或硬件来实现FFT算法。
软件实现需要编写相应的程序代码来实现FFT算法，而硬件实现则可以利用专门的FFT芯片或者FPGA来加速计算。
除了傅里叶变换，单片机还可以实现其他的频域分析方法，如离散余弦变换（DCT）、小波变换等。
这些方法在不同的应用领域中有着广泛的应用，可以对信号进行更加深入的分析。
在实际应用中，单片机实现傅里叶变换可以用于音频信号处理、图像处理、通信系统等领域。
例如，在音频信号处理中，可以通过傅里叶变换来实现音频信号的频谱分析、滤波等操作。
在图像处理中，可以利用傅里叶变换来进行图像增强、去噪等操作。
在通信系统中，可以利用傅里叶变换来进行信号调制、解调等操作。
C8051F系列单片机的FFT实现
Rev. 1.1 12/03Copyright © 2003 by Silicon LaboratoriesAN142-DS11AN142FFT R O U T I N E S F O R T H E C8051F12X F A M I L Y Relevant DevicesThis application note applies to the following devices:C8051F124, C8051F125, C8051F126, and C8051F127.IntroductionThe Fast Fourier Transform (FFT) is an efficient method for calculating the Discrete Fourier Trans-form (DFT) of a signal. This note provides a brief introduction to the FFT, and describes two example FFT routines written in ‘C’ that have been opti-mized for execution time and RAM storage space on Silicon Labs microcontrollers. The example routines use the 10 or 12-bit ADC to collect theinput data for the FFT routine, and the results are sent out through the UART, where they can be dis-played using terminal software on a PC.Only the very basic aspects of the FFT that are nec-essary to describe the algorithms are presented here. A more detailed explanation of the DFT and the FFT can be found in References [1] and [2].Radix-2 FFT AlgorithmsThe output of the FFT is identical to the output of the DFT, but a number of redundant calculations have been eliminated to allow for faster computa-tion. For an N-point DFT, the required number of complex multiplications is N 2. For an N-point FFT,the number of complex multiplications required is:This optimization leads to a drastic speed improve-ment over the DFT as N becomes large. For exam-ple, a 64-point DFT requires 4096 complexmultiplications, while the corresponding FFT requires only 192. Using this and other optimiza-tions, an FFT can be calculated in a relatively short amount of time on a Silicon Labs 8051 processor.The FFT routines presented in this note are both Radix-2 Decimation-in-Time algorithms. Radix-2algorithms operate by separating the original DFT into a number of 2-point DFT computations. First,the original N-point DFT is split into two DFTs of N/2 points each. The resulting N/2-point DFTs are then each split into two N/4-point DFTs, and so on,until the number of points in each smaller DFT is reduced to two. This method requires that the FFT size be a power of two.The basic 2-point DFT performed in the Radix-2Decimation-in-Time algorithm is shown in Figure 1. This structure, named a “butterfly”, isused to perform all of the computations necessaryfor the FFT. The inputs (A and B) and the outputs (A’ and B’) of the butterfly are complex numbers containing the data that is being processed. W n rep-resents a complex sinusoidal value that is applied ateach stage of the FFT.N2---N 2log⋅Figure 1. Radix-2 Decimation-in-TimeButterfly StructureAN1422Rev. 1.1Index Bit ReversalThe FFT algorithms presented here are performed on the data in-place, to minimize the amount of temporary storage space required for intermediate data. To perform these algorithms in-place, either the input data or the output data of the FFT routine will be sorted in bit reversed order. To change between normal order and bit reversed order, each data point is swapped with another location in the data set determined by reversing the order of the bits in the sample index. For example, in a 16-point FFT, the sample stored at index 0001b (1 decimal) would swap locations with the sample stored at index 1000b (8 decimal). Locations where the bit reversed index are equal to the not bit-reversed index, such as 0110b (6 decimal) are not swapped. The order of operations in the FFT computation is determined by whether the inputs or the outputs of the FFT are sorted in bit-reversed order. Windowing Input DataThe FFT algorithm operates on a data set that rep-resents a finite length of time, but makes the assumption that the data set is periodic and repeated infinitely. When the sample set is repeated in this way, the last sample (index[N-1]) is adja-cent to the first sample (index[0]). As shown in Figure2, this can lead to a discontinuity in the sig-nal that the FFT “sees” when the data is not peri-odic over the sample set. Because of this, data is normally windowed before it is processed by an FFT routine. Windowing makes the data periodic over the sample set and removes any discontinuity between the first and last samples in the set. Because windowing changes the input data set, it produces some artifacts in the frequency domain. Windowing “spreads” the signal energy among multiple bins, as shown in Figure3. This energy spreading has the effect of attenuating the peak value of the signal. Most of the signal’s original content is stored in the “main lobe”, while a small amount leaks into the “side lobes”. The width of the main lobe and the height of the side lobes are dependent on the window algorithm that is applied to the data. Some common windows and their prop-erties are summarized in Table1 . Some equations for computing the window coefficients for an N-point FFT are listed in Table2 . More extensive information on window algorithms and their parameters can be found in Reference [3].AN142Rev. 1.13Figure 2. Time Domain WindowingAN1424Rev. 1.1Table 1. Window PropertiesMain Lobe Width Processing Gain Side Lobe Height None (Rectangular) 1 Bin 1.0-13 dBTriangular 3 Bins0.5-27 dBHanning 3 Bins0.5-32 dBHamming 3 Bins0.54-43 dBBlackman 5 Bins0.42-58 dBTable 2. Window Coefficient EquationsWindow EquationNone (Rectangular)W(n) = 10 ≤ n < NTriangular W(n) = n / (N / 2)W(n) = 2 - n / (N / 2)0 ≤ n ≤ N / 2 N / 2 < n < NHanning W(n) = 0.5 - 0.5cos(2πn / N)0 ≤ n < N Hamming W(n) = 0.54 - 0.46cos(2πn / N)0 ≤ n < NBlackman W(n) = 0.42 - 0.5cos(2πn / N) +0.08cos(4πn / N)0 ≤ n <N Figure 3. Window Effects in the Frequency DomainAN142Rev. 1.15Software ExamplesThis application note contains two FFT routines that are very similar. Both examples collect data from the ADC, window the collected data, perform an FFT, and print the results to the UART running at 115200 baud. The difference between the two example files is whether the inputs or the outputs of the FFT routine are in bit-reversed order. The first routine, “IntFFT_BRIN.c”, accepts input data sorted in bit-reversed order and produces outputsorted in normal order. The second routine,“IntFFT_BROUT.c”, accepts data sorted in normalorder and produces output sorted in bit-reversed order. The two examples produce identical results, but each has its advantages. Performing a bit-reversal sort takes some processing time and is not always necessary for every application. For example, an application that is performing peak frequency detection needs only to know which output bin of the FFT has the highest magnitude. Once this is determined, the frequency for that particular bin can be calculated by bit reversing only the bin of interest. In this example, the routine that accepts input data in normal order and produces outputs in bit-reversed order (IntFFT_BROUT.c) can be used without having to re-sort the output data after-wards.If the application requires that the output of the FFT appears in normal order, and the input to the routine is a real signal (not a complex input), the routine that requires bit-reversed input data (IntFFT_BRIN.c) may be a better choice. The rea-son for this is that the imaginary locations in the real input set can be assumed to be zero, and some processor time can be saved by only sorting the real data. In this example, the complex output of the FFT appears in normal order.Parameter SpecificationThe file “FFT_Code_Tables.h” contains code tables for FFT sizes from 4 points to 1024 points,and four different window types. Conditional com-pilation is used so that only the necessary tables for the desired FFT size and window type are included.The parameter NUM_FFT should contain the size of the FFT to perform. Valid FFT sizes for this example code are 4, 8, 16, 32, 64, 128, 256, 512,and 1024 points. The parameter WINDOW_TYPE should contain a number from 0 to 4, which speci-fies the desired window algorithm, according to Table 3 .One other parameter that affects the program oper-ation is the RUN_ONCE definition, which is con-tained in the main source file. If RUN_ONCE is zero, the program will continue to collect new data sets and perform FFTs on them. If RUN_ONCE is a non-zero value, the program will stop after one iter-ation.Data CollectionADC0 collects the data to process using Timer 3 as a start-of-conversion source. ADC0 is configured to sample at a speed determined by the SAMPLE_RATE constant, in single-ended mode on channel AIN0.0. Samples are collected and stored in the Real[] array. The 12-bit ADC data is left-jus-tified and stored as 16-bit data with trailing zeros.Once NUM_FFT samples have been collected,ADC0 interrupts are turned off, and the data is win-dowed.Table 3. Window Selection 0None (Rectangular)1Triangular 2Hanning 3Hamming 4BlackmanAN1426Rev. 1.1Window ImplementationAfter the data has been collected, it is windowed using the window selected by the WINDOW_TYPE constant (see Table3 ). The window information is stored in code space as a table of unsigned integer (16-bit) values. The integers represent a fractional number which can be computed by dividing the stored integer by 65536. For example, the value 32768 represents a multiplication value of 32768 / 65536, or 0.5. Because the windows are all sym-metrical about NUM_FFT/2, the window tables only contain values for one half of the window to save storage space. The window value at index NUM_FFT/2 is assumed to be equal to 1.0, and is not stored. The WindowCalc() function performs a multiplication of each input sample with its corre-sponding value in the window table. Table4 shows how the window table is indexed when multiply-ing.The WindowCalc() function also changes data that has been stored in single-ended (unsigned) format into differential (signed 2’s complement) format. This centers the data about 0x0000 to remove the DC bias. When the SE_data input variable is non-zero, each sample is XORed with the value 0x8000. This inverts the MSB of the data, which has the same effect as subtracting 0x8000 from all samples. If the SE_data variable is zero, the data is assumed to already be in 2’s complement format, and is not changed.FFT OptimizationsA number of optimizations were implemented in the example routines. The primary goals of these optimizations were to maximize the speed of the routine and to minimize the amount of RAM needed for data storage. The specifics of each opti-mization are detailed in the sections that follow. Integer Storage and Computation The FFT algorithm has an inherent processing gain. For a real input, the processing gain is equal to N/2, where N is the number of points in the FFT (the processing gain for complex inputs is equal to N). Because of this gain, more bits are required to store the output information of the FFT than are required to store the input information. Instead of using additional space, the example FFT routines use 16-bit signed integer values to store all input, output, and intermediate data. To compensate for the FFT processing gain, and to store all data as 16-bit num-bers, the computed values are divided by 2 after each butterfly calculation. Using this method, the overall gain of the FFT routine, not including win-dow gain, is 1/2 for a real input, and 1 for a com-plex input.Integer math is used to perform all calculations. This improves the speed of the FFT routine over fixed-point or floating-point math. Multiplication is avoided whenever possible, and all divide opera-tions are implemented using right shifts. Whenever this type of divide is performed, the result is rounded to counter the asymmetrical effects of truncating a 2’s complement number. If the number is negative, and the deleted bits are non-zero, a one is added to the result so that both negative and pos-itive numbers are always rounded towards zero. Sinusoid Table StorageTo reduce computation time, sine and cosine values are not calculated in real-time. Instead, they are pre-calculated and stored in code space. The SinTable[] array declared in “FFT_Code_Tables.h”Table 4. Window Index Decoding0 to (N / 2 - 1)MN / 2No Multiplication (N / 2 + 1) to (N -1)N - MN is represented by NUM_FFT in the software examples.AN142Rev. 1.17contains the data for 1⁄4 period of a sine function.The values stored in the SinTable[] array are 2’s complement integer (16-bit) values. The integers represent fractional numbers between -1.0 and 1.0(though negative values and 1.0 are not stored).The actual multiplication value can be computed by dividing the integer value by 32768. For exam-ple, a value of 8192 in the SinTable[] array repre-sents a multiplication value of 8192⁄32768, or 0.25. The FFT routine needs sine and cosine values for 1⁄2 of the sinusoid period. The 1⁄4 sine wave stored in code space is indexed according to Table 5 to generate the necessary values.Real Input OnlyThe example routines have been optimized to com-pute an FFT on a real input. During the first stage of the FFT, it is assumed that all imaginary loca-tions are equal to zero, which eliminates a number of calculations. The real input is modified as needed, and the imaginary locations are set to zero.The FFT examples can easily be modified to oper-ate on complex input data, as detailed in "Using Complex Inputs‚" on page 11.Bit Reversal SortingInstead of fully computing a bit-reversed value for each index, the Bit_Reverse() function uses the BRTable[] array to look up the bit-reversed index locations. Bit-reversed index values for the first half of the indices (0 through NUM_FFT /2-1)are stored in the table directly. Bit-reversed index values for the second half of the indices (NUM_FFT /2 through NUM_FFT -1) are calcu-lated by adding 1 to the table value stored at Index -NUM_FFT /2. For example, the bit-reversed value for index NUM_FFT /2 is calcu-lated as:BRTable [NUM_FFT /2-NUM_FFT /2]+1.Unnecessary MultiplicationsEach iteration of the butterfly calculation requires one complex multiplication (a total of four long integer multiplications). Long integer multiplica-tions take many processor cycles to complete.Many of these multiplications are not necessary and can be eliminated. Specifically, when multiply-ing by zero, the result will be zero, and when multi-plying by one, the result will be the original number. The example code checks for the cases where the sine and cosine values are equal to zero or one, and takes these shortcuts when it can. TheTable 5. Sine Table Index DecodingM = 0010 < M < N / 4SinTable[M]SinTable[(N / 4) - M]M = N / 41N / 4 < M < N / 2SinTable[(N / 2) - M]SinTable[M - (N / 4)]M and N are represented by the following parameters in the software examples:N = NUM_FFT M = sin_indexAN1428Rev. 1.1number of complex multiplications saved with this optimization is equal to:where N is the number of points in the FFT.N 2---2ii 1=N ()1–2log ∑+AN142Rev. 1.19PerformanceThe approximate number of clock cycles the processor requires for each of the key routines in the exam-ple code is shown in Table 6 and Table 7 , for FFT sizes of 16, 64, 256, and 1024 points. For larger FFT sizes, the number of clock cycles required for the WindowCalc() and BitReverse() routines increases lin-early, while the number of clock cycles required for the IntFFT() routine increases according to the num-ber of complex multiplications required, as discussed in "Radix-2 FFT Algorithms‚" on page 1. The total execution time for these routines when the controller is operating with a 49.7664MHz clock is shown in Table 8 .Table 6. Approximate Timing in SYSCLK Cycles (IntFFT_BRIN.c)166,5001,60038,00046,0006427,0007,500260,000294,000256113,00030,000 1.5 Million 1.64 Million 1024452,000145,0007.8 Million8.4 MillionTable 7. Approximate Timing in SYSCLK Cycles (IntFFT_BROUT.c)166,5002,90038,00047,0006427,00013,000260,000299,000256113,00055,000 1.5 Million 1.66 Million 1024452,000244,0007.8 Million8.5 MillionTable 8. Total Time With a 49.7664 MHz System Clock16924 µs 945 µs 64 5.9 ms 6 ms25633 ms 33.4 ms 1024169 ms171 msAN14210Rev. 1.1Code and XRAM SizeTable9 and Table10 list the code space requirements for the key functions and code tables used in the example software. The total amount of XRAM required for the storage of the Real[] and Imag[] integer arrays is equal to the FFT size times 4. For example, a 1024-point FFT requires 4096 bytes (4 x 1024) of XRAM storage space.Table 9. Code Space Requirements in Bytes (IntFFT_BRIN.c)4560195229183054 85601972291163064 165601972291323080 325601972291643112 6456019722911283176 12856019722912563304 25661719722915123617 512619234241710244294 1024619252241925605850Table 10. Code Space Requirements in Bytes (IntFFT_BROUT.c)4560353229283213 85603552292163223 165603552292323239 325603552292643271 6456035522921283335 12856035522922563463 25661735522925123776 512619406242610244475 1024619424243725606040Possible ModificationsThe FFT routines can be expanded or optimized further depending on the application. The follow-ing sections describe a few of the possible modifi-cations.FFT SizeThe example routines are limited to between 4 and 1024 samples with the included sinusoid, window, and bit reversal tables. Larger tables can be defined to allow the software to perform larger FFTs. As written, the routines will perform FFTs of up to 32768 points, provided that there is enough code space to store the tables and enough RAM to store the Real[] and Imag[] arrays.Using Complex InputsThe example sfotware is optimized for real input data only. With some minor changes, the FFT rou-tines can be used to compute the FFT of complex data sets as well. The following things must be changed to perform FFTs of complex input data: 1.Remove or comment out the section of codelabeled “FIRST STAGE” in the Int_FFT() rou-tine.2.For “IntFFT_BRIN.c”, the Bit_Reverse() rou-tine should be changed to perform a bit reversal on both the Real and Imaginary data. The Bit_Reverse() routine in “IntFFT_BROUT.c”does the bit reversal sort on both arrays, and can be used here. Alternately, the Bit_Reverse() function in “IntFFT_BRIN.c” can be called twice: once to sort the Real[] array, and once to sort the Imag[] array.3.For “IntFFT_BRIN.c”, the variable group inthe function Int_FFT() should be initialized to ‘1’, and the variable stage should be initialized to NUM_FFT⁄2.4.For “IntFFT_BROUT.c”, the variable group inthe function Int_FFT() should be initialized to NUM_FFT⁄2, and the variable stage should be initialized to ‘1’.5.If the data is to be windowed, the function Win-dowCalc() must be modified to window the Imaginary data as well as the Real data. Alter-nately, the WindowCalc() function can be called twice: once to window the Real[] array, and once to window the Imag[] array. Increasing Output PrecisionA more involved software modification would entail changing the way that data is computed and stored. To help maximize the speed of the FFT, the routines presented here store data as 16-bit integer values. This limits the output precision of the FFT. If larger variables or floating-point numbers are used to store the data, the output precision of the routine can be increased to make it more suitable for measurements such as Signal-to-Noise or Sig-nal-to-Distortion.References[1]Lyons, R. G. Understanding Digital Signal Pro-cessing, Addison Wesley Longman, Inc., Reading, Massachusetts, 1997, pp. 49-154.[2]Oppenheim, A. V. and Schafer, R. W. Discrete-Time Signal Processing, 2nd Ed. Prentice Hall, Inc., Upper Saddle River, New Jersey, 1999, pp. 541-669.[3]Harris, F. J. “On the Use of Windows for Har-monic Analysis with the Discrete Fourier Trans-form,” Proceedings of the IEEE, V ol. 66, No. 1, January 1978.IntFFT_BRIN.c//-----------------------------------------------------------------------------// IntFFT_BRIN.c//-----------------------------------------------------------------------------// Copyright 2003 Cygnal Integrated Products, Inc.//// AUTH: BD// DATE: 30 JAN 03//// This program collects data using ADC0 at <SAMPLE_RATE> Hz and performs// an FFT on the data. The Real and Imaginary parts of the results are then// sent to the UART peripheral at <BAUDRATE> bps, where they can be displayed // or captured using a terminal program.//// Note that the FFT performed in this software is optimized for storage space // (RAM). The resulting Frequency-domain data is not suitable for analyzing// Signal-to-noise or distortion performance.//// This program uses a 22.1184 MHz crystal oscillator multiplied by (9/4)// for an effective SYSCLK of 49.7664 Mhz. This program also initializes and// uses UART0 at <BAUDRATE> bits per second.//// Target: C8051F12x// Tool chain: KEIL C51 6.03////-----------------------------------------------------------------------------// Includes//-----------------------------------------------------------------------------#include <c8051f120.h> // SFR declarations#include <stdio.h>#include “FFT_Code_Tables.h” // Code Tables for FFT routines//-----------------------------------------------------------------------------// 16-bit SFR Definitions for ‘F12x//-----------------------------------------------------------------------------sfr16 DP = 0x82; // data pointersfr16 ADC0 = 0xbe; // ADC0 datasfr16 ADC0GT = 0xc4; // ADC0 greater than windowsfr16 ADC0LT = 0xc6; // ADC0 less than windowsfr16 RCAP2 = 0xca; // Timer2 capture/reloadsfr16 RCAP3 = 0xca; // Timer3 capture/reloadsfr16 RCAP4 = 0xca; // Timer4 capture/reloadsfr16 TMR2 = 0xcc; // Timer2sfr16 TMR3 = 0xcc; // Timer3sfr16 TMR4 = 0xcc; // Timer4sfr16 DAC0 = 0xd2; // DAC0 datasfr16 DAC1 = 0xd2; // DAC1 datasfr16 PCA0CP5 = 0xe1; // PCA0 Module 5 capturesfr16 PCA0CP2 = 0xe9; // PCA0 Module 2 capturesfr16 PCA0CP3 = 0xeb; // PCA0 Module 3 capturesfr16 PCA0CP4 = 0xed; // PCA0 Module 4 capturesfr16 PCA0 = 0xf9; // PCA0 countersfr16 PCA0CP0 = 0xfb; // PCA0 Module 0 capturesfr16 PCA0CP1 = 0xfd; // PCA0 Module 1 capture//-----------------------------------------------------------------------------// Global CONSTANTS and Variable Type Definitions//-----------------------------------------------------------------------------#define NUM_BITS 16 // Number of Bits in Data#define DATA_BEGIN 0x0000 // Beginning of XRAM Data#define EXTCLK 22118400 // External oscillator frequency in Hz #define SYSCLK 49760000 // Output of PLL derived from// (EXTCLK*9/4)#define BAUDRATE 115200 // Baud Rate for UART0#define SAMPLE_RATE 10000 // Sample frequency in Hz#define RUN_ONCE 1 // Setting to a non-zero value will// cause the program to stop after one // data set.typedef union IBALONG { // Integer or Byte-addressable LONGlong l; // long: Var.lunsigned int i[2]; // u int: Var.i[0]:Var.i[1]unsigned char b[4]; // u char: Var.b[0]:Var.b[1]:// Var.b[2]:Var.b[3]} IBALONG;typedef union BAINT { // Byte-addressable INTint i; // int: Var.iunsigned char b[2]; // u char: Var.b[0]:Var.b[1]} BAINT;//-----------------------------------------------------------------------------// Function PROTOTYPES//-----------------------------------------------------------------------------void WindowCalc(int Win_Array[], unsigned char SE_data);void Int_FFT(int ReArray[], int ImArray[]);void Bit_Reverse(int BR_Array[]);void SYSCLK_Init (void);void PORT_Init (void);void UART0_Init (void);void ADC0_Init (void);void TIMER3_Init (int counts);void ADC0_ISR (void);//-----------------------------------------------------------------------------// Global Variables//-----------------------------------------------------------------------------// XRAM storage of FFT: requires NUM_FFT*4 Bytes after DATA_BEGIN addressint xdata Real[NUM_FFT] _at_ DATA_BEGIN;int xdata Imag[NUM_FFT] _at_ (DATA_BEGIN + (NUM_FFT * 2));// NUM_FFT is defined in the “FFT_Code_Tables.h” header file#if (NUM_FFT >= 256)unsigned int index, ADC_Index;#endif#if (NUM_FFT < 256)unsigned char index, ADC_Index;#endifunsigned int BinNum;bit Conversion_Set_Complete; // This indicates when the data has been // stored, and is ready to be processed // using the FFT routines//-----------------------------------------------------------------------------// MAIN Routine//-----------------------------------------------------------------------------void main(){// disable watchdog timerWDTCN = 0xde;WDTCN = 0xad;SYSCLK_Init(); // initialize external clock and PLLPORT_Init (); // set up Port I/OUART0_Init (); // initialize UART0TIMER3_Init (SYSCLK/SAMPLE_RATE); // initialize Timer3 to overflow at// <SAMPLE_RATE>ADC0_Init (); // init ADC0EA = 1; // globally enable interruptswhile (1){ADC_Index = 0;Conversion_Set_Complete = 0;EIE2 |= 0x02; // enable ADC interruptsSFRPAGE = LEGACY_PAGE;while(!Conversion_Set_Complete);SFRPAGE = UART0_PAGE;printf(“\nCollected Data\nSample\tValue\n”);for (BinNum = 0; BinNum < NUM_FFT; BinNum++){// Print Data in the format: Sample <tab> Value <tab>printf(“%d\t%u\n”, BinNum, Real[BinNum]);}WindowCalc(Real, 1); // Window Real Data, and convert to// differential if it is single-endedBit_Reverse(Real); // Sort Real (Input) Data in bit-reverse // orderInt_FFT(Real, Imag); // Perform FFT on dataSFRPAGE = UART0_PAGE;printf(“\nBin\tReal\tImag\n”);// Output the FFT data to the UARTfor (BinNum = 0; BinNum < NUM_FFT; BinNum++){// Print Data in the format: Bin <tab> Real <tab> Imaginaryprintf(“%d\t%d\t%d\n”, BinNum, Real[BinNum], Imag[BinNum]);}if (RUN_ONCE)while(1);}} // END MAIN//-----------------------------------------------------------------------------// WindowCalc//-----------------------------------------------------------------------------//// Uses the values in WindowFunc[] to window the stored data.//// The WindowFunc[] Array contains window coefficients between samples// 0 and (NUM_FFT/2)-1, and samples from NUM_FFT/2 to NUM_FFT-1 are the mirror // image of the other side.// Window values are interpreted as a fraction of 1 (WindowFunc[x]/65536).// The value at NUM_FFT/2 is assumed to be 1.0 (65536).//// If SE_data = 1, the input data is assumed to be single-ended, and is// converted to a 2’s complement, differential representation, to cancel the DC // offset.//void WindowCalc(int Win_Array[], unsigned char SE_data){#if (WINDOW_TYPE != 0) // Use this section if a window has been specified IBALONG NewVal;if (SE_data) // If data is single-ended,Win_Array[0] ^= 0x8000; // convert it to differentialNewVal.l = (long)Win_Array[0] * WindowFunc[0];if ((NewVal.l < 0)&&(NewVal.i[1]))Win_Array[0] = NewVal.i[0] + 1;else Win_Array[0] = NewVal.i[0];if (SE_data) // If data is single-ended,Win_Array[NUM_FFT/2] ^= 0x8000; // convert it to differentialfor (index = 1; index < NUM_FFT/2; index++){// Array positions 1 to (NUM_FFT/2 - 1)if (SE_data) // If data is single-ended,Win_Array[index] ^= 0x8000; // convert it to differentialNewVal.l = (long)Win_Array[index] * WindowFunc[index];if ((NewVal.l < 0)&&(NewVal.i[1]))Win_Array[index] = NewVal.i[0] + 1;。
单片机fft运算详解
单片机fft运算详解摘要：I.单片机简介- 定义和特点- 应用领域II.FFT 运算介绍- FFT 的含义和作用- FFT 算法的基本原理III.单片机FFT 运算详解- 单片机FFT 运算的基本流程- 单片机FFT 运算的实现方法- 单片机FFT 运算的优缺点分析IV.单片机FFT 运算的应用- 在信号处理领域的应用- 在图像处理领域的应用V.单片机FFT 运算的发展趋势- 技术的发展- 应用的拓展正文：单片机是一种集成电路，具有一个或多个处理器内核，用于控制和执行特定任务。
单片机的特点是体积小、成本低、功耗低、功能强大，因此被广泛应用于各个领域，如家电、工业控制、通信等。
在数字信号处理领域，快速傅里叶变换(FFT) 是一种重要的算法。
它可以在短时间内将信号从时域转换到频域，进而分析信号的频率成分和能量分布。
FFT 算法的基本原理是将信号序列分解成一组正交基函数的线性组合，从而实现频域分析。
单片机FFT 运算是一种基于单片机的FFT 算法，它可以在单片机上实现FFT 运算。
单片机FFT 运算的基本流程包括输入信号的预处理、FFT 算法的实现和输出结果的解析。
在实现单片机FFT 运算时，通常采用蝶形算法或旋转基2 算法等高效算法。
单片机FFT 运算的优点是运算速度快、精度高、功耗低，适用于实时信号处理和低功耗应用。
缺点是运算量较大，对单片机的处理能力和存储容量要求较高。
单片机FFT 运算在信号处理和图像处理领域有广泛的应用。
例如，在通信系统中，它可以实现数字信号的快速傅里叶变换，从而实现数字信号的调制和解调；在图像处理领域，它可以实现图像的频域滤波和边缘检测，从而提高图像的质量和清晰度。
单片机fft计算频率
单片机fft计算频率在单片机上进行FFT（快速傅里叶变换）可以用于计算信号的频谱，从而分析信号中的不同频率分量。
在单片机上执行FFT通常需要一些数学库或专门的FFT库，因为FFT是一个计算密集型的算法，需要较多的计算资源。
以下是在单片机上执行FFT的一般步骤：选择FFT库：选择适合你的单片机架构的FFT库。
一些流行的嵌入式系统FFT库包括CMSIS-DSP库、FFTW（Fastest Fourier Transform in the West）等。
这些库提供了实现FFT所需的函数和算法。
获取信号数据：从传感器、模拟输入或其他来源获取信号数据。
这可以是一维的时间域信号，比如声音或振动信号。
初始化FFT库：使用FFT库的初始化函数进行设置，包括FFT 的大小和其他参数。
输入信号数据：将获取到的信号数据输入到FFT库中。
这可能涉及到对数据进行预处理，例如去直流分量、窗函数处理等。
执行FFT：调用FFT库的计算函数执行FFT算法。
获取频谱数据：从FFT库中获取计算得到的频谱数据。
频谱数据表示信号在不同频率上的幅度。
分析频谱数据：使用频谱数据进行进一步的分析。
你可以识别信号中的主要频率分量，了解信号的频谱特性。
下面是一个简单的例子，使用CMSIS-DSP库在ARM Cortex-M微控制器上执行FFT的步骤：cCopy code#include "arm_math.h"#define FFT_SIZE 1024#define SAMPLING_FREQ 1000// 定义信号数组float32_t signal[FFT_SIZE];// 定义FFT输出数组float32_t fftOutput[FFT_SIZE * 2];int main(void) {// 初始化CMSIS-DSP库SystemInit();// 初始化FFT库arm_cfft_radix4_instance_f32 fftInstance;arm_cfft_radix4_init_f32(&fftInstance, FFT_SIZE, 0, 1);// 获取信号数据（示例中使用随机数据）for (int i = 0; i < FFT_SIZE; i++) {signal[i] = rand() % 100;}// 执行FFTarm_cfft_radix4_f32(&fftInstance, signal);// 获取频谱数据arm_cmplx_mag_f32(signal, fftOutput, FFT_SIZE);// 在这里进行进一步的频谱数据分析while (1) {// 主循环}}这只是一个简单的例子，实际上需要根据具体的硬件和需求进行更详细的配置和处理。
基于单片机的快速傅里叶算法实现
１ 音 乐 频谱 显 示 介绍
的平方次复数乘法和 Ｎ ／ ２ （ Ｎ ／ ２ — １ ） 次复数加法 。
由此 可 见 ， 仅 仅 经 过 一 次 分 解 ，就 使 运 算 量 减 少 近 一 半 ，所 以可 以对 Ｎ／ ２点 ＤＦ Ｔ 作 进 一 步
数 据采 集测 试：将 音频 信号 输入 至单 片 机 ，再通过 串口观察数据 的变化，观 察工具为
系统进行测试 以保证方案 的可 行性 。测试 可分
Ｆ Ｆ Ｔ算 法 测 试 ： 该 测 试 为 验 证 Ｆ Ｆ Ｔ算 法
两次复数加法。经过一次分解后，计算 １ 个Ｎ 为Ｆ Ｆ Ｔ 算法的测试 、 数据采 集测试和整体测试 。
得 出的结果 。主要测试环境为 Ｍｉ ｃ ｒ ｏ ｓ ｏ ｆ ｔ Ｖ ｉ ｓ ｕ ａ ｌ
２ 硬 件 设 计
其大 体设 计思 路为 ：首先 对音 频信 号进 行采集 ，经过数据处理后输 出至端 口，最后驱 动屏幕显示。最小系统 ，是指芯片能够正常工
作 所 需 的最 小 电路 。在 最 小 系 统 的基 础 上 进 行
பைடு நூலகம்
隔为 ２ 点 的 ２ 个蝶形 。
总 结上述 运算 规律 ，便可 采用 下述 运算 方法 。先从输入端 （ 第 ｌ 级 ）开始 ，逐 级进行 ， 共进 行 Ｍ 级运 算。在进 行第 Ｌ级运 算时 ，依
Ｓ Ｃ Ｍ Ｔ ｅ ｃ ｈ ｎ ｏ ｌ ｏ ｇ ｙ・ 单 片机技术
基于单片机 的快速傅里叶算法 实现
文／ 周 小丽
个管脚 的控制 ，极大地节省 了 Ｉ Ｏ 口的使用 。
本 文 介 绍 了一 种 成 本 相 对 较
（ Ｊ Ｎ ／ ２ ），则 先 将 最 高 位 变 成 ０ （ Ｊ 乍Ｊ — Ｎ／ ２ ），
stm32 傅里叶谐波计算
stm32 傅里叶谐波计算一、傅里叶变换简介傅里叶变换是一种将信号从时域转换到频域的数学方法。
它可以将复杂的信号分解成一系列简单的正弦和余弦函数，即谐波。
这种方法在许多领域都有广泛的应用，如信号处理、图像处理等。
二、STM32硬件介绍STM32是意法半导体（STMicroelectronics）公司推出的一款基于ARM Cortex-M内核的微控制器。
它具有高性能、低功耗、多功能等特点，广泛应用于各种嵌入式系统和物联网领域。
三、傅里叶谐波计算原理傅里叶谐波计算是基于傅里叶变换的理论，通过将信号分解成无数个正弦波和余弦波，从而得到信号的频谱。
在每个频率成分上，可以根据幅度和相位信息还原原始信号。
四、STM32实现傅里叶谐波计算的方法在STM32上实现傅里叶谐波计算，可以采用FFT（快速傅里叶变换）算法。
FFT是一种高效计算离散傅里叶变换（DFT）的算法，时间复杂度为O(n log n)，相较于直接计算DFT的O(n^2)时间复杂度，具有更高的计算效率。
STM32可以通过编程实现FFT算法，或者使用现有的库（如OpenMP、KissFFT等）来进行傅里叶谐波计算。
在实际应用中，首先需要对输入信号进行采样和窗函数处理，以减少频谱泄漏和旁瓣干扰。
接下来，按照FFT算法进行计算，得到信号的频谱。
最后，根据频谱信息进行谐波分析，提取有用信号。
五、应用场景及优势傅里叶谐波计算在许多领域都有广泛的应用，如电力系统中的谐波分析、音频信号处理、图像处理等。
通过STM32实现傅里叶谐波计算，可以充分发挥其高性能、低功耗的特点，满足实时性和精度要求。
六、总结本文介绍了傅里叶变换及其在STM32上的实现方法。
通过傅里叶谐波计算，可以有效地分析信号的频谱成分，为各种应用场景提供有用信息。
STM32-FFT代码说明
FFT 代码说明FFT 为Fast Fourier Transformation ，即快速傅里叶变换，本工程中，FFT 的目标是识别频率为形如式(1.1)的一个正弦信号：()sin(2)f x A ft π=•(1.1)其中，21Hz 5%f =±；因为单片机通过ADC 接口读取该正弦信号的电压值，而12位精度的ADC 的值范围在0-4096之间，如信号经过放大器后映射到之间，那么振幅A 的取值0-4096之间。
假设，信号经过放大器后，其电压值最大为，最小为0V ，那么此信号的振幅为，对应A=2048，即该信号为：()20482048sin(221)f x t π=+•⨯•(1.2)本文中给出的例程即通过FFT 识别式(1.2)这种正弦信号。
假设采样频率为Fs ，信号频率Fn ，采样点数为N 。
那么FFT 之后结果就是一个为N 点的复数。
每一个点就对应着一个频率点。
这个点的模值，就是该频率值下的幅度特性。
具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A ，那么FFT 的结果的每个点〔除了第一个点直流分量之外〕的模值就是A 的N/2倍。
而第一个点就是直流分量，它的模值就是直流分量的N 倍。
第一个点表示直流分量〔即0Hz 〕，而最后一个点N 的再下一个点〔实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后〕那么表示采样频率Fs ，这中间被N -1个点平均分成N 等份，每个点的频率依次增加。
例如某点n 所表示的频率为：()1*/Fn n Fs N =-(1.3)频率分辨率〔/f Fs N ∆=〕等于采样时间的倒数。
例如要分辨，那么需要采集10s 。
1 基于STM32官方DSP 库的FFT 算法工程文件中包含三个函数库，分别为：分别对应数据点数为64,256和1024时的FFT算法。
下面将以数据点数是1024为例，说明FFT的实现过程。
1.1通过函数生成原始数据for(i=0;i<NPT;i++){Fx=2048+2048*sin(PI2*i*21/Fs+20)+1000*sin(PI2*i*15/Fs)+100*sin(PI2*i*25/Fs) +rand()%100;//振幅2048，频率21HZ，为主要的谐波分量，21HZ也是所需信号的频率，rand()%100为0-100之间的随机噪声lBUFIN[i] = ((s16)fx)<<16;//高位为实部，低位为虚部}Fs为采样频率，此处设为Hz，NPT为数据长度〔即前文中提到的N〕，为了能够分辨，数据长度NPT设为1024，那么频率分辨率=102.4/1024=Hz。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
这次的任务是要用单片机做FFT(快速离散傅利叶变换),我觉得单片机越来越无敌了...于是,连续几天被埋在时间与空间的斗争中.因为单片机的硬件资源和条件所限,不能像PC上的MatLab那样方便地进行FFT,我们单片机 (SPCE061A)的时钟频率最快约为50MHz,其RAM 为2K字,ROM为32K字.从前学数据结构和算法的时候,千方百计想着用时间换空间还是用空间换时间,现在才发现,那时是站着说话不腰疼,因为无论拿什么换什么,都是行得通的,损失的效率对于PC应用来说(就我现在的水平),是可以忽略不计的,毕竟不是ACM(此理论的唯一一次失手就是在某次ACM比赛上).整体来说,单片机最大的弱点还是在于速度太慢,根本不可能进行实时FFT,而手动启动来进行运算的速度恐怕也慢得让人难以接受.为了加速,我几天苦苦挣扎于算法的优化上.若不是因为数据精度还有一定的要求(汗~),不能采用8Bit数据,我肯定会将两个数据压缩到一个字(16位)中去以节省空间.坚决榨干单片机的每一滴油!我将所有的冗余计算全部转化为变量表示独立计算然后重复利用,而常量,哪怕是局部的,我也将它写成常数表(虽然还不知道如何在C中将常数表写入ROM区,再汗一个~).现在唯一的希望就在,单片机对浮点数运算的支持能力上...从前碰到浮点数我都会转化为整数,舍不得拿单片机来算浮点数,而这次计算量实在太大了一点儿... 到目前为止,只在GCC上编译通过,并且结果是正确的(与MabLab对照),接下来还是任重而道远啊~贴出主要部分代码:/************************************************************//FFT(快速离散傅利叶变换)//本例为8级256点数据//HUST STI OE0402 POPAPPLE//popapple@//2007-8-9 18:10:20************************************************************/#include "stdio.h"#define LEVEL 8 //级数#define LENGTH 256 //点数//用结构体来表示复数typedef struct _JComplex{float real;float image;}JComplex,*pJComplex;typedef unsigned int uint;//WN为复常数码表,为提高计算速度//其值为exp(-j*2*pi/LENGTH)^n//n为[0..LENGTH/2-1]JComplex WN[LENGTH/2]={{ 1.0000,-0.0000},{ 0.9997,-0.0245},{ 0.9988,-0.0491},{ 0.9973,-0.0736},{ 0.9952 ,-0.0980},{ 0.9925,-0.1224},{ 0.9892,-0.1467},{ 0.9853,-0.1710},{ 0.9808,-0.1951},{ 0.9757 ,-0.2191},{ 0.9700,-0.2430},{ 0.9638,-0.2667},{ 0.9569,-0.2903},{ 0.9495,-0.3137},{ 0.9415 ,-0.3369},{ 0.9330,-0.3599},{ 0.9239,-0.3827},{ 0.9142,-0.4052},{ 0.9040,-0.4276},{ 0.8932 ,-0.4496},{ 0.8819,-0.4714},{ 0.8701,-0.4929},{ 0.8577,-0.5141},{ 0.8449,-0.5350},{ 0.8315 ,-0.5556},{ 0.8176,-0.5758},{ 0.8032,-0.5957},{ 0.7883,-0.6152},{ 0.7730,-0.6344},{ 0.7572 ,-0.6532},{ 0.7410,-0.6716},{ 0.7242,-0.6895},{ 0.7071,-0.7071},{ 0.6895,-0.7242},{ 0.6716 ,-0.7410},{ 0.6532,-0.7572},{ 0.6344,-0.7730},{ 0.6152,-0.7883},{ 0.5957,-0.8032},{ 0.5758 ,-0.8176},{ 0.5556,-0.8315},{ 0.5350,-0.8449},{ 0.5141,-0.8577},{ 0.4929,-0.8701},{ 0.4714 ,-0.8819},{ 0.4496,-0.8932},{ 0.4276,-0.9040},{ 0.4052,-0.9142},{ 0.3827,-0.9239},{ 0.3599 ,-0.9330},{ 0.3369,-0.9415},{ 0.3137,-0.9495},{ 0.2903,-0.9569},{ 0.2667,-0.9638},{ 0.2430 ,-0.9700},{ 0.2191,-0.9757},{ 0.1951,-0.9808},{ 0.1710,-0.9853},{ 0.1467,-0.9892},{ 0.1224 ,-0.9925},{ 0.0980,-0.9952},{ 0.0736,-0.9973},{ 0.0491,-0.9988},{ 0.0245,-0.9997},{ 0.0000 ,-1.0000},{-0.0245,-0.9997},{-0.0491,-0.9988},{-0.0736,-0.9973},{-0.0980,-0.9952},{-0.1224 ,-0.9925},{-0.1467,-0.9892},{-0.1710,-0.9853},{-0.1951,-0.9808},{-0.2191,-0.9757},{-0.2430 ,-0.9700},{-0.2667,-0.9638},{-0.2903,-0.9569},{-0.3137,-0.9495},{-0.3369,-0.9415},{-0.3599 ,-0.9330},{-0.3827,-0.9239},{-0.4052,-0.9142},{-0.4276,-0.9040},{-0.4496,-0.8932},{-0.4714 ,-0.8819},{-0.4929,-0.8701},{-0.5141,-0.8577},{-0.5350,-0.8449},{-0.5556,-0.8315},{-0.5758 ,-0.8176},{-0.5957,-0.8032},{-0.6152,-0.7883},{-0.6344,-0.7730},{-0.6532,-0.7572},{-0.6716 ,-0.7410},{-0.6895,-0.7242},{-0.7071,-0.7071},{-0.7242,-0.6895},{-0.7410,-0.6716},{-0.7572 ,-0.6532},{-0.7730,-0.6344},{-0.7883,-0.6152},{-0.8032,-0.5957},{-0.8176,-0.5758},{-0.8315 ,-0.5556},{-0.8449,-0.5350},{-0.8577,-0.5141},{-0.8701,-0.4929},{-0.8819,-0.4714},{-0.8932 ,-0.4496},{-0.9040,-0.4276},{-0.9142,-0.4052},{-0.9239,-0.3827},{-0.9330,-0.3599},{-0.9415 ,-0.3369},{-0.9495,-0.3137},{-0.9569,-0.2903},{-0.9638,-0.2667},{-0.9700,-0.2430},{-0.9757 ,-0.2191},{-0.9808,-0.1951},{-0.9853,-0.1710},{-0.9892,-0.1467},{-0.9925,-0.1224},{-0.9952 ,-0.0980},{-0.9973,-0.0736},{-0.9988,-0.0491},{-0.9997,-0.0245}};//复数乘法,进出参数均为结构体指针void JC_Mul(pJComplex pjc_1,pJComplex pjc_2,pJComplex pjc_r){pjc_r->real = (pjc_1->real)*(pjc_2->real) - (pjc_1->image)*(pjc_2->image);pjc_r->image = (pjc_1->real)*(pjc_2->image) + (pjc_1->image)*(pjc_2->real);//复数加法,进出参数均为结构体指针void JC_Add(pJComplex pjc_1,pJComplex pjc_2,pJComplex pjc_r) {pjc_r->real = (pjc_1->real) + (pjc_2->real);pjc_r->image = (pjc_1->image) + (pjc_2->image);}//复数减法,进出参数均为结构体指针void JC_Sub(pJComplex pjc_1,pJComplex pjc_2,pJComplex pjc_r) {pjc_r->real = (pjc_1->real) - (pjc_2->real);pjc_r->image = (pjc_1->image) - (pjc_2->image);}//FFT主程序void popFFT(pJComplex dataIN){uint i,j;uint mask_1,mask_2;uint tempIndex_1,tempIndex_2;uint t,m_1,m_2,n_1;JComplex buf[LENGTH/2];for(i=0;i<LENGTH;i++){mask_1 = 1;mask_2 = 1<<(LEVEL-1);tempIndex_1 = 0;for(j=0;j<LEVEL;j++){if(i&mask_1)tempIndex_1 |= mask_2;mask_1 <<= 1;mask_2 >>= 1;}dataIN[i].real = dataIN[tempIndex_1].image;dataIN[tempIndex_1].image = 0;}for(t=1;t<=LEVEL;t++){m_1 = 1<<t;//蝶形区内数据个数m_2 = m_1>>1;n_1 = LENGTH/m_1;//蝶形区个数for(i=0;i<n_1;i++){for(j=0;j<m_2;j++) //前半部索引{tempIndex_1 = (i<<t)+j; //总索引(前半部)tempIndex_2 = tempIndex_1 + m_2;//总索引(后半部)JC_Mul(&dataIN[tempIndex_2],&WN[j*n_1],&buf[j]);JC_Sub(&dataIN[tempIndex_1],&buf[j],&dataIN[tempIndex_2]); //必须先算后半部JC_Add(&dataIN[tempIndex_1],&buf[j],&dataIN[tempIndex_1]); //再算前半部 }}}return;}//主程序,测试int main(){int i;JComplex dataIN[LENGTH];for(i=0;i<LENGTH;i++)dataIN[i].image = i;popFFT(dataIN);for(i=0;i<LENGTH;i++)printf("%f\t%f\n",dataIN[i].real,dataIN[i].image);}。