C#编写的FFT

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.IO;

namespace TM_Device

{

/*********************************************************************

快速福利叶变换C#函数

函数简介：此函数是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依赖硬件。此函数采用联合体的形式表示一个复数，输入为自然顺序的复

数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的

复数

使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的应该为2的N次方，不满足此条件时应在后面补0

作者：彭天明

QQ:469270355

Email:pengtianming2008@https://www.wendangku.net/doc/d26285860.html,

函数调用：FFT(s);

时间：2010-2-20

版本：Ver1.0

参考文献：

**********************************************************************/ class TM_fft

{

///

public class complex

{

///

public double r;

///

public double i;

};

///

///

///

///

private complex complex_subtract(complex a, complex b)

{

complex sum = new complex();

sum.r = a.r - b.r;

sum.i = a.i - b.i;

return sum;

}

///

///

///

///

private complex complex_plus(complex a, complex b)

{

complex sum = new complex();

sum.r = a.r + b.r;

sum.i = a.i + b.i;

return sum;

}

/*******************************************************************

函数原型：struct compx EE(struct compx b1,struct compx b2)

函数功能：对两个复数进行乘法运算

输入参数：两个以联合体定义的复数a,b

输出参数：a和b的乘积，以联合体的形式输出

*******************************************************************/

///

///

///

///

private complex Diexing_operator(complex a, complex b)

{

complex c = new complex();

c.r = a.r * b.r - a.i * b.i;

c.i = a.r * b.i + a.i * b.r;

return c;

}

/*****************************************************************

函数原型：void FFT(struct compx *xin,int N)

函数功能：对输入的复数组进行快速傅里叶变换（FFT）

输入参数：*xin复数结构体组的首地址指针，struct型

*****************************************************************/

///

///要变化的数组

///定义福利叶变换的点数(必须确定为2的整数倍)

///

public complex[] function_fft(complex[] xin, int FFT_N)

{

double PI = 3.1415926535897932384626433832795028841971; //定义圆周率值

//complex[] Da_S = new complex[FFT_N]; //FFT输入和输出：从Da_S[0]开始存放，根据大小自己定义

int f, m, nv2, nm1, i, k, l, j = 0;

complex u = new complex();

complex w = new complex();

complex t = new complex();

nv2 = FFT_N / 2; //变址运算，即把自然顺序变成倒位序，采用雷德算法

nm1 = FFT_N - 1;

for (i = 0; i < nm1; i++)

{

if (i < j) //如果i

{

t = xin[j];

xin[j] = xin[i];

xin[i] = t;

}

k = nv2; //求j的下一个倒位序

while (k <= j) //如果k<=j,表示j的最高位为1

{

j = j - k; //把最高位变成0

k = k / 2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0

}

j = j + k; //把0改为1

}

{

int le, lei, ip; //FFT运算核，使用蝶形运算完成FFT运算 f = FFT_N;

for (l = 1; (f = f / 2) != 1; l++) //计算l的值，即计算蝶形级数 ;

for (m = 1; m <= l; m++) // 控制蝶形结级数

{ //m表示第m级蝶形，l为蝶形级总数l=log（2）N

le = 2 << (m - 1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点

lei = le / 2; //同一蝶形结中参加运算的两点的距离

u.r = 1.0; //u为蝶形结运算系数，初始值为1

u.i = 0.0;

w.r = Math.Cos(PI / lei); //w为系数商，即当前系数与前一个系数的商

w.i = -Math.Sin(PI / lei);

for (j = 0; j <= lei - 1; j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结

{

for (i = j; i <= FFT_N - 1; i = i + le) //控制同一蝶形结运算，即计算系数相同蝶形结

{

ip = i + lei; //i，ip分别表示参加蝶形运算的两个节点

t = Diexing_operator(xin[ip], u); //蝶形运算，详见公式

xin[ip].r = xin[i].r - t.r;

xin[ip].i = xin[i].i - t.i;

xin[i].r = xin[i].r + t.r;

xin[i].i = xin[i].i + t.i;

}

u = Diexing_operator(u, w); //改变系数，进行下一个蝶形运算

}

}

}

return xin;

}

///

///输入的文件

///输出的文件

public double[,] Juzhen_Readdatafromfile(string inputfile)

{

FileStream stream1 = new FileStream(inputfile, FileMode.Open, FileAccess.Read);

byte[] buffer1 = new byte[(int)stream1.Length];

stream1.Read(buffer1, 0, buffer1.Length);

stream1.Close();

string myTxt = Encoding.UTF8.GetString(buffer1);

string[] Rows_data = myTxt.Split('\n');

string[] Row_data = Rows_data[0].Split(',');

int Rows = Rows_data.Length;

int Cols = Row_data.Length-1;

double[,] double_data = new double[Rows, Cols];

for (int ix = 0; ix < Rows_data.Length; ix++)

{

try

{

string[] Each_data = Rows_data[ix].Split(',');

for (int iy = 0; iy < Each_data.Length-1; iy++)

{

double_data[ix, iy] = Convert.ToDouble(Each_data[iy]);

}

}

catch (Exception ex)

{

ex.ToString();

}

}

return double_data;

}

///

///输入的文件

///输出的文件

public double[,] Juzhen_Readdatafromfile(string inputfile,char FG_char)

{

FileStream stream1 = new FileStream(inputfile, FileMode.Open, FileAccess.Read);

byte[] buffer1 = new byte[(int)stream1.Length];

stream1.Read(buffer1, 0, buffer1.Length);

stream1.Close();

string myTxt = Encoding.UTF8.GetString(buffer1);

string[] Rows_data = myTxt.Split('\n');

string[] Row_data = Rows_data[0].Split( FG_char );

int Rows = Rows_data.Length;

int Cols = Row_data.Length - 1;

double[,] double_data = new double[Rows, Cols];

for (int ix = 0; ix < Rows_data.Length; ix++)

{

try

{

string[] Each_data = Rows_data[ix].Split(FG_char);

for (int iy = 0; iy < Each_data.Length - 1; iy++)

{

double_data[ix, iy] = Convert.ToDouble(Each_data[iy]);

}

}

catch (Exception ex)

{

ex.ToString();

}

}

return double_data;

}

///

public double[] YW_to_Mean(double[] SZ)

{

int count = SZ.Length;

double [] C_SZ = new double [count];

///-------求和

for (int i = 0; i < count; i++)

{

C_SZ[i] = SZ[i] / count;

}

return C_SZ;

}

///

public double[] YW_Sumto_Mean(double[] SZ)

{

int count = SZ.Length;

double[] C_SZ = new double[count];

///-------求和

double Sum_All=new double() ;

for (int i = 0; i < count; i++)

{

Sum_All +=SZ[i];

}

for (int i = 0; i < count; i++)

{

C_SZ[i] = SZ[i] / Sum_All;

}

return C_SZ;

}

///

public double[ ,] Juzhen_zhuanzhi(double[,] MatrixA)

{

int M = MatrixA.GetLength(0);

int N = MatrixA.GetLength(1);

double[,] MatrixB = new double[N,M];

for (int i = 0; i < N; i++)

{

for (int j = 0; j < M; j++)

{

MatrixB[i,j] = MatrixA[j,i];

}

}

return MatrixB;

}

///

///

public complex[] YW_to_ComplexYW(double[] SZ )

{

int count=SZ.Length ;

complex[] C_SZ = new complex[count];

for (int i = 0; i < count;i++ )

{

complex d = new complex();

d.r = SZ[i];

d.i = 0;

C_SZ[i] = d;

}

return C_SZ;

}

///

///

public double [] ComplexYW_to_YW(complex[] SZ)

{

int count = SZ.Length;

double [] C_SZ = new double [count];

for (int i = 0; i < count; i++)

{

C_SZ[i] = SZ[i].r;

}

return C_SZ;

}

} }

Matlab中的FFT使用说明

FFT是Fast Fourier Transform（快速傅里叶变换）的简称，FFT算法在MATLAB 中实现的函数是Y=fft（x,n）。刚接触频谱分析用到FFT时，几乎都会对MATLAB 的fft函数产生一些疑惑，下面以看一个例子（根据MATLA帮助修改）。 Fs = 2000; % 设置采样频率 T = 1/Fs; % 得到采用时间 L = 1000; % 设置信号点数，长度1 秒 t = (0:L-1)*T; % 计算离散时间， % 两个正弦波叠加 f1 = 80; A1 = 0.5; % 第一个正弦波100Hz，幅度0.5 f2 = 150; A2 = 1.0 ; % 第2个正弦波150Hz，幅度 1.0 A3 = 0.5; % 白噪声幅度； x = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t); % 产生离散时间信号； y = x + A3*randn(size(t)); % 叠加噪声； % 时域波形图 subplot(2,1,1) plot(Fs*t(1:50),x(1:50)) title('Sinusoids Signal') xlabel('time (milliseconds)') subplot(2,1,2) plot(Fs*t(1:50),y(1:50)) title('Signal Corrupted with Zero-Mean Random Noise') xlabel('time (milliseconds)') NFFT = 2A nextpow2(L); % 设置FFT点数，一般为2 的N次方，如1024，512 等Y = fft(y,NFFT)/L; % 计算频域信号， f = Fs/2*linspace(0,1,NFFT/2+1); %频率离散化，fft后对应的频率是-Fs/2到Fs/2,由NFFT个离散频点表示 % 这里只画出正频率； % Plot single-sided amplitude spectrum. figure; plot(f,2*abs(Y(1:NFFT/2+1))); % fft 后含幅度和相位，一般观察幅度谱，并把负频率加上去， title('Single-Sided Amplitude Spectrum of y(t)') xlabel('Frequency (Hz)')

Simulink下的频谱分析方法及matlab的FFT编程

Simulink下的频谱分析方法 实现功能： 信号发生器一个信号输入，实时显示其频谱分析 调用模块： 信号源（Signal Processing Blockset -> Signal Processing Sources -> Sine Wave） Tip 1：不能用连续的信号源 频谱观察窗（Signal Processing Blockset -> Signal Processing Sources -> Spectrum Scope）Tip 2: 不能用普通的观察窗 Tip 3：必须构上设置中的Buffer input. Buffer size 越大越精细。 Tip 4: 剩下的tips读帮助。 连接关系： 如下图所示 原理框图实验结果：

输出示意图－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 实现功能： 从Workspace读取一组数，进行频谱分析 调用模块： From Workspace Tip 1: 采样时间不能用0，即必须使用离散模式 Tip 2: 从其他模型中Scope保存出来的“Structure with time”的数据可以直接用频谱观察窗（同上一功能） －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ －－－－－－－－－－－－－－－－－－－－－－－－－－－－－－ 实现功能： 从dSPACE读取一组数，进行频谱分析 实现方法：

1. 从dSPACE读数保存成文件，数据导入Workspace（过程略） 2. 采用从其他模型的Scope保存数据为“Structure with time”的方式构建一个结构变量ScopeData1 3. 使用以下代码将dSPACE数据dscapture拷贝到结构变量ScopeData1中 %% =[0::]; %纯粹为占位，19157为dSPACE保存数据长度 for i=1:19157 end %% 4. 采用下图中的模型进行频谱分析 实验结果： 通过以上方法对单轴压电加速度传感器进行灵敏度分析，下图分别为采用dSPACE和直接利用示波器分析的结果对比。

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

按时间抽取的基2FFT 算法分析及MATLAB 实现 一、DIT-FFT 算法的基本原理 基2FFT 算法的基本思想是把原始的N 点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT ，再进行适当的组合，得到原N 点序列的DFT ，最终达到减少运算次数，提高运算速度的目的。 按时间抽取的基2FFT 算法，先是将N 点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT 运算，求出与之对应的X1(k)和X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N 点序列的DFT 。只要N 是2的整数次幂，这种分解就可一直进行下去，直到其DFT 就是本身的1点时域序列。 图1 DIT-FFT 蝶形运算流图 二、DIT-FFT 算法的运算规律及编程思想 1.原位计算 对N=M 2点的FFT 共进行M 级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，经过M 级运算后，原来存放输入序列数据的N 个存储单元中可依次存放X(k)的N 个值，这种原位(址)计算的方法可节省大量内存。 2.旋转因子的变化规律 N 点DIT ―FFT 运算流图中，每个蝶形都要乘以旋转因子p W N ，p 称为旋转因子的指数。例如N ＝8 ＝3 2 时各级的旋转因子： 第一级：L=1， 有1个旋转因子：p W N =J /4W N =J 2L W J=0 第二级：L=2，有2个旋转因子：p W N =J /2W N =J 2L W J=0,1 第三级：L=3，有4个旋转因子：p W N =J W N =J 2L W J=0,1,2,3 对于N ＝M 2的一般情况，第L 级共有1 -L 2个不同的旋转因子： p W N =J 2L W J=0,1,2,… ,1 -L 2－1 L 2=M 2×M -L 2 = N ·M -L 2 故： 按照上面两式可以确定第L 级运算的旋转因子

Matlab编程实现FFT变换.

Matlab编程实现FFT变换及频谱分析的程序代码 内容 1．用Matlab产生正弦波,矩形波,以及白噪声信号，并显示各自时域波形图 2．进行FFT变换，显示各自频谱图，其中采样率，频率、数据长度自选 3．做出上述三种信号的均方根图谱,功率图谱,以及对数均方根图谱 4．用IFFT傅立叶反变换恢复信号，并显示恢复的正弦信号时域波形图 源程序 %*************************************************************** **********% % FFT实践及频谱分析% %*************************************************************** **********% %*************************************************************** **********% %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值'); title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱

MATLAB中FFT结果的物理意义

FFT结果的物理意义 最近正在做一个音频处理方面的项目，以前没有学过fft，只是知道有这么个东西，最近这一用才发现原来欠缺这么多，最基本的，连fft的输入和输出各自代表什么都不知道了，终于在网上查到这样的一点资料，得好好保存了，也欢迎大家分享。 FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。现在圈圈就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此罗嗦了。 采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N（ps：横坐标第n个点对应的频率值Fn的计算公式。整个横坐标代表了采样频率Fs，被分为N点。故其频率分辨率为Fs/N）。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。 假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。 好了，说了半天，看着公式也晕，下面圈圈以一个实际的信号来做说明。假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)式中cos参数为弧度，所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？我们来看看FFT的结果的模值如图所示。

