基于stratix 8192点FFT算法的实现

合集下载

数字信号处理实验FFT快速傅里叶变换C语言

数字信号处理实验<第八章FFT）一、实验内容针对典型序列，用C语言编程实现基2-FFT，并与MATLAB自带的FFT函数的计算结果进行比较。

二、实验工具1.VC++6.02.matlab2018三、实验涉及知识图1按时间抽选的基2—FFT算法根据蝶形运算图，如图1，可以看出，一个点基2-FFT运算由n级蝶形运算构成，每一个蝶形结构完成下属基本迭代运算：式中m表示第m列迭代，k,j为数据所在行数。

上式的蝶形运算如图2所示。

图2四、实验设计思路首先，根据基2-FFT的算法特点，可以整理出程序设计的大致思路。

基2-FFT主要运用的就是蝶形运算来降低进行DFT运算时的运算量。

因为是基2，所以进行DFT计算的点数N必须的。

程序设计中主要解决3个问题，变址运算，复数运算，节点距离运算，及旋转因子的计算。

下面对这三个模块进行说明。

1.变址运算由蝶形图我们可以看出，FFT的输出X(k>是按正常顺序排列在存储单元中，即按X(0>,X(1>,…,X(7>的顺序排列，但是这时输入x(n>却不是按自然顺序存储的，而是按x(0>,x(4>,…,x(7>的顺序存入存储单元，所以我们要对输入的按正常顺序排列的数据进行变址存储，最后才能得到输出的有序的X(K>。

通过观察，可以发现，如果说输出数据是按原位序排列的话，那么输入数据是按倒位序排列的。

即如果输入序列的序列号用二进制数,则到位序就为。

我们需将输入的数据变为输出的倒位序存储，这里用雷德算法来实现。

下面给出雷德算法。

假如使用A[I]存的是顺序位序，而B[J]存的是倒位序。

例如 N = 8 的时候，倒位序顺序二进制表示倒位序顺序0 0 000 0004 1 100 0012 2 010 0106 3 110 0111 4 001 1005 5 101 1013 6 011 1107 7 111 111由上面的表可以看出，按自然顺序排列的二进制数，其下面一个数总是比其上面一个数大1，即下面一个数是上面一个数在最低位加1并向高位进位而得到的。

FFT实现的C语言代码（转载）-（基2FFT及IFFT算法C语言实现）

FFT实现的C语言代码（转载）-（基2FFT及IFFT算法C语言实现）FFT实现的C语言代码- -（基2FFT及IFFT算法C语言实现）Given two images A and B, use image B to cover image A. Where would we put B on A, so that the overlapping part of A and B has the most likelihood? To simplify the problem, we assume that A and B only contain numbers between 0 and 255. The difference between A and B is defined as the square sum of the differences of corresponding elements in the overlapped parts of A and B.For example, we haveA (3 * 3): a1 a2 a3B (2 * 2): b1 b2a4 a5 a6 b4 b5a7 a8 a9When B is placed on position a5, the difference of them is ((b1-a5)^2 + (b2-a6)^2 + (b4-a8)^2 + (b5-a9)^2). Now we hope to have the position of the top left corner of B that gives the minimum difference. (B must completely reside on A) It is clear that a simple solution will appear with very low efficiency when A and B have too many elements. But we can use 1-dimensional repeat convolution, which can be computed by Fast Fourier Transform (FFT), to improve the performance.A program with explanation of FFT is given below:/*** Given two sequences {a1, a2, a3.. an} and {b1, b2, b3... bn}, * their repeat convolution means:* r1 = a1*b1 + a2*b2 + a3*b3 + ... + an*bn* r2 = a1*bn + a2*b1 + a3*b2 + ... + an*bn-1* r3 = a1*bn-1 + a2*bn + a3*b1 + ... + an*bn-2* ...* rn = a1*b2 + a2*b3 + a3*b4 + ... + an-1*bn + an*b1* Notice n >= 2 and n must be power of 2.*/#include#include#include#define for if (0); else forusing namespace std;const int MaxFastBits = 16;int **gFFTBitTable = 0;int NumberOfBitsNeeded(int PowerOfTwo) {for (int i = 0;; ++i) {if (PowerOfTwo & (1 << i)) {return i;}}}int ReverseBits(int index, int NumBits) {int ret = 0;for (int i = 0; i < NumBits; ++i, index >>= 1) {ret = (ret << 1) | (index & 1);}return ret;}void InitFFT() {gFFTBitTable = new int *[MaxFastBits];for (int i = 1, length = 2; i <= MaxFastBits; ++i, length <<= 1) {gFFTBitTable[i - 1] = new int[length];for (int j = 0; j < length; ++j) {gFFTBitTable[i - 1][j] = ReverseBits(j, i);}}}inline int FastReverseBits(int i, int NumBits) {return NumBits <= MaxFastBits ? gFFTBitTable[NumBits - 1][i] : ReverseBits(i, NumBits);}void FFT(bool InverseTransform, vector >& In, vector >& Out) {if (!gFFTBitTable) { InitFFT(); }// simultaneous data copy and bit-reversal ordering into outputsint NumSamples = In.size();int NumBits = NumberOfBitsNeeded(NumSamples);for (int i = 0; i < NumSamples; ++i) {Out[FastReverseBits(i, NumBits)] = In[i];}// the FFT processdouble angle_numerator = acos(-1.) * (InverseTransform ? -2 : 2);for (int BlockEnd = 1, BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {double delta_angle = angle_numerator / BlockSize;double sin1 = sin(-delta_angle);double cos1 = cos(-delta_angle);double sin2 = sin(-delta_angle * 2);double cos2 = cos(-delta_angle * 2);for (int i = 0; i < NumSamples; i += BlockSize) {complex a1(cos1, sin1), a2(cos2, sin2);for (int j = i, n = 0; n < BlockEnd; ++j, ++n) {complex a0(2 * cos1 * a1.real() - a2.real(), 2 * cos1 * a1.imag() - a2.imag());a2 = a1;a1 = a0;complex a = a0 * Out[j + BlockEnd];Out[j + BlockEnd] = Out[j] - a;Out[j] += a;}}BlockEnd = BlockSize;}// normalize if inverse transformif (InverseTransform) {for (int i = 0; i < NumSamples; ++i) {Out[i] /= NumSamples;}}}vector convolution(vector a, vector b) {int n = a.size();vector > s(n), d1(n), d2(n), y(n);for (int i = 0; i < n; ++i) {s[i] = complex(a[i], 0);}FFT(false, s, d1);s[0] = complex(b[0], 0);for (int i = 1; i < n; ++i) {s[i] = complex(b[n - i], 0);}FFT(false, s, d2);for (int i = 0; i < n; ++i) {y[i] = d1[i] * d2[i];}FFT(true, y, s);vector ret(n);for (int i = 0; i < n; ++i) {ret[i] = s[i].real();}return ret;}int main() {double a[4] = {1, 2, 3, 4}, b[4] = {1, 2, 3, 4};vector r = convolution(vector(a, a + 4), vector(b, b + 4));// r[0] = 30 (1*1 + 2*2 + 3*3 + 4*4)// r[1] = 24 (1*4 + 2*1 + 3*2 + 4*3)// r[2] = 22 (1*3 + 2*4 + 3*1 + 4*2)// r[3] = 24 (1*2 + 2*3 + 3*4 + 4*1)return 0;}InputThe first line contains n (1 <= n <= 10), the number of test cases.For each test case, the first line contains four integers m, n, p and q, where A is a matrix of m * n, B is a matrix of p * q (2 <= m, n, p, q <= 500, m >= p, n >= q). The following m lines are the elements of A and p lines are the elements of B.OutputFor each case, print the position that gives the minimum difference (the top left corner of A is (1, 1)). You can assume that each test case has a unique solution.Sample Input22 2 2 21 23 42 31 43 3 2 20 5 50 5 50 0 05 55 5Sample Output1 11 2Author:DU, Peng。

范例FFT算法程序及分析

FFT 算法程序及分析摘要：《FFT 的算法程序分析》主要分析了按时间抽取（DIT)的快速傅立叶变换的基2FFT 算法，通过对基2FFT 算法的原理的分析及与DFT 算法运算量的比较，进一步推导出了基rFFT 算法，重点是基rFFT 算法的推导。

在具体的实例中，我们重点分析了FFT 过程中幅值大小与FFT 选用点数N 的关系,验证FFT 变换的可靠性,考察在FFT 中数据样本的长度与DFT 的点数对频谱图的影响。

关键字：基2FFT 算法,基rFFT 算法，样本长度,选用点数要求：● 学习书上第六节的内容，自己编程实现FFT 算法 . ● 给出典型信号的时域和频域图，并加以分析。

● 可尝试实现分段卷积程序。

● 论文内容含原程序、运行结果，理论分析和典型信号时域图。

一．快速傅立叶变换(FFT ）简介离散傅立叶变换（DFT ）是信号分析与处理中的一种重要的变换。

因直接计算DFT 的计算量与变换区间长度N 的平方成正比，当N 较大时，计算量太大。

所以在快速傅立叶变换(FFT ）出现以前，直接用DFT 算法进行频谱分析和信号的实时处理是不切实际的。

1965年,库利（J.W.Cooley ）和图基(J.W.Tukey ）在《计算数学》杂志上发表了“机器计算傅立叶级数的一种算法”的文章，这是一篇关于计算DFT 的一种快速有效的计算方法的文章。

它的思路建立在对DFT 运算内在规律的认识之上。

这篇文章的发表使DFT 的计算量大大减少,并导致了许多计算方法的发现。

这些算法统称为快速傅立叶变换（Fast Fourier Transform),简称FFT.1984年,法国的杜哈梅尔（P.Dohamel ）和霍尔曼（H.Hollmann ）提出的分裂基快速算法，使运算效率进一步提高。

快速傅立叶变换（FFT ）不是一种新的变换,而是离散傅立叶变换（DFT ）的一种快速算法。

FFT 分成两大类,即按时间抽取（decimation —in —time,缩写为DIT ）法和按频率抽取（decimation —in —frequency ，缩写为DIF)法。

按时间抽取的基2FFT算法分析及MATLAB实现

电子技术研发Electronics R&D电一子一技一术～按时问抽取的基2FFT 算法分析及MATLAB 实现张登奇李宏民李丹 (湖南理工学院信息与通信工程学院)摘要：DFT 是一种应用广泛的数学变换工具，MATLAB 是～款功能强大的科学计算语言。

MATLAB 提供的FFT 函数解决了DFT 的快速计算问题，但由于它是内建函数而不能了解到软件实现的过程。

文章以按时间抽取的基2FFT 算法为例，根据快速傅里叶变换的原理和规律，绘出了算法实现的程序框图，列出TIVIAT L&B 环境下软件实现的程序，建立了从算法理论到程序实现的完整概念。

关键词：数字信号处理；离散傅里叶变换；快速傅里叶变换：按时间抽取；MATL Pd3Analysis of Decimation ．in ．time R adi x 2 FF T Algorithmand MATLAB Based Realizatio n Zh a ng D e ng q i Li Ho ngr nin Li Da n(Colleg e of Informat ion&C ommu nica tion Engineering ，Hunan Institute ofS ci en ce a nd Technology) Abs tra ct ：DFT i S a wi dely used mathem atica I transfol'nl tool 。

and MATLAB iS a powe rf ul scientific c om p u t in gl an g u ag e ． The FFT functionof MATLAB c a n solve the problems of fast c alc ul at io n of D FT ,bu t it can’t understand the pr o c e s s o f s o f t w a r e realization b ecause i t i s a built-in f un c t io n ．T a k in g the algori thm of the ra di x 一2DI T F F T as the exam ple andbased o n the p ri nc ip le a nd rules of FF T,this paper draws the bl ock diagrams of the alg o r i t hm realization ，lists the p r o c e d u r es of the s of t w ar erealization u nd e r MATLAB environment ．and also est ab li sh es the complete concept from al g o ri t h m theories to p ro gr am realization ．Ke y w o rd s ：di g it a l signal processing(DSP)；di screte Fou r ie r transform(DFT)；fastF ou r i er transform(FFT)；decimation in time(DIT)：Ⅳ队TLAB变换工具，自1965年库利和图基提出基2FFT 算法[1]以来，主：珲嚣裟嚣篙∞0引言离散傅里变换(DFT)是数字信号处理中重要的数学图1 DIT 蝶形运算流图符号现已有多种快速算法，且还在不断地进行研究和探索。

基于STM32芯片的128点FFT

Main函数/******************** (C) COPYRIGHT 2008 STMicroelectronics ******************** * File Name : main.c* Author : MCD Application Team* Version : V2.0.1* Date : 06/13/2008* Description : Main program body******************************************************************************* ** THE PRESENT FIRMWARE WHICH IS FOR GUIDANCE ONL Y AIMS A T PROVIDING CUSTOMERS* WITH CODING INFORMA TION REGARDING THEIR PRODUCTS IN ORDER FOR THEM TO SA VE TIME.* AS A RESULT, STMICROELECTRONICS SHALL NOT BE HELD LIABLE FOR ANY DIRECT,* INDIRECT OR CONSEQUENTIAL DAMAGES WITH RESPECT TO ANY CLAIMS ARISING FROM THE* CONTENT OF SUCH FIRMW ARE AND/OR THE USE MADE BY CUSTOMERS OF THE CODING* INFORMATION CONTAINED HEREIN IN CONNECTION WITH THEIR PRODUCTS.******************************************************************************* //* Includes ------------------------------------------------------------------*/#include "stm32f10x_lib.h"#include "stdio.h"#include "math.h"#include "stm32_dsp.h"/* Private typedef -----------------------------------------------------------*//* Private define ------------------------------------------------------------*/#define ADC1_DR_Address ((u32)0x4001244C)#define USART1_DR_Base ((u32)0x40013804)#define LEN 128/* Private macro -------------------------------------------------------------*//* Private variables ---------------------------------------------------------*/ErrorStatus HSEStartUpStatus;/* Private function prototypes -----------------------------------------------*/void RCC_Configuration(void);void GPIO_Configuration(void);void NVIC_Configuration(void);void TIM_BaseConfiguration(void);void DMA_Configuration(void);void USART_Configuration(void);void ADC_Configuration(void);void ADC_ALL_Init(void);void FFT_Configuration(double *pr,double *pi,int n,int k,double*fr,double*fi,int l,int il);void FFT_solve(void);void Delay(u32 counter);/* Private functions ---------------------------------------------------------*/u16 data_buff[LEN];u16 px_buff[LEN];u16 pz_buff[LEN];double x[128];double pr[128];double pi[128];double fr[128];double fi[128];double mo[128];//int *px = (int*)0x1f00;//int *pz = (int*)0x1f80;int xm,zm,i,t = 0;u8 Flag = 0;/****************************************************************************** ** Function Name : main* Description : Main program* Input : None* Output : None* Return : None******************************************************************************* /int main(void){#ifdef DEBUGdebug();#endif/* System clocks configuration ---------------------------------------------*/RCC_Configuration();/* NVIC configuration ------------------------------------------------------*/NVIC_Configuration();/* GPIO configuration ------------------------------------------------------*/GPIO_Configuration();/* ADC1 configuration ------------------------------------------------------*/ADC_Configuration();/*USART Configuration--------------------------------------------------------*/USART_Configuration();ADC_ALL_Init();while (1){FFT_solve();}}/****************************************************************************** ** Function Name : RCC_Configuration* Description : Configures the different system clocks.* Input : None* Output : None* Return : None******************************************************************************* /void RCC_Configuration(void){/* RCC system reset(for debug purpose) */RCC_DeInit();/* Enable HSE */RCC_HSEConfig(RCC_HSE_ON);/* Wait till HSE is ready */HSEStartUpStatus = RCC_WaitForHSEStartUp();if(HSEStartUpStatus == SUCCESS){/* Enable Prefetch Buffer */FLASH_PrefetchBufferCmd(FLASH_PrefetchBuffer_Enable);/* Flash 2 wait state */FLASH_SetLatency(FLASH_Latency_2);/* HCLK = SYSCLK */RCC_HCLKConfig(RCC_SYSCLK_Div1);/* PCLK2 = HCLK */RCC_PCLK2Config(RCC_HCLK_Div1);/* PCLK1 = HCLK/2 */RCC_PCLK1Config(RCC_HCLK_Div2);/* ADCCLK = PCLK2/4 */RCC_ADCCLKConfig(RCC_PCLK2_Div4);/* PLLCLK = 8MHz * 7 = 56 MHz */RCC_PLLConfig(RCC_PLLSource_HSE_Div1, RCC_PLLMul_7);/* Enable PLL */RCC_PLLCmd(ENABLE);/* Wait till PLL is ready */while(RCC_GetFlagStatus(RCC_FLAG_PLLRDY) == RESET){}/* Select PLL as system clock source */RCC_SYSCLKConfig(RCC_SYSCLKSource_PLLCLK);/* Wait till PLL is used as system clock source */while(RCC_GetSYSCLKSource() != 0x08){}}/* Enable peripheral clocks --------------------------------------------------*//* Enable DMA1 clock */RCC_AHBPeriphClockCmd(RCC_AHBPeriph_DMA1, ENABLE);/* Enable USART2 clock */RCC_APB1PeriphClockCmd(RCC_APB2Periph_USART1, ENABLE);/*Enable TIM2 clock*/RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE);/* Enable ADC1 and GPIOC clock */RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1 | RCC_APB2Periph_GPIOC | RCC_APB2Periph_GPIOA, ENABLE);}/****************************************************************************** ** Function Name : GPIO_Configuration* Description : Configures the different GPIO ports.* Input : None* Output : None* Return : None******************************************************************************* /void GPIO_Configuration(void){GPIO_InitTypeDef GPIO_InitStructure;/* Configure PC.04 (ADC Channel14) as analog input -------------------------*/GPIO_InitStructure.GPIO_Pin = GPIO_Pin_4;GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AIN;GPIO_Init(GPIOC, &GPIO_InitStructure);/* Configure PA.03 (ADC Channel14) as analog input -------------------------*/GPIO_InitStructure.GPIO_Pin = GPIO_Pin_3;GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AIN;GPIO_Init(GPIOA, &GPIO_InitStructure);/* Configure USART1 Rx (PA.10) as input floating */GPIO_InitStructure.GPIO_Pin = GPIO_Pin_10;GPIO_InitStructure.GPIO_Mode = GPIO_Mode_IN_FLOATING;GPIO_Init(GPIOA, &GPIO_InitStructure);/* Configure USART1 Tx (PA.09) as alternate function push-pull */GPIO_InitStructure.GPIO_Pin = GPIO_Pin_9;GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AF_PP;GPIO_Init(GPIOA, &GPIO_InitStructure);}/****************************************************************************** ** Function Name : NVIC_Configuration* Description : Configures Vector Table base location.* Input : None* Output : None* Return : None*******************************************************************************/void NVIC_Configuration(void){NVIC_InitTypeDef NVIC_InitStructure;#ifdef VECT_TAB_RAM/* Set the Vector Table base location at 0x20000000 */NVIC_SetVectorTable(NVIC_VectTab_RAM, 0x0);#else /* VECT_TAB_FLASH *//* Set the Vector Table base location at 0x08000000 */NVIC_SetVectorTable(NVIC_VectTab_FLASH, 0x0);#endif/* Enable the TIM2 gloabal Interrupt */NVIC_InitStructure.NVIC_IRQChannel = TIM2_IRQChannel;NVIC_InitStructure.NVIC_IRQChannelPreemptionPriority = 0;NVIC_InitStructure.NVIC_IRQChannelSubPriority = 2;NVIC_InitStructure.NVIC_IRQChannelCmd = ENABLE;NVIC_Init(&NVIC_InitStructure);/* Enable the ADC1_EOC Interrupt */NVIC_InitStructure.NVIC_IRQChannel = ADC1_2_IRQChannel;NVIC_InitStructure.NVIC_IRQChannelPreemptionPriority = 0;NVIC_InitStructure.NVIC_IRQChannelSubPriority = 1;NVIC_InitStructure.NVIC_IRQChannelCmd = ENABLE;NVIC_Init(&NVIC_InitStructure);}/****************************************************************************** ** Function Name : TIM_BaseConfiguration* Description : Configures Vector Table base location.* Input : None* Output : None* Return : None******************************************************************************* /void TIM_BaseConfiguration(void){TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure;/* Time base configuration */TIM_TimeBaseStructure.TIM_Period = 65535;TIM_TimeBaseStructure.TIM_Prescaler = 0;TIM_TimeBaseStructure.TIM_ClockDivision = 0;TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up;TIM_TimeBaseInit(TIM2, &TIM_TimeBaseStructure);/*Enable TIM2*/TIM_Cmd(TIM2,ENABLE);/* TIM IT enable */TIM_ITConfig(TIM2, TIM_IT_Update , ENABLE);}/****************************************************************************** ** Function Name : ADC_Configuration* Description : Configures Vector Table base location.* Input : None* Output : None* Return : None******************************************************************************* /void ADC_Configuration(void){ADC_InitTypeDef ADC_InitStructure;ADC_InitStructure.ADC_Mode = ADC_Mode_Independent; /*工作在独立模式*/ADC_InitStructure.ADC_ScanConvMode = ENABLE; /*规定了模数转换工作在扫描模式（多通道）还是单次（单通道）模式*/ADC_InitStructure.ADC_ContinuousConvMode = ENABLE; /*规定了模数转换工作在连续还是单次模式*/ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_None; /*转换由软件而不是外部触发启动*/ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right; /*规定了ADC 数据向左边对齐还是向右边对齐*/ADC_InitStructure.ADC_NbrOfChannel = 1; /*规定了顺序进行规则转换的ADC通道的数目*/ADC_Init(ADC1, &ADC_InitStructure);}void ADC_ALL_Init(void){/* ADC1 regular channel14 configuration */ADC_RegularChannelConfig(ADC1, ADC_Channel_14, 1, ADC_SampleTime_7Cycles5);/*设置指定ADC的规则组通道，设置它们的转化顺序和采样时间*//* Enable ADC1 EOC interrupts */ADC_ITConfig(ADC1, ADC_IT_EOC , ENABLE);/* Enable ADC1 */ADC_Cmd(ADC1, ENABLE);/* Enable ADC1 reset calibaration register */ADC_ResetCalibration(ADC1); /*重置指定的ADC的校准寄存器*//* Check the end of ADC1 reset calibration register */while(ADC_GetResetCalibrationStatus(ADC1));/*获取ADC重置校准寄存器的状态*//* Start ADC1 calibaration */ADC_StartCalibration(ADC1);/* Check the end of ADC1 calibration */while(ADC_GetCalibrationStatus(ADC1));/* Start ADC1 Software Conversion */ADC_SoftwareStartConvCmd(ADC1, ENABLE);/*使能指定的ADC的软件转换启动功能*/}/****************************************************************************** ** Function Name : USART_Configuration* Description : Configures Vector Table base location.* Input : None* Output : None* Return : None******************************************************************************* /void USART_Configuration(void){USART_InitTypeDef USART_InitStructure;USART_ART_BaudRate = 9600;USART_ART_WordLength = USART_WordLength_8b;USART_ART_StopBits = USART_StopBits_1;USART_ART_Parity = USART_Parity_No;USART_ART_HardwareFlowControl = USART_HardwareFlowControl_None;USART_ART_Mode = USART_Mode_Rx | USART_Mode_Tx;/* Configure USART1 */USART_Init(USART1, &USART_InitStructure);/*Enable USART1 */USART_Cmd(USART1, ENABLE);}/******************************************************************************* Function Name : FFT_Configuration* Description :* Input : None* Output : None* Return : None******************************************************************************* /void FFT_Configuration(double *pr,double *pi,int n,int k,double*fr,double*fi,int l,int il)//int n,k,l,il; double pr[],pi[],fr[],fi[];/*pr(实部)，pi(虚部)，n(点数)，k(阶数)*/{int it,m,is,i,j,nv,l0;double p,q,s,vr,vi,poddr,poddi;/*雷德(Rader)算法，实现数据的倒位序排列-----------*/for (it=0; it<=n-1; it++){m=it;is=0;for (i=0; i<=k-1; i++){j=m/2;is=2*is+(m-2*j);m=j;}fr[it]=pr[is];fi[it]=pi[is];}pr[0]=1.0;pi[0]=0.0;p=6.283185306/(1.0*n);pr[1]=cos(p);pi[1]=-sin(p);if (l!=0)pi[1]=-pi[1];for (i=2; i<=n-1; i++){p=pr[i-1]*pr[1];q=pi[i-1]*pi[1];s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);pr[i]=p-q;pi[i]=s-p-q;}for (it=0; it<=n-2; it=it+2){vr=fr[it];vi=fi[it];fr[it]=vr+fr[it+1];fi[it]=vi+fi[it+1];fr[it+1]=vr-fr[it+1];fi[it+1]=vi-fi[it+1];}m=n/2;nv=2;for (l0=k-2; l0>=0; l0--){m=m/2;nv=2*nv;for (it=0; it<=(m-1)*nv; it=it+nv)for (j=0; j<=(nv/2)-1; j++){p=pr[m*j]*fr[it+j+nv/2];q=pi[m*j]*fi[it+j+nv/2];s=pr[m*j]+pi[m*j];s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);poddr=p-q;poddi=s-p-q;fr[it+j+nv/2]=fr[it+j]-poddr;fi[it+j+nv/2]=fi[it+j]-poddi;fr[it+j]=fr[it+j]+poddr;fi[it+j]=fi[it+j]+poddi;}}if (l!=0)for (i=0; i<=n-1; i++){ fr[i]=fr[i]/(1.0*n);fi[i]=fi[i]/(1.0*n);}if (il!=0)for (i=0; i<=n-1; i++){pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);if (fabs(fr[i])<0.000001*fabs(fi[i])){if ((fi[i]*fr[i])>0) pi[i]=90.0;elsepi[i]=-90.0;}elsepi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;}}/****************************************************************************** ** Function Name : FFT_solve* Description :* Input : None* Output : None* Return : None******************************************************************************* /void FFT_solve(void){if(Flag){/*寄存器接收满后Flag = 1,关闭ADC中断,进行FFT运算--------*/ADC_Cmd(ADC1, DISABLE);//px = (int*)0x1f00;/*转换存储位置-------------*/for(i=0;i<128;i++){//*px=data_buff[i];//px++;px_buff[i] = data_buff[i];}//px = (int*)0x1f00;/*将ADC采来的数据转换为实际的幅度值--------------*/for (i=0; i<=128; i++){//xm=*px;xm = px_buff[i];x[i]=xm/32768.0;pr[i]=x[i];/*pr数组存放函数点的实部-------------*/pi[i]=0; /*pi数组存放函数点的虚部-------------*///px++;}/*进行FFT运算----------------------*/FFT_Configuration(pr,pi,128,7,fr,fi,0,1);//pz = (int*)0x1f80;/*将最终计算出的值放入pz 中----------------------*/for (i=0;i<=128;i++){mo[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i]);zm = (int)(mo[i]*1000.0);//*pz = zm;//pz++;pz_buff[i] = zm;}}Flag = 0;}/****************************************************************************** ** Function Name : Delay* Description :* Input : None* Output : None* Return : None******************************************************************************* /void Delay(u32 counter){for(;counter > 0;counter--);}#ifdef DEBUG/****************************************************************************** ** Function Name : assert_failed* Description : Reports the name of the source file and the source line number* where the assert_param error has occurred.* Input : - file: pointer to the source file name* - line: assert_param error line source number* Output : None* Return : None******************************************************************************* /void assert_failed(u8* file, u32 line){/* User can add his own implementation to report the file name and line number, ex: printf("Wrong parameters value: file %s on line %d\r\n", file, line) *//* Infinite loop */while (1){}}#endif/******************* (C) COPYRIGHT 2008 STMicroelectronics *****END OF FILE****/stm32f10x_it.c/******************** (C) COPYRIGHT 2008 STMicroelectronics ******************** * File Name : stm32f10x_it.c* Author : MCD Application Team* Version : V2.0.1* Date : 06/13/2008* Description : Main Interrupt Service Routines.* This file provides template for all exceptions handler* and peripherals interrupt service routine.******************************************************************************* ** THE PRESENT FIRMWARE WHICH IS FOR GUIDANCE ONL Y AIMS A T PROVIDING CUSTOMERS* WITH CODING INFORMA TION REGARDING THEIR PRODUCTS IN ORDER FOR THEM TO SA VE TIME.* AS A RESULT, STMICROELECTRONICS SHALL NOT BE HELD LIABLE FOR ANY DIRECT,* INDIRECT OR CONSEQUENTIAL DAMAGES WITH RESPECT TO ANY CLAIMS ARISING FROM THE* CONTENT OF SUCH FIRMW ARE AND/OR THE USE MADE BY CUSTOMERS OF THE CODING* INFORMATION CONTAINED HEREIN IN CONNECTION WITH THEIR PRODUCTS.******************************************************************************* //* Includes ------------------------------------------------------------------*/#include "stm32f10x_it.h"#include "stdio.h"/* Private typedef -----------------------------------------------------------*//* Private define ------------------------------------------------------------*/#define ADC1_DR_Address ((u32)0x4001244C)#define LEN 128/* Private macro -------------------------------------------------------------*//* Private variables ---------------------------------------------------------*//* Private function prototypes -----------------------------------------------*//* Private functions ---------------------------------------------------------*/extern u16 data_buff[LEN];extern int t ;extern u8 Flag;u32 *P = (u32*)0x4001244C;/****************************************************************************** ** Function Name : NMIException* Description : This function handles NMI exception.* Input : None* Output : None* Return : None******************************************************************************* /void NMIException(void){}/****************************************************************************** ** Function Name : HardFaultException* Description : This function handles Hard Fault exception.* Input : None* Output : None* Return : None******************************************************************************* /void HardFaultException(void){/* Go to infinite loop when Hard Fault exception occurs */while (1){}}/****************************************************************************** ** Function Name : MemManageException* Description : This function handles Memory Manage exception.* Input : None* Output : None* Return : None******************************************************************************* /void MemManageException(void){/* Go to infinite loop when Memory Manage exception occurs */while (1){}}/****************************************************************************** ** Function Name : BusFaultException* Description : This function handles Bus Fault exception.* Input : None* Output : None* Return : None******************************************************************************* /void BusFaultException(void){/* Go to infinite loop when Bus Fault exception occurs */while (1){}}/****************************************************************************** ** Function Name : UsageFaultException* Description : This function handles Usage Fault exception.* Input : None* Output : None* Return : None******************************************************************************* /void UsageFaultException(void){/* Go to infinite loop when Usage Fault exception occurs */while (1){}}/****************************************************************************** ** Function Name : DebugMonitor* Description : This function handles Debug Monitor exception.* Input : None* Output : None* Return : None******************************************************************************* /void DebugMonitor(void){}/****************************************************************************** ** Function Name : SVCHandler* Description : This function handles SVCall exception.* Input : None* Output : None* Return : None******************************************************************************* /void SVCHandler(void){}/****************************************************************************** ** Function Name : PendSVC* Description : This function handles PendSVC exception.* Input : None* Output : None* Return : None******************************************************************************* /void PendSVC(void){}/****************************************************************************** ** Function Name : SysTickHandler* Description : This function handles SysTick Handler.* Input : None* Output : None* Return : None******************************************************************************* /void SysTickHandler(void){}/****************************************************************************** ** Function Name : WWDG_IRQHandler* Description : This function handles WWDG interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void WWDG_IRQHandler(void){}/****************************************************************************** ** Function Name : PVD_IRQHandler* Description : This function handles PVD interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void PVD_IRQHandler(void){}/****************************************************************************** ** Function Name : TAMPER_IRQHandler* Description : This function handles Tamper interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void TAMPER_IRQHandler(void){}/****************************************************************************** ** Function Name : RTC_IRQHandler* Description : This function handles RTC global interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void RTC_IRQHandler(void){}/****************************************************************************** ** Function Name : FLASH_IRQHandler* Description : This function handles Flash interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void FLASH_IRQHandler(void){}/****************************************************************************** ** Function Name : RCC_IRQHandler* Description : This function handles RCC interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void RCC_IRQHandler(void){}/****************************************************************************** ** Function Name : EXTI0_IRQHandler* Description : This function handles External interrupt Line 0 request.* Input : None* Output : None* Return : None******************************************************************************* /void EXTI0_IRQHandler(void){}/****************************************************************************** ** Function Name : EXTI1_IRQHandler* Description : This function handles External interrupt Line 1 request.* Input : None* Output : None* Return : None******************************************************************************* /void EXTI1_IRQHandler(void){}/****************************************************************************** ** Function Name : EXTI2_IRQHandler* Description : This function handles External interrupt Line 2 request.* Input : None* Output : None* Return : None******************************************************************************* /void EXTI2_IRQHandler(void){}/****************************************************************************** ** Function Name : EXTI3_IRQHandler* Description : This function handles External interrupt Line 3 request.* Input : None* Output : None* Return : None******************************************************************************* /void EXTI3_IRQHandler(void){}/****************************************************************************** ** Function Name : EXTI4_IRQHandler* Description : This function handles External interrupt Line 4 request.* Input : None* Output : None* Return : None******************************************************************************* /void EXTI4_IRQHandler(void){}/****************************************************************************** ** Function Name : DMA1_Channel1_IRQHandler* Description : This function handles DMA1 Channel 1 interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void DMA1_Channel1_IRQHandler(void){}/****************************************************************************** ** Function Name : DMA1_Channel2_IRQHandler* Description : This function handles DMA1 Channel 2 interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void DMA1_Channel2_IRQHandler(void){}/******************************************************************************** Function Name : DMA1_Channel3_IRQHandler* Description : This function handles DMA1 Channel 3 interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void DMA1_Channel3_IRQHandler(void){}/****************************************************************************** ** Function Name : DMA1_Channel4_IRQHandler* Description : This function handles DMA1 Channel 4 interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void DMA1_Channel4_IRQHandler(void){}/****************************************************************************** ** Function Name : DMA1_Channel5_IRQHandler* Description : This function handles DMA1 Channel 5 interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void DMA1_Channel5_IRQHandler(void){}/****************************************************************************** ** Function Name : DMA1_Channel6_IRQHandler* Description : This function handles DMA1 Channel 6 interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void DMA1_Channel6_IRQHandler(void){}/****************************************************************************** ** Function Name : DMA1_Channel7_IRQHandler* Description : This function handles DMA1 Channel 7 interrupt request.* Input : None* Output : None* Return : None******************************************************************************* /void DMA1_Channel7_IRQHandler(void){}/****************************************************************************** ** Function Name : ADC1_2_IRQHandler* Description : This function handles ADC1 and ADC2 global interrupts requests.* Input : None* Output : None* Return : None******************************************************************************* /void ADC1_2_IRQHandler(void){if(ADC_GetITStatus(ADC1, ADC_IT_EOC) == SET){ADC_ClearITPendingBit(ADC1, ADC_IT_EOC);//*P = ADC1_DR_Address;data_buff[t] = *P;t++;while(t == 128){Flag = 1;t = 0;}}}。

FFT的C语言算法实现

FFT的C语言算法实现程序如下：/************FFT***********/#include <stdio.h>#include <math.h>#include <stdlib.h>#define N 1000typedef struct{double real;double img;}complex;void fft(); /*快速傅里叶变换*/void ifft(); /*快速傅里叶逆变换*/void initW();void change();void add(complex ,complex ,complex *); /*复数加法*/void mul(complex ,complex ,complex *); /*复数乘法*/void sub(complex ,complex ,complex *); /*复数减法*/void divi(complex ,complex ,complex *);/*复数除法*/void output(); /*输出结果*/complex x[N], *W;/*输出序列的值*/int size_x=0;/*输入序列的长度，只限2的N次方*/ double PI;int main(){int i,method;system("cls");PI=atan(1)*4;printf("Please input the size of x:\n");/*输入序列的长度*/scanf("%d",&size_x);printf("Please input the data in x[N]:(such as:5 6)\n");/*输入序列对应的值*/for(i=0;i<size_x;i++)scanf("%lf %lf",&x[i].real,&x[i].img);initW();/*选择FFT或逆FFT运算*/printf("Use FFT(0) or IFFT(1)?\n");scanf("%d",&method);if(method==0)fft();elseifft();output();return 0;}/*进行基-2 FFT运算*/void fft(){int i=0,j=0,k=0,l=0;complex up,down,product;change();for(i=0;i< log(size_x)/log(2) ;i++) {l=1<<i;for(j=0;j<size_x;j+= 2*l ){for(k=0;k<l;k++){mul(x[j+k+l],W[size_x*k/2/l],&product); add(x[j+k],product,&up);sub(x[j+k],product,&down);x[j+k]=up;x[j+k+l]=down;}}}}void ifft(){int i=0,j=0,k=0,l=size_x;complex up,down;for(i=0;i< (int)( log(size_x)/log(2) );i++) /*蝶形运算*/{l/=2;for(j=0;j<size_x;j+= 2*l ){for(k=0;k<l;k++){add(x[j+k],x[j+k+l],&up);up.real/=2;up.img/=2;sub(x[j+k],x[j+k+l],&down);down.real/=2;down.img/=2;divi(down,W[size_x*k/2/l],&down);x[j+k]=up;x[j+k+l]=down;}}}change();}void initW(){int i;W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;i<size_x;i++){W[i].real=cos(2*PI/size_x*i);W[i].img=-1*sin(2*PI/size_x*i);}}void change(){complex temp;unsigned short i=0,j=0,k=0;double t;for(i=0;i<size_x;i++){k=i;j=0;t=(log(size_x)/log(2));while( (t--)>0 ){j=j<<1;j|=(k & 1);k=k>>1;}if(j>i){temp=x[i];x[i]=x[j];x[j]=temp;}}}void output() /*输出结果*/{int i;printf("The result are as follows\n"); for(i=0;i<size_x;i++){printf("%.4f",x[i].real);if(x[i].img>=0.0001)printf("+%.4fj\n",x[i].img);else if(fabs(x[i].img)<0.0001)printf("\n");elseprintf("%.4fj\n",x[i].img);}}void add(complex a,complex b,complex *c) {c->real=a.real+b.real;c->img=a.img+b.img;}void mul(complex a,complex b,complex *c) {c->real=a.real*b.real - a.img*b.img;c->img=a.real*b.img + a.img*b.real;}void sub(complex a,complex b,complex *c) {c->real=a.real-b.real;c->img=a.img-b.img;}void divi(complex a,complex b,complex *c) {c->real=( a.real*b.real+a.img*b.img )/(b.real*b.real+b.img*b.img);c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.i mg);}。

TI库实现2048点或更大点数FFT

在DSP运算中，经常需要把输入时域信号在频域进行处理之后，再还原为时域信号，这样就需要进行FFT 和IFFT运算：x(n) -> FFT -> X(f) -> 频域处理-> Y(f) -> IFFT -> y(n)而一般的DSP芯片只支持整数运算，也就是说只能进行定点小数计算。

N点FFT计算出0… N-1，N个复数：0，A，N/2，A*，A为(N/2-1)个复数，A*为A的共轭复数。

FFT的公式为：NX(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N)、1 < = k < = N.n = 1IFFT的公式为：Nx(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N)、1 < = n < = N.k = 1假设我们对ADC转换器转换的数字信号进行FFT运算，若输入数据为16bit的短整型数，我们可以把它看作Q15的从-1到1之间的小数。

根据FFT的公式我们可以知道，FFT变换之后的结果将超出这个范围。

例如在matlab中输入fft(sin([1:8]*0.5))，可以看到结果：2.8597，-0.8019 -3.0216i，0.4312 - 0.8301i，0.5638 - 0.3251i，0.5895，0.5638 + 0.3251i，0.4312 + 0.8301i，-0.8019 + 3.0216i实际上，FFT变换之后的数据的范围在-N到N之间，N为FFT的点数。

为了正确地表达-N到N之间的数值，输出数据的Q值将变小，例如若N=1024，输入数据为Q15的话，那么输出数据则必须为Q5才能够确保结果不会溢出。

这样的结果将丢失很多信息，以至于IFFT无法还原为原来的数据。

如下图所示：Q15 -> 1024FFT -> Q5 -> 1024IFFT -> Q5这样，经过FFT和IFFT变换之后，数据从Q15变成了Q5，丢失了10bit的信息。

快速傅里叶变换FFT原理与实现

FFT原理与实现2010-10-07 21:10:09| 分类：数字信号处理 | 标签：fft dft |举报|字号订阅在数字信号处理中常常需要用到离散傅立叶变换(DFT)，以获取信号的频域特征。

尽管传统的DFT算法能够获取信号频域特征，但是算法计算量大，耗时长，不利于计算机实时对信号进行处理。

因此至DFT被发现以来，在很长的一段时间内都不能被应用到实际的工程项目中，直到一种快速的离散傅立叶计算方法——FFT，被发现，离散傅立叶变换才在实际的工程中得到广泛应用。

需要强调的是，FFT并不是一种新的频域特征获取方式，而是DFT的一种快速实现算法。

本文就FFT的原理以及具体实现过程进行详尽讲解。

DFT计算公式本文不加推导地直接给出DFT的计算公式：其中x(n)表示输入的离散数字信号序列，WN为旋转因子，X(k)为输入序列x(n)对应的N个离散频率点的相对幅度。

一般情况下，假设x(n)来自于低通采样，采样频率为fs，那么X(k)表示了从-fs/2率开始，频率间隔为fs/N,到fs/2-fs/N截至的N个频率点的相对幅度。

因为DFT计算得到的一组离散频率幅度值实际上是在频率轴上从成周期变化的，即X(k+N)=X(k)。

因此任意取连续的N个点均可以表示DFT的计算效果，负频率成分比较抽象，难于理解，根据X(k)的周期特性，于是我们又可以认为X(k)表示了从零频率开始，频率间隔为fs/N,到fs-fs/N 截至的N个频率点的相对幅度。

N点DFT的计算量根据(1)式给出的DFT计算公式，我们可以知道每计算一个频率点X(k)均需要进行N次复数乘法和N-1次复数加法，计算N各点的X(k)共需要N^2次复数乘法和N*(N-1)次复数加法。

当x(n)为实数的情况下，计算N点的DFT需要2*N^2次实数乘法，2*N*(N-1)次实数加法。

旋转因子WN的特性1.WN的对称性2.WN的周期性3.WN的可约性根据以上这些性质，我们可以得到式(5)的一系列有用结果基-2 FFT算法推导假设采样序列点数为N=2^L，L为整数，如果不满足这个条件可以人为地添加若干个0以使采样序列点数满足这一要求。

用 VC++6.0实现离散傅立叶（DFT）变换

用 VC++6.0实现离散傅立叶（DFT）变换马晓敏【摘要】The Discrete Fourier Transform (DFT), which can convert signal-In-Time to signal-In-Frequency, is a very important method in signal analysis. This paper mainly introduces two algorithms realizing the DFT, i.e. Decimation-In-Time and Decimation-In-Frequency. The comparison between them is also made.%离散傅立叶变换是信号分析中的一种重要手段，它将时域信号转换成频域信号。

本文主要介绍离散傅立叶变换的两种算法，即时域抽取法，频域抽取法，并对两种算法进行比较。

【期刊名称】《安徽电子信息职业技术学院学报》【年(卷),期】2014(000)003【总页数】4页(P5-7,31)【关键词】离散傅立叶变换;信号分析;时域抽取法;频域抽取法【作者】马晓敏【作者单位】中国电子科技集团公司第 41 研究所信息中心，安徽蚌埠 233006【正文语种】中文【中图分类】TN919.1在电子工程上，最为熟悉的表示信号的方法是时域表示法，但人们在很多情况下采用频域表示法来描述信号与系统，因此，常常需要将时域信号转换成频域信号来分析。

实现这一变换的方法很多，最为常用的是傅立叶变换。

很显然，要通过计算机来实现这一变换就需要先通过数据采样将原始连续信号转换成离散信号，并将计算范围收缩到一个有限区间。

因此，在一定的误差范围内，我们将使用离散傅立叶变换将时域信号转换成频域信号。

要计算一个N点的离散傅立叶变换需要同一个N×N点的W矩阵（关于W矩阵请参阅信号与系统方面或数学方面的书籍）相运算，随着N值的增大，运算次数显著上升，当点数达到1024时，需要进行复数乘法运算1048576次。

LTE小区搜索算法研究王延双

分类号 TN929.5 密级重庆邮电大学硕士学位论文论文题目LTE-TDD系统中小区搜索算法研究及FPGA实现英文题目 The Research of Cell Search Algorithms in LTETDD Systems and FPGA implementation 硕士研究生王延双指导教师陈发堂副教授学科专业通信与信息系统论文提交日期 2011.4 论文答辩日期2011.5论文评阅人杨浩副教授重庆大学柴蓉副教授重庆邮电大学答辩委员会主席陈前斌教授重庆邮电大学2011年 5 月 28 日摘要LTE 作为3G 技术的演进，在频带利用率，带宽的可配置性和数据传输速率上都有很好的优越性，从而受到全球的关注。

而小区搜索是指移动台UE 在初始接入小区时或移动台UE 在进行小区切换过程时，找到服务小区ID 号以及与服务小区取得时间和频率的同步过程。

作为UE 读取基站信息的前提，小区搜索过程的好坏直接影响整个通信系统的正常通信，只有完成成功的小区搜索，才能保证移动台UE 与基站之间实现正常的数据交换。

本文首先介绍OFDM 和MIMO 技术在LTE 系统中的应用，接着介绍与LTE-TDD 系统小区搜索算法和搜索过程机制相关的物理层知识。

然后重点对小区搜索进行了方案设计与具体分析，在此基础上提出了，UE 在初始接入小区时的符号同步算法，组内小区ID 号检测算法和频率同步算法：在符号同步过程中，采用分级结构进行检测，首先对本地预先生成的三组主同步信号进行128点IFFT 运算，然后用基于互相关的同步算法，分别与经过降采样的接收数据进行128点滑动相关，然后从三组相关结果中找出最大相关值，由最大值来确定组内小区ID 号和粗定时同步的位置，完成粗符号同步和组内小区ID 号检测，接着通过基于CP 的最大似然检测算法来进行符号精同步操作，从而完成小区搜索算法中的整个符号同步过程。

(2)ID N 在频率同步算法分析中，首先利用基于CP 的最大似然检测算法来进行小数倍频偏估计，然后利用主同步信号良好的相关特性，将接收端的信号经FFT 变换后与本地的主同步信号进行频域相关运算，如果不存在整数倍频偏时，本地的主同步信号序列与接收信号做相关后，峰值应该出现在相关后序列的正中间的位置；而当存在n 个整数倍频偏时，则做完相关运算后，序列峰值应该出现在距离正中间为n 个点的位置上，即偏离了中间位置n 个点。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

3 仿真与验证
本设计采用自顶向下的设计思路 ,完成系统模块划分和各个功能模块的 Verilog代码编写。选用 Quar2 tusⅡ作为软件开发平台 ,利用该软件宏模块的调用功能产生 RAM、ROM、乘法器等模块 ,以节省设计时间 , 提高系统工作效率。得出整个系统占用 30% FPGA 的逻辑单元和 90%片上 RAM , Clock最快为 75 MHz。系统在工作主频率 75 MHz情况下通过板级验证 ,完成全部 1 024点复数 FFT运算时间大约需要 250μs。
关键词 : FFT;频率抽取 ;蝶形运算 ; FPGA 中图分类号 : TN402
0 引言
随着数字技术的飞速发展 , DSP (数字信号处理器 )已广泛应用于通信、多媒体、医疗仪器和军事等领域。 FFT (快速傅里叶变换 )是 DSP的核心技术之一和 DFT (离散傅里叶变换 )的快速算法 ,作为时域和频域转换的基本运算 ,是数字谱分析的必要前提 ,在雷达、观测、跟踪、高速图像处理、保密无线通信和数字通信等领域有广泛应用。而 FPGA (现场可编程门阵列 )内部含有大量逻辑单元和高速 RAM 模块 ,使 FFT算法可以灵活快速地实现 ,并具有很高的性能。
本系统选用 GE02 开发板 (处理器为东南大学研 ·13·
·微电子与基础产品 ·
电子工程师
2007年 2月
发的基于 ARM7 的 SEP3203 芯片 )为硬件平台 ,通过 DMA 方式传送 1 024个数据到 FPGA 的 RAM 模块 ,数据传输完成后开始 FFT运算 ,处理器等到 done 信号后 ,通过 CodeW arrior for ARM Developer Suite 将存储器中的运算结果导出 ,在 MATLAB 中画出图形 ,图 4 为理想值与运算值的对比。也可在 M emory区中观测具体数值 ,判断由于截断和定点运算带来的误差。
参照 FFT结构 ,旋转因子要参与蝶形运算以得到第 2个运算结果。旋转因子都是小于 1 的小数 ,所以设计中需将其定点化 ,定点化过程是将旋转因子扩大 16 384 ( 214 )倍 ,取整数部分定点化 ,存入内置 ROM。蝶形运算后将运算结果右移 14 位后作为第 2 个运算结果存入 RAM ,作为下一级的输入数据。所需 ROM 的容量应为 1 024 点 FFT所需的全部旋转因子数 ,共 512个 16 bit实部 ROM , 512个 16 bit虚部 ROM。 2. 3 蝶形运算单元
2 FFT 的硬件实现
整个设计主要包含存放数据的 RAM、存放旋转因子的 ROM、蝶形运算单元、地址产生单元以及使各个模块协调工作的控制单元 5个模块。如图 2所示。 2. 1 RAM 数据存储单元
RAM 是存储输入数据、中间运算结果以及最终运算结果的单元。每个蝶形运算的输入、输出数据均要经过 RAM 的读写操作 ,因此 , RAM 的频繁读写操作对 FFT的处理速度影响较大。为了加快 FFT的运算速度 ,需要构造双端口 RAM 来加快数据传输的吞吐量。 EP1C6Q240c内置了 20 块 4 kbit的 RAM ,将 RAM 设置在 FPGA 内部不存在驱动和 Pad延时 ,速度极快 ,而
功能仿真结果如图 3所示。从图中可以看到 ,第 1个节拍输出第 1 个操作数的地址 ,并寄存上一个蝶形运算得到的第 2个操作数 ,将其送入 RAM 输入数据总线上 ;第 2个节拍得到第 1个操作数 ,并输出第 2个操作数的地址 ;第 3个节拍得到第 2个操作数 ,并计算出第 1 个结果 ; 第 4 个节拍寄存第 1 个结果 , 送入 RAM 输入数据总线上。
[ 5 ] SUNG C H , LEE K B , JEN C W. Design and imp lementation of a scalable fast Fourier transform core [ C ] / / Proceedings of 2002 IEEE ASIA 2Pacific Conference on ASIC, Aug 628, 2002, Taipei, China. Piscataway, NJ, USA: IEEE, 2002: 2952298.
摘要 :介绍了一种用低成本 Cyclone系列 FPGA (现场可编程门阵列 )实现基于按 D IF (频率抽取 ) radix 2结构 1 024点 FFT (快速傅里叶变换 )算法的方法。本设计采用 Verilog语言编程实现 ,利用 EDA (电子设计自动化 )工具对设计进行了仿真、综合 ,并在开发板上实现板级验证 ,最后分析了整个设计的性能 ,说明在低成本 Cyclone系列上可以实现高速 FFT算法。
器 ,第 i级的地址取计数器的第 [M - 2: i ]位即可。
2. 5 控制单元
控制单元是整个系统的控制核心 ,各个功能模块
之间的地址、数据传递均通过控制单元协调工作。控
制单元记录当前蝶形运算所处的级数和个数 ,并产生 ctl信号传递给地址产生单元 ,以产生操作数和旋转因子的地址 ,地址产生单元将旋转因子的地址送入 ROM 模块的地址总线 ,直接读出旋转因子的值 ,送入蝶形运算单元 ;将操作数的地址送入控制单元 ,控制单元将两个操作数的地址连续送入 RAM 的地址总线 ,连续读取两个操作数 ,将其寄存后送入蝶形运算单元。运算结果在中央控制单元控制下连续写入 RAM。 FFT完成后 ,产生 done信号通知外部读取处理后的数据。
图 3 功能仿真结果
图 4 理想值与实际值的对比
4 结束语
本文研究用低成本 FPGA 实现 FFT的重要意义在
于电子设计师可以采用相对廉价的大规模可编程集成电路来设计灵活多变的专用芯片而不必受传统 ASIC (专用集成电路 )的固有设计模式制约。本设计全部由 Verilog实现 ,部分采用 Pipeline 并行结构 ,使用了 FPGA 内置的双端口 RAM、ROM 单元 ,加快了系统总体速度。全部电路已由功能仿真、综合、映射、布局布线和时序仿真验证 , 并成功下载到 Cyclone EP1C6Q240c。验证结果良好 ,最大误差不超过 2% , 在 75 MHz条件下 ,完成 1 024 点 FFT的运行时间为 267μs,可实现实时处理。
蝶形单元是 FFT的核心部分 ,主要包括乘法器和加法器。每级蝶形运算的输入数据和运算结果存储在相同地址单元中。在蝶形运算单元中 ,首先将两个操作数相加减 ,加法结果为第 1个运算结果 ;减法结果乘以旋转因子 ,再右移 14位 ,得出第 2个运算结果 ,同时将第 1个操作数存入 RAM ,在下 1个节拍时存入第 2 个运算结果 ,可实现部分 Pipeline结构。 2. 4 地址产生单元
仿真实验表明 ,随着可编程器件规模、速度的不断提高和成本的相对低廉 ,采用 FPGA 实现高速数字信号处理的算法具有可行性和优越性。
参考文献
[ 1 ] 胡广书. 数字信号处理 ———理论、算法与实现 [M ]. 2版. 北京 :清华大学出版社 , 2003.
[ 2 ] IFEACHOR E C, JERV IS B W. 数字信号处理实践方法 [M ]. 2版. 罗鹏飞 ,等译. 北京 :电子工业出版社 , 2004.
第 33卷第 2期 2007年 2月
EL
电子工 ECTRON IC
程 EN
师 G IN
EER
V oFl. e3b.3
No. 2 2007
基于 Cyclone系列 FPGA的 1 024点 FFT算法的实现
钱文明 ,刘新宁 ,张艳丽
(东南大学国家专用集成电路系统工程技术研究中心 ,江苏省南京市 210096)
add r_B
=
add r_A
+
N/ 2i
2
第 0级 : addr_A 随 ctl信号从 0递增至 N /2 - 1,此
时本级的地址已产生完全 ,之后再加 1溢出变成 0,成
为下一级的首地址。
第 1～M - 1 级 : 判断 addr_A 的第 [M - i - 1: 0 ]
位
,若小于
N/ 2i
地址产生单元根据控制单元输出的 ctl信号产生下一级蝶形运算所需的两个操作数和旋转因子的地址。观察 RAM 地址生成规律 :序列 x ( n)的长度为 N , M = log2 N ,则所需的级数和地址位宽都为 M ( 0～M 1 ) 。第 i级蝶形运算可分为 2i 组 ,每组有 (N /2) /2i个蝶形运算。每个蝶形运算的两个操作数的地址都有一定关系 :
1 FFT 算法
本文采用的算法是将频域 X ( k)按奇偶分开 ,故称 D IF (频率抽取 )基 2 FFT算法。
设序列 x ( n)的长度为 N ,且满足 N = 2M , M 为整数。将 N 点 DFT先分成 2 个 N /2 点 DFT,再是 4 个 N /4点 DFT,进而 8 个 N /8 点 DFT,直至 N /2 个 2 点 DFT。每分一次 ,称为一级运算。因为 M = log2 N ,所以 N 点 DFT可分成 M 级。
FFT的实现结构如图 1所示。图中有大量蝶形运算单元 ,它是 FFT算法的基本运算单位。
收稿日期 : 2006207231; 修回日期 : 2006209227。
·12·
图 1 FFT实现结构
这种结构的优点是 :在同一级运算中 ,每个蝶形的两个输入数据只对计算本蝶形有用 ,而且每个蝶形的输入、输出数据节点又同在一条水平线上 ,这就意味着计算完一个蝶形运算后 ,所得输出数据可以立即存入原输入数据所占用的存储单元 ,即同址运算 ,因此在工程实现时可以节省化的 ,无法在程序中进行循环控制 ,所以不便于扩展。