按频率抽取基2-快速傅里叶逆变换算法_MATLAB代码

合集下载

matlab自行编写fft傅里叶变换

matlab自行编写fft傅里叶变换

傅里叶变换(Fourier Transform)是信号处理中的重要数学工具,它可以将一个信号从时域转换到频域。在数字信号处理领域中,傅里叶变换被广泛应用于频谱分析、滤波、频谱估计等方面。MATLAB作为一个功能强大的数学软件,自带了丰富的信号处理工具箱,可以用于实现傅里叶变换。

在MATLAB中,自行编写FFT(Fast Fourier Transform)的过程需要以下几个步骤:

1. 确定输入信号

我们首先需要确定输入信号,可以是任意时间序列数据,例如声音信号、振动信号、光学信号等。假设我们有一个长度为N的信号x,即x = [x[0], x[1], ..., x[N-1]]。

2. 生成频率向量

在进行傅里叶变换之前,我们需要生成一个频率向量f,用于表示频域中的频率范围。频率向量的长度为N,且频率范围为[0, Fs),其中Fs 为输入信号的采样频率。

3. 实现FFT算法

FFT算法是一种高效的离散傅里叶变换算法,它可以快速计算出输入信号的频域表示。在MATLAB中,我们可以使用fft函数来实现FFT 算法,其调用方式为X = fft(x)。其中X为输入信号x的频域表示。

4. 计算频谱

通过FFT算法得到的频域表示X是一个复数数组,我们可以计算其幅

度谱和相位谱。幅度谱表示频率成分的强弱,可以通过abs(X)得到;

相位谱表示不同频率成分之间的相位差,可以通过angle(X)得到。

5. 绘制结果

我们可以将输入信号的时域波形和频域表示进行可视化。在MATLAB 中,我们可以使用plot函数来绘制时域波形或频谱图。

05-第五章-快速傅里叶变换(蝶形运算)

05-第五章-快速傅里叶变换(蝶形运算)

倒位序顺序数 nˆ
0 4 2 6 1 5 3 7
28
倒位序的变址处理(N=8)
29
同址运算(原位运算)
同址运算(原位运算)
某一列任何两个节点k 和j 的节点变量进行蝶形运算 后,得到结果为下一列k、j两节点的节点变量,而和其他 节点变量无关。这种原位运算结构可以节省存储单元, 降低设备成本。
例 运算前
11
因此,X (k ) X 1 (k ) W N kX 2 (k )只能计算出X(k)的前一半值。 后一半X(k) 值, N/2 , N/2 +1, …,N ?
利用 可得到
W W r(N 2k) N2
rk N2
N
X1( 2
k)
N 21
N 21
x1(r)WNr(2N 2k) x1(r)WNrk2
X 3(1 )x 3(0 ) W 2 1 x 3(1 )x(0)W 21x(4)x(0)WN0x(4) 这说明,N=2M的DFT可全部由蝶形运算来完成。
20
以8点为例第三次按奇偶分解
N=8按时间抽取法FFT信号流图
21
5.3.2 按时间抽取基2-FFT算法与直接计算DFT运算量的比较
由按时间抽取法FFT的信号流图可知,当N=2L时,共有 L 级 蝶形运算;每级都由 N/2 个蝶形运算组成,而每个蝶形有
23
FFT算法与直接DFT算法运算量的比较

时间抽取的基2快速傅里叶变换FFT分析与算法实现

时间抽取的基2快速傅里叶变换FFT分析与算法实现

离散时间信号的基2快速傅里叶变换FFT (时间抽取)蝶形算法实现

一、一维连续信号的傅里叶变换

连续函数f(x)满足Dirichlet (狄利克雷)条件下,存在积分变换:

正变换:2()()()()j ux F u f x e dx R u jI u π+∞

--∞==+⎰ 反变换:2()()j ux f x F u e du π+∞

-∞

=⎰

其中

()()cos(2)R u f t ut dt π+∞

-∞

=⎰

()()sin(2)I u f t ut dt π+∞

-∞

=-⎰

定义幅值、相位和能量如下:

幅度:1

22

2

()()()F u R u I u ⎡⎤⎡⎤=+⎣⎦⎣⎦ 相位:()arctan(()/())u I u R u ϕ= 能量:2

2

()()(E u R u I u =+)

二、一维离散信号的傅里叶变换