MATLAB关于FFT频谱分析的程序

MATLAB关于FFT频谱分析的程序 %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值');

title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y); figure(1); subplot(233); plot(f,sq); xlabel('频率(Hz)'); ylabel('均方根谱'); title('正弦信号y=2*pi*10t均方根谱'); grid; %求功率谱 power=sq.^2; figure(1); subplot(234); plot(f,power); xlabel('频率(Hz)'); ylabel('功率谱'); title('正弦信号y=2*pi*10t功率谱'); grid; %求对数谱 ln=log(sq); figure(1); subplot(235); plot(f,ln);

利用MATLAB实现信号DFT的计算

07级电信（2）班 刘坤洋 24 实验一 利用MATLAB 实现信号DFT 的计算 一、实验目的： 1、熟悉利用MATLAB 计算信号DFT 的方法 2、掌握利用MATLAB 实现由DFT 计算线性卷积的方法 二、实验设备：电脑、matlab 软件 三、实验内容： 1、练习用matlab 中提供的内部函数用于计算DFT （1） fft （x ），fft （x ，N ），ifft （x ），ifft （x ，N ）的含义及用法 （2） 在进行DFT 时选取合适的时域样本点数N 请举例，并编程实现 题目： 源程序： >> N=30; %数据的长度 >>L=512; %DFT 的点数 >>f1=100; f2=120; >>fs=600; %抽样频率 >>T=1/fs; %抽样间隔 >>ws=2*pi*fs; >>t=(0:N-1)*T; >>f=cos(4*pi*f1*t)+cos(4*pi*f2*t); >>F=fftshift(fft(f,L)); >>w=(-ws/2+(0:L-1)*ws/L)/(2*pi); >>hd=plot(w,abs(F)); >>ylabel('幅度谱') >> xlabel('频率/Hz') 的频谱 分析利用)π4cos()π4cos()(DFT 21t f t f t x +=Hz 600,Hz 120,Hz 10021===s f f f

>> title('my picture') 结果图： （3） 在对信号进行DFT 时选择hamming 窗增加频率分辨率 请举例，并编程实现 题目： 源程序：>> N=50; %数据的长度 >>L=512; %DFT 的点数 >>f1=100;f2=150; >>fs=600; %抽样频率 >>T=1/fs; %抽样间隔 >>ws=2*pi*fs; >>t=(0:N-1)*T; >>f=cos(4*pi*f1*t)+0.15*cos(4*pi*f2*t); 的频谱 分析利用)π4cos(15.0)π4cos()(DFT 21t f t f t x +=Hz 600,Hz 150,Hz 10021===s f f f

实验二 FFT算法的MATLAB实现

班级：学号：姓名 实验二FFT算法的MATLAB实现 （一）实验目的： （1）掌握用matlab进行FFT在数字信号处理中的高效率应用。 （2）学习用FFT对连续信号和时域离散信号进行谱分析。 （二）实验内容及运行结果： 题1：若x(n)=cos(nπ/6)是一个N=12的有限序列，利用MATLAB计算它的DFT 并进行IDFT变换同时将原图与IDFT变换后的图形进行对比。当求解IFFT变换中，采样点数少于12时，会产生什么问题。 程序代码： N=12; n=0:11; Xn=cos(n*pi/6); k=0:11; nk=n'*k; WN=exp(-j*2*pi/N) WNnk=WN.^nk XK=Xn*WNnk; figure(1) stem(Xn) figure(2) stem(abs(XK)) 运行结果：

IFFT变换中，当采样点数少于12时图像如下图显示：

分析：由图像可以看出，当采样点数小于12时，x(n)的频谱不变，周期为6，而XK 的频谱图发生改变。 题2：对以下序列进行谱分析 132()()103()8470x n R n n n x n n n =+≤≤?? =-≤≤??? 其他n 选择FFT 的变换区间N 为8和16点两种情况进行频谱分析，分别打印其幅频特 性曲线并进行对比、分析和讨论。 ㈠ 程序代码： x=ones(1,3);nx=0:2; x1k8=fft(x,8); F=(0:length(x1k8)-1)'*2/length(x1k8); %进行对应的频率转换 stem(f,abs(x1k8));%8点FFT title('8点FFTx_1(n)'); xlabel('w/pi'); ylabel('幅度'); N=8时：

MATLAB中FFT的使用方法

MATLAB中FFT的使用方法 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) →Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数

n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果：

用MATLAB进行FFT频谱分析

用MATLAB 进行FFT 频谱分析 假设一信号： ()()292.7/2cos 1.0996.2/2sin 1.06.0+++=t t R ππ 画出其频谱图。 分析： 首先，连续周期信号截断对频谱的影响。 DFT 变换频谱泄漏的根本原因是信号的截断。即时域加窗，对应为频域卷积，因此，窗函数的主瓣宽度等就会影响到频谱。 实验表明，连续周期信号截断时持续时间与信号周期呈整数倍关系时，利用DFT 变换可以得到精确的模拟信号频谱。举一个简单的例子： ()ππ2.0100cos +=t Y 其周期为0.02。截断时不同的持续时间影响如图一.1：（对应程序shiyan1ex1.m ） 图 错误!文档中没有指定样式的文字。.1 140.0160.0180.02 截断时，时间间期为周期整数倍，频谱图 0.0250.03 20 40 60 80 100 截断时，时间间期不为周期整数倍，频谱图

其次，采样频率的确定。 根据Shannon 采样定理，采样带限信号采样频率为截止频率的两倍以上，给定信号的采样频率应>1/7.92，取16。 再次，DFT 算法包括时域采样和频域采样两步，频域采样长度M 和时域采样长度N 的关系要符合M ≧N 时，从频谱X(k)才可完全重建原信号。 实验中信号R 经采样后的离散信号不是周期信号，但是它又是一个无限长的信号，因此处理时时域窗函数尽量取得宽一些已接近实际信号。 实验结果如图一.2：其中，0点位置的冲激项为直流分量0.6造成（对应程序为shiyan1.m ） 图 错误!文档中没有指定样式的文字。.2 ?ARMA （Auto Recursive Moving Average ）模型： 将平稳随机信号x(n)看作是零均值，方差为σu 2的白噪声u(n)经过线性非移变系统H(z)后的输出，模型的传递函数为 020406080100120140160180200 0.4 0.50.60.7 0.800.050.10.150.20.250.30.350.40.450.5 50100 150