将连续信号对自变量进行抽样得到离散信号(理想冲击抽样脉冲),利用连续信号的傅里叶变换公式得到离散时间傅里叶变换DTFT ;再利用周期性进行频域抽样,得离散傅里叶变换DFT (详情参考任何一本《数字信号处理》教材)。 DFT 变换如下:

正变换:1

2/0

()(),0,1,2,1N j ux N

x F u f x e

u N π--==

=-∑。

反变换:1

2/0

1

()(),0,1,2,1N j ux N

u f x F u e

x N N

π-==

=-∑。

DFT 是信号分析与处理中的一种重要变换,因为计算机等数字设备只能存储和处理离散数据(时域、频域)。因直接计算DFT 的计算量大(与变换区间长度N 的平方成正比,当N 较大时,计算量太大),所以在快速傅里叶变换(简称FFT)出现以前,直接用DFT 算法进行谱分析和信号的实时处理是不切实际的。直到1965年发现了DFT 的一种快速算法(快速傅里叶变换,即FFT )以后,情况才发生了根本的变化。FFT 有时间抽取和频率抽取两种,下面介绍时间抽取FFT 。

FFT快速傅里叶变换(蝶形算法)详解解析

FFT快速傅里叶变换(蝶形算法)详解解析

4
5.2.1 DFT的运算量
(2)计算全部N个X(k) 值的运算量
复数乘法次数: N2
复数加法次数: N(N-1)
(3)对应的实数运算量
N 1
N 1
X (k) x(n)WNnk [Re x(n) j Im x(n)][ReWNnk j ImWNnk ]
n0
n0
N 1
{[Re x(n) ReWNnk Im x(n) ImWNnk ] n0
蝶形运算式
蝶形运算信 号流图符号
因此,只要求出2个N/2点的DFT,即X1(k)和X2(k),再 经过蝶形运算就可求出全部X(k)的值,运算量大大减少。
14
以8点为例第一次按奇偶分解
以N=8为例,
分解为2个4点
的DFT,然后
做8/2=4次蝶形
运算即可求出
WN0
所有8点X(k)的
值。
WN1
WN2
WN3
第五章 快速傅里叶变换
本章目录
直接计算DFT的问题及改进的途径 按时间抽取的基2-FFT算法 按频率抽取的基2-FFT算法 快速傅里叶逆变换(IFFT)算法 Matlab实现
2
5.1 引言
DFT在实际应用中很重要: 可以计算信号的频 谱、功率谱和线性卷积等。
直接按DFT变换进行计算,当序列长度N很 大时,计算量非常大,所需时间会很长。

(完整word版)按时间抽取的基2FFT算法分析

(完整word版)按时间抽取的基2FFT算法分析

第四章 快速傅里叶变换

有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT). 1965年,Cooley 和Tukey 提出了计算离散傅里叶变换(DFT )的快速算法,将DFT 的运算量减少了几个数量级。从此,对快速傅里叶变换(FFT )算法的研究便不断深入,数字信号处理这门新兴学科也随FFT 的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT 的多种算法,基本算法是基2DIT 和基2DIF 。FFT 在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。

快速傅里叶变换(FFT )是计算离散傅里叶变换(DFT )的快速算法。 DFT 的定义式为

)(k X =)()(1

0k R W n x N N n kn

N

∑-= 在所有复指数值kn

N W 的值全部已算好的情况下,要计算一个)(k X 需要N 次复数乘法和N -1次复数加法。算出全部N 点)(k X 共需2

N 次复数乘法和)1(-N N 次复数加法。即计算量是与2

N 成正比的。

FFT 的基本思想:将大点数的DFT 分解为若干个小点数DFT 的组合,从而减少运算量。

N W 因子具有以下两个特性,可使DFT 运算量尽量分解为小点数的DFT

运算:

(1) 周期性:k N n N kn N n

N k N W W W )()(++== (2) 对称性:k N N k N

W W -=+)