基于matlab的FFT算法程序设计

数字通信课程设计报告书 课题名称 基于matlab 的FFT 算法程序设计 姓 名 学 号 院 系 物理与电信工程系 专 业 电子信息工程 指导教师 2010年 01 月15日 ※※※※※※※※※ ※ ※ ※※ ※※ ※※ ※※※※※ ※※ 2007级数字通信 课程设计

基于matlab的FFT算法程序设计 0712401-36 李晔 （湖南城市学院物理与电信工程系通信工程专业,益阳,413000） 一、设计目的 1.通过该设计，进一步了解MATLAB软件。 2.通过该设计，进一步熟悉MATLAB的语法规则和编辑方式。 3.通过该设计，掌握傅里叶变换的含义和方法。 二、设计的主要要求 掌握Fourier变换，解了关于MATLAB软件在数字信号处理方面的应用，熟悉MATLAB的语法规则和编程。用MATLAB实现快速Fourier变换。 三、整体设计方案 对信号x=sin(2*pi*f0*t)进行频谱分析，用MATLAB仿真。选取抽样频率为fs=100Hz,依照下列条件用MATLAB软件对信号xt进行傅里叶变换y=fft(xt,N)并绘制频谱图，观察所产生的六幅频谱图进行对比，并进行分析。 四、程序设计 fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率

%生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(321); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 m=length(y); f=(0:m/2-1)'*fs/m;%进行对应的频率转换 figure(1); subplot(322); plot(f,mag(1:m/2));%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值'); title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y); figure(1); subplot(323); plot(f,sq(1:m/2)); xlabel('频率(Hz)'); ylabel('均方根谱');

MATLAB中FFT使用详解

MATLAB中FFT使用详解 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MA TLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) → Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅

xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果： fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz

利用MATLAB编写FFT快速傅里叶变换

一、实验目的 1.利用MATLAB 编写FFT 快速傅里叶变换。 2.比较编写的myfft 程序运算结果与MATLAB 中的FFT 的有无误差。 二、实验条件 PC 机，MATLAB7.0 三、实验原理 1. FFT （快速傅里叶变换）原理： 将一个N 点的计算分解为两个N/2点的计算，每个N/2点的计算再进一步分解为N/4点的计算，以此类推。根据DFT 的定义式，将信号x[n]根据采样号n 分解为偶采样点和奇采样点。设偶采样序列为y[n]=x[2n]，奇采样序列为z[n]=x[2n+1]。 上式中的k N W -为旋转因子N k j e /2π-。下式则为y[n]与z[n]的表达式： 2. 蝶形变换的原理：

下图给出了蝶形变换的运算流图，可由两个N/2点的FFT （Y[k]和Z[k]得出N 点FFT X[k]）。同理，每个N/2点的FFT 可以由两个N/4点的FFT 求得。按这种方法，该过程可延迟后推到2点的FFT 。 下图为N=8的分解过程。图中最右边的为8个时域采样点的8点FFTX[k]，由偶编号采样点的4点FFT 和奇编号采样点的4点得到。这4点偶编号又由偶编号的偶采样点的2点FFT 和奇编号的偶采样点的2点FFT 产生。相同的4点奇编号也是如此。依次往左都可以用相同的方法算出，最后由偶编号的奇采样点和奇编号的偶采样点的2点FFT 算出。图中没2点FFT 成为蝶形，第一级需要每组一个蝶形的4组，第二级有每组两个蝶形的两组，最后一级需要一组4个蝶形。 四、实验内容 1.定义函数disbutterfly ，程序根据FFT 的定义：]2 [][][N n x n x n y + +=、n N W N n x n x n z -+ -=])2 [][(][，将序列x 分解为偶采样点y 和奇采样点z 。

用matlab进行fft谐波分析

FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 现在就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此罗嗦了。 采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。 假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs 为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。 假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。 由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。 好了，说了半天，看着公式也晕，下面以一个实际的信号来做说明。 假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V 的交流信号，以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下： S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180) 式中cos参数为弧度，所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？我们来看看FFT的结果的模值如图所示。

用matlab实现fft算法

A1=str2double(get(handles.edit8,'String')); A2=str2double(get(handles.edit9,'String')); F1=str2double(get(handles.edit10,'String')); F2=str2double(get(handles.edit11,'String')); Fs=str2double(get(handles.edit12,'String')); N=str2double(get(handles.edit13,'String')); t=[0:1/Fs:(N-1)/Fs]; x=A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t); %信号x的离散值 axes(handles.axes1) %在axes1中作原始信号图 plot(x); grid on m=nextpow2(x);N=2^m; % 求x的长度对应的2的最低幂次m if length(x)