2/(

利用这两个性质,可以使DFT 运算中有些项合并,以减少乘法次数。例子:求当N =4时,X(2)的值

FFT快速傅里叶变换(蝶形算法)详解

FFT快速傅里叶变换(蝶形算法)详解

r0
l 0
l0
N 41
N 41
x3(l)WNlk 4 WNk 2 x4 (l)WNlk 4
l0
l0
X 3(k ) WNk / 2 X 4 (k )
k=0,1,…,
N 1 4
第17页,共53页。
17

X1
N 4
k
X 3 (k )
WNk/ 2 X 4 (k )
k=0,1,…,
N 1 4
直接按DFT变换进行计算,当序列长度N很大时, 计算量非常大,所需时间会很长。
FFT并不是一种与DFT不同的变换,而是DFT的一 种快速计算的算法。
第3页,共53页。
3
5.2 直接计算DFT的问题及改进的途径
DFT的运算量
设复序列x(n) 长度为N点,其DFT为
N 1
X (k) x(n)WNnk n0
FFT算法与直接DFT算法运算量的比较
N
N2
N
计算量
2 log2 N 之比M
N
N2
N
计算量
2 log2 N 之比M
24
1
4 16
4
8 64
12
16 256
32
32 1028 80
4.0 128
16 384
448 36.6
4.0 256

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

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

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

基2FFT算法是一种快速傅里叶变换(Fast Fourier Transform,FFT)的算法,在信号处理、图像处理等领域有着广泛的应用。该算法通过将N

个输入值分解成两个长度为N/2的DFT(离散傅里叶变换)来实现快速的

计算。本文将对基2FFT算法进行分析,并给出MATLAB实现。

基2FFT算法的主要思路是将输入序列分解成奇偶两个子序列,然后

分别对这两个子序列进行计算。具体步骤如下:

1.将输入序列拆分成奇数位和偶数位两个子序列。比如序列

x[0],x[1],x[2],x[3]可以拆分成x[0],x[2]和x[1],x[3]两个子序列。

2. 对两个子序列分别进行DFT计算。DFT的定义为:X[k] = Σ(x[n] * exp(-i * 2π * k * n / N)),其中k为频率的索引,N为序列长度。

3.对得到的两个DFT结果分别进行合并。将奇数位子序列的DFT结果

和偶数位子序列的DFT结果合并成一个长度为N的DFT结果。

4.重复以上步骤,直到计算结束。

基2FFT算法的时间复杂度为O(NlogN),远远小于直接计算DFT的时

间复杂度O(N^2)。这是因为基2FFT算法将问题的规模逐步减半,从而实

现了快速的计算。

下面是MATLAB中基2FFT算法的简单实现:

```matlab

function X = myFFT(x)

N = length(x);

if N == 1

X=x;%递归结束条件

return;

end

even = myFFT(x(1:2:N)); % 偶数位子序列的FFT计算

FFT_基2按时间抽取倒位序算法_C程序

FFT_基2按时间抽取倒位序算法_C程序

FFT--基2按时间抽取倒位序算法

--C程序

/*fWaveR[SAMPLENUMBER]数组初始化欲离散傅里叶变换的离散序列,SAMPLENUMBER为序列长度且SAMPLENUMBER须满足=2^m的条件,此程序最多做128个点的运算*/

#include "stdio.h"

#include

#define PI 3.1415926

#define SAMPLENUMBER 16 /*此处需要改成你实际计算序列的长度值*/

void InitForFFT();

void FFT();

float

fWaveR[SAMPLENUMBER]={1.0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},fWaveI[SAMPLENUMBER] ; /*存放存放待变换的值,离散数列一般只有实部有值,虚部为零,之所以定义虚部,因为FFT计算过程中会产生虚部*/

float w[SAMPLENUMBER]; /*存放FFT变换后的浮点型的值*/

float sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER]; /*存放 W 的实部和虚部*/

main()

{ int i;

InitForFFT(); /*计算W=e^(-j*2*PI*i/N)=0.5*[cos(2*PI*i/N)-j*sin (2*PI*i/N)]*/

FFT(fWaveR,fWaveI); /*做FFT运算*/

for ( i=0;i

{

printf(" %f",w[i]);

}

}

void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER])

6.3.3 按时间抽取的基-2 FFT算法[共4页]

6.3.3 按时间抽取的基-2 FFT算法[共4页]

数字图像处理与机器视觉——Visual C++与Matlab 实现

– 188 – 2,N r N r r r N N N N W W W W ++==- (6-30)

式(6-29)是W 矩阵中元素的某些特殊值,而式(6-30)则说明了W 矩阵元素的周期性和对称性。利用W 的周期性,DFT 运算中的某些项就可以合并;而利用W 的对称性,则可以仅计算半个W 序列。而根据这两点,我们就可以将一个长度为N 的序列分解成两个长度为N /2的序列并分别计算DFT ,这样就能节省大量的运算量。我们将在讲述常见的FFT 算法后分析节省的运算量。

这正是快速傅立叶变换(FFT: Fast Fourier Transform )的基本思路——通过将较长的序列转换成相对短得多的序列来大大减少运算量。

6.3.2 常见的FFT 算法

目前流行的大多数成熟的FFT 算法的基本思路大致可以分为两大类,一类是按时间抽取的快速傅立叶算法(Decimation In Time ,DIT-FFT ),另一类是按频率抽取的快速傅立叶算法(Decimation In Freqency ,DIF-FFT )。这两种算法思路的基本区别如下。

按时间抽取的FFT 算法是基于将输入序列f (x )分解(抽取)成较短序列,然后从这些序列的DFT 中求得输入序列F (u )的方法。由于抽取后的较短序列仍然可分,所以最终仅仅需要计算一个很短序列的DFT 。在这种算法中,我们主要关注的是当序列长度是2的整数次幂时,如何高效地进行抽取和运算的方法。

(完整word版)按时间抽取的基2FFT算法分析

(完整word版)按时间抽取的基2FFT算法分析

第四章 快速傅里叶变换

有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT). 1965年,Cooley 和Tukey 提出了计算离散傅里叶变换(DFT )的快速算法,将DFT 的运算量减少了几个数量级。从此,对快速傅里叶变换(FFT )算法的研究便不断深入,数字信号处理这门新兴学科也随FFT 的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT 的多种算法,基本算法是基2DIT 和基2DIF 。FFT 在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。

快速傅里叶变换(FFT )是计算离散傅里叶变换(DFT )的快速算法。 DFT 的定义式为

)(k X =)()(1

0k R W n x N N n kn

N

∑-= 在所有复指数值kn

N W 的值全部已算好的情况下,要计算一个)(k X 需要N 次复数乘法和N -1次复数加法。算出全部N 点)(k X 共需2

N 次复数乘法和)1(-N N 次复数加法。即计算量是与2

N 成正比的。

FFT 的基本思想:将大点数的DFT 分解为若干个小点数DFT 的组合,从而减少运算量。

N W 因子具有以下两个特性,可使DFT 运算量尽量分解为小点数的DFT

运算:

(1) 周期性:k N n N kn N n

N k N W W W )()(++== (2) 对称性:k N N k N

W W -=+)

2/(

利用这两个性质,可以使DFT 运算中有些项合并,以减少乘法次数。例子:求当N =4时,X(2)的值

基--2按频率抽取的FFT算法Decimation-in-Frequency(DIF)

基--2按频率抽取的FFT算法Decimation-in-Frequency(DIF)

(6)DIF与DIT比较1
• • • • 相同之处: (1)DIF与DIT两种算法均为原位运算。 (2)DIF与DIT运算量相同。 它们都需要 N N mF log 2 次复乘 2 N aF N log 2 次复加
DIF 与DIT 是两种等价的FFT算法
(6)DIF与DIT比较2
• 不同之处: (1)DIF与DIT两种算法结构倒过来。 DIF为输入顺序,输出乱序。运算完毕再运行 “二进制倒读”程序。 DIT为输入乱序,输出顺序。先运行“二进制 倒读”程序,再进行求DFT。 (2)DIF与DIT根本区别:在于蝶形结不同。 DIT的复数相乘出现在减法之前。 DIF的复数相乘出现在减法之后。
2.代入DFT中
X(k ) x(n)W x(n)W
n 0 N 1 2 N 1 2 n 0 nk N N 1 nk N
[ x(n)W
n 0 ( n
N 1 2
nk N
N x(n )W N 2
(n
N )k 2
]
N x(n )W N 2 n 0
同理:当k 奇数时,频率的奇数部分 N W 1, 令k 2k '1,k ' 0,1, , 1 2 N / 2 1 N nk X (k ) [ x(n) x(n )]WN k 1,3,5,7 , N 1 2 n 0 令k 2k '1代入 X (2k '1) 又 W

快速傅立叶变换(FFT)的计算机实现

快速傅立叶变换(FFT)的计算机实现

信号与系统课程设计

——FFT的计算机实现

快速傅里叶变换(FFT )的计算机实现

赖智鹏

华中科技大学电气与电子工程学院0809班U200811806

Email: 592425891@

摘要:本文是信号与系统课程的课程设计,旨在熟悉FFT 的计算过程,结合DFT 物理意义和实验结果加深对傅立叶变换的理解。文章首先用MATLAB 对一个简单信号进行FFT 仿真,得出频谱图;其次完成了FFT 的C 语言实现,结合MATLAB 作图及数据处理功能得出了C 实现下的FFT 结果;最后,讨论分析实验结果。 关键词:DFT 、基--2按时间抽取FFT 算法、MATLAB 、C 、频谱、物理意义

1. 算法描述

1) DFT 的运算量

-==

=-∑=⨯-≈⎡⎤⎣⎦1

()()()

(0,1,2...1)

2

2复数乘:次,复数加:(1)次

N nk

X k DFT x n x n W k N N n N

N N N

2) 减少运算的方法:

化长序列为短序列。如将长度为N 的序列分解为两个长度为N/2的序列

-==

=-∑=⨯-≈

⎡⎤⎣⎦ 1

2()()()(0,1,,1)02

22复数乘:次,复数加:

1)次4

2

2

4

N

nk

X k DFT x n x n W k

N N n N

N

N

N

利用

nk

N

W 的性质(注:本文中的C 程序未用到此性质)

()

⨯⨯+rN

周期性:=*

-共轭对称性:=r 可约性:=r nk nk W

W N N

n n W W N N

n n

W

W N N

3) C 程序采用基--2按时间抽取的FFT 算法

设输入序列长度为2

M

N

= (M 为正整数),将该序列按时间顺序的奇偶分解为越来越短

matlab的fft变换

matlab的fft变换

在MATLAB中进行FFT(Fast Fourier Transform)变换,可以使用fft函数。该函数将离散傅里叶变换(DFT)经过一系列变换得到简化式,使运算次数由原来的n^2次降为nlogn。

fft函数可以接受两个参数,第一个参数是待变换的序列y,第二个参数是序列的长度N。如果y为一向量,则fft返回值是y的快速傅里叶变换,与y具有相同的长度;如果y为一矩阵,则fft对矩阵的每一列向量进行快速傅里叶变换。

需要注意的是,fft变换能分辨的最高频率为采样频率的一半(即Nyquist频率),函数fft返回值是以Nyqusit频率为轴对称的,Y的前一半与后一半是复数共轭关系。

以上信息仅供参考,如果还有疑问,建议查阅专业书籍或咨询专业人士。

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

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

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

基2FFT算法(Fast Fourier Transform)是一种高效的离散傅里叶

变换(DFT)算法,其时间复杂度为O(NlogN),其中N为输入数据的大小。该算法利用了DFT的对称性质以及分治的思想,通过将一个大规模的DFT

问题分解为若干个规模较小的DFT问题来加快计算速度。

基2FFT算法的实现步骤如下:

1.输入N个复数序列x(n),其中n取值范围为0到N-1

2.如果N为1,直接返回x(n)。

3. 将x(n)分为两个子序列,分别为x_odd(n)和x_even(n),其中

x_odd(n)包含所有奇数索引的元素,x_even(n)包含所有偶数索引的元素。

4. 对x_odd(n)和x_even(n)分别进行基2FFT变换,递归地计算它们

的DFT结果。

5. 根据DFT的对称性,计算出x(k)的前半部分和后半部分的值,其

中k为频率索引,取值范围为0到N/2-1、具体计算方法是将x_odd(k)和

x_even(k)与旋转因子W_N^k相乘,可通过以下公式计算:

x(k) = x_even(k) + W_N^k * x_odd(k)

x(k+N/2) = x_even(k) - W_N^k * x_odd(k)

其中W_N^k = e^(-j*2*pi*k/N)为旋转因子。

6.返回x(k)作为输出结果。

基2FFT算法的MATLAB实现如下:

```matlab

function X = myfft(x)

N = length(x);

if N == 1

X=x;%如果序列长度为1,直接返回原始序列

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

function x=MyIFFT_FB(y)

%MyIFFT_TB:My Inverse Fast Fourier Transform Time Based

%按频率抽取基2-傅里叶逆变换算法

%input:

% y -- 傅里叶正变换结果,1*N的向量

%output:

% x -- 逆变换结果,1*N的向量

%参考文献:

% /view/fea1e985b9d528ea81c779ee.html

N=length(y);

x=conj(y); %求共轭

x=MyFFT_FB(x);%求FFT

x=conj(x);%求共轭

x=x./N;%除以N

end

%% 内嵌函数====================================================== function y=MyFFT_FB(x,n)

%MYFFT_TB:My Fast Fourier Transform Frequency Based

%按频率抽取基2-fft算法

%input:

% x -- 输入的一维样本

% n -- 变换长度,缺省时n=length(x) 当n小于x数据长度时,x数据被截断到第n个数据% 当n大于时,x数据在尾部补0直到x 含n个数据

%output:

% y -- 1*n的向量,快速傅里叶变换结果

%variable define:

% N -- 一维数据x的长度

% xtem -- 临时储存x数据用

% m,M -- 对N进行分解N=2^m*M,M为不能被2整除的整数

% two_m -- 2^m

% adr -- 变址,1*N的向量

% l -- 当前蝶形运算的级数

% W -- 长为N/2的向量,记录W(0,N),W(1,N),...W(N/2-1,N)

% d -- 蝶形运算两点间距离

% t -- 第l级蝶形运算含有的奇偶数组的个数

% mul -- 标量,乘数

% ind1,ind2 -- 标量,下标

% tem -- 标量,用于临时储存

%参考文献:

% /view/fea1e985b9d528ea81c779ee.html

%% 输入参数个数检查

msg=nargchk(1,2,nargin);

error(msg);

%% 输入数据截断或加0

N=length(x);

if nargin==2

if N

xtem=x;

x=zeros(1,n);

x(1:N)=xtem;

N=n;

else % 截断

xtem=x;

x=xtem(1:n);

N=n;

end

end

%% 对N进行分解N=2^m*M

[m,M]=factorize(N);

two_m=N/M;

%% 变换

if m~=0

%% 如果N可以被2整除

adr=address(m,M,two_m);

y=x; % 蝶形运算级数l=m 时

%% 计算W向量

W=exp(-2*pi*i* ( 0:N/2-1 ) /N);

%% 蝶形运算

d=N/2;

t=1;

for l=1:m

% 加

for ii=0:t-1

ind1=ii*2*d+1;

ind2=ind1+d;

for r=0:d-1

tem=y(ind1)+y(ind2);

y(ind2)=y(ind1)-y(ind2);

y(ind1)=tem;

ind1=ind1+1;

ind2=ind2+1;

end

end

% 乘

for r=0:d-1

mul=W(r*t+1);

for ii=0:t-1

y(ii*2*d+d+1+r) = y(ii*2*d+d+1+r)*mul;

end

end

d=d/2;t=t*2;

end

%% 直接傅里叶变换

if M~=1 % N 分解含有非2因数M时,对y中每M个数据做直接傅里叶变换for ii=1:two_m

y((ii-1)*M+1 : ii*M ) = DDFT( y((ii-1)*M+1 : ii*M ) );

end

end

%% 变址输出

y=y(adr+1);

else

%% 如果N 不能被2整除

y=DDFT(x);

end

end

function y=DDFT(x)

%% 直接离散傅里叶变换

%input:

% x -- 样本数据,N维向量

%output:

% y -- N维向量

%参考文献:

% 结构动力学,克拉夫,P82

% variable define

% s -- sum,用于求和

N=length(x);

y=zeros(size(x));

for n=1:N

s=0;

for m=1:N

s=s+x(m)*exp( -i*2*pi*(m-1)*(n-1)/N );

end

y(n)=s;

end

end

function [m,M]=factorize(N)

%% 对N分解

m=0;

while true

if mod(N,2)==0

m=m+1;

N=N/2;

else

M=N;

break;

end

end

end

function adr=address(m,M,two_m)

%% 变址

% b -- 2^m * m 的矩阵,用来存储二进制数据

% ds -- 数,公差

adr=zeros(two_m,M);

b=de2bi(0:two_m-1,m);%转换为2进制注:matlab中二进制[0 1 1]=6 b=b(:,end:-1:1);% 逆序

adr(:,1)=bi2de(b);%2进制转换为10进制

if M~=1

ds=two_m;

adr=adr(:,1)*ones(1,M);

adr=adr+ds*ones(size(adr,1),1)*(0:M-1);

adr=reshape(adr',1,[]);

end

end

联系方式:matrixsuper@

相关文档
最新文档