实验三用FFT对信号进行频谱分析和MATLAB程序

实验三 用FFT 对信号进行频谱分析 一 实验目的 1 能够熟练掌握快速离散傅立叶变换的原理及应用FFT 进行频谱分析的基本方法； 2了解用FFT 进行频谱分析可能出现的分析误差及其原因； 二 实验原理 1.用DFT 对非周期序列进行谱分析 单位圆上的Z 变换就是序列的傅里叶变换，即 ()()j j z e X e X z ωω== （3-1） ()j X e ω是ω的连续周期函数。对序列()x n 进行N 点DFT 得到()X k ，则()X k 是在区间[]0,2π上对()j X e ω的N 点等间隔采样，频谱分辨率就是采样间隔 2N π。因此序列的傅里叶变换可利用DFT （即FFT ）来计算。 用FFT 对序列进行谱分析的误差主要来自于用FFT 作频谱分析时，得到的是离散谱，而非周期序列的频谱是连续谱，只有当N 较大时，离散谱的包络才能逼近连续谱，因此N 要适当选择大一些。 2.用DFT 对周期序列进行谱分析 已知周期为N 的离散序列)(n x ，它的离散傅里叶级数DFS 分别由式（3-2）和（3-3） 给出： DFS ： ∑-=-=1 2)(1N n kn N j k e n x N a π ， n =0,1,2,…,N -1 （3-2） IDFS ： ∑-==1 02)(N k kn N j k e a n x π ， n =0,1,2,…,N -1 （3-3） 对于长度为N 的有限长序列x (n )的DFT 对表达式分别由式（3-4）和（3-5）给出： DFT ： ∑-=-=1 02)()(N n kn N j e n x k X π ， n =0,1,2,…,N -1 （3-4） IDFT ： ∑-==1 02)(1)(N k kn N j e k X N n x π ， n =0,1,2,…,N -1 （3-5） FFT 为离散傅里叶变换DFT 的快速算法，对于周期为N 的离散序列x (n )的频谱分析便可由式（3-6）和（3-7）给出：

基于MATLAB的FFT算法实现(论文)

基于MATLAB的FFT算法实现 摘要 MATLAB软件是目前全世界范围内非常流行的具有很强的科学计算和图形界面的软件系统。利用MATLAB的强大运算功能，可以解决数字信号处理过程中遇到的许多问题。本文给出了基于MATLAB软件实现信号DFT变换和FFT频谱分析的方法。利用MATLAB软件方法，使得设计方便、快捷，大大减轻了工作量。并且，在信号DFT变换中可以清楚得看到DFT变换结果和截取长度之间的关系。通过编程仿真可以得到序列的幅频特性曲线，便于对信号进行谱分析。 FFT（Fast Fourier Transformation），即为快速傅立叶变换，是离散傅立叶变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 在实际应用中，FFT是最常见的数字信号处理算法，它在各种数字信号处理系统中扮演重要的角色。在信号处理过程中。频域分析往往比时域分析方便和高效，FFT是时域和频域转换的基本运算。 关键词：FFT算法；MATLAB；数字信号处理；频谱分析

THE FFT ALGORITHM BASED ON MATLAB ABSTRACT MATLAB software is very popular around the world have a strong scientific computing and graphic interface of the software system. Using the powerful operation function of MATLAB, can solve many problems encountered in the process of digital signal processing. In this paper, based on DFT MATLAB software to realize signal transform and FFT spectrum analysis method. Using MATLAB software method, makes the design of convenient, quick, greatly reduce the workload. And the signal DFr transform can be clearly seen in the DFT transform results and clipping of the relationship between the length. Sequences can be obtained by programming the simulation of the amplitude frequency characteristic curve, facilitate the signal spectrum analysis. FFT (Fast Fourier changed), which is Fast Fourier transform, is a Fast algorithm of discrete Fourier transform, it is according to the odd and even of discrete Fourier transform, virtual and real features, the discrete Fourier transform algorithm was improved. It the theory of Fourier transform and found nothing new, but the application in computer systems or digital system discrete Fourier transform, can be said to be into a big step. In practical applications, the FFT is the most common form of digital signal processing algorithm, it is play an important role in all kinds of digital signal processing system. In the process of signal processing. Frequency domain analysis than time-domain analysis is convenient and efficient, FFT is the basic operation of the time domain and frequency domain transformation. Key words:FFT algorithm; MATLAB; Digital signal processing; Spectrum analysis