一 : 用快速傅里叶变换对信号进行频谱分析
实验二 用快速傅里叶变换对信号进行频谱分析
一、实验目的
1．理解离散傅里叶变换的意义；
2．掌握时域采样率的确定方法；
3．掌握频域采样点数的确定方法；
4．掌握离散频率与模拟频率之间的关系；
5．掌握离散傅里叶变换进行频谱分析时，各参数的影响。(www.61k.com)
二、实验原理
序列的傅里叶变换结果为序列的频率响应，但是序列的傅里叶变换是频率的连续函数，而且在采用计算机计算时，序列的长度不能无限长，为了便于计算机处理，作如下要求：序列x(n)为有限长，n从0～N－1，再对频率ω在0～2π范围内等间隔采样，采样点数为N，采样间隔为2π/N。第k个采样点对应的频率值为2πk/N。可得离散傅里叶变换及其逆变换的定义为
N?1
n?02?knNX(k)??x(n)e?j (1)
(2) j1N?1
x(n)??X(k)eNk?02?nkN
如果把一个有限长序列看作是周期序列的一个周期，则离散傅里叶变换就是傅里叶级数。离散傅里叶变换也是周期的，周期为N。
数字频率与模拟频率之间的关系为
??2?f/fs，即f?
则第k个频率点对应的模拟频率为 ?fs?? (3) 2?2?Ts
fk?kf2?k1k???s (4) N2?TsNTsN
在用快速傅里叶变换进行频谱分析时，要确定两个重要参数：采样率和频域采样点数，采样率可按奈奎斯特采样定理来确定，采样点数可根据序列长度或频率分辨率△f来确定
fsf??f，则N?s (5) N?f
用快速傅里叶变换分析连续信号的频谱其步骤可总结如下：
（1）根据信号的最高频率，按照采样定理的要求确定合适的采样频率fs；
（2）根据频谱分辨率的要求确定频域采样点数N，如没有明确要求频率分辨率，则根据实际需要确定频率分辨率；
（3）进行N点的快速傅里叶变换，最好将纵坐标根据帕塞瓦尔关系式用功率来表示，
傅里叶变换对 用快速傅里叶变换对信号进行频谱分析
横坐标根据式(7-21)转换为模拟频率Hz；
（4）根据所得结果进行分析。[www.61k.com）
三、实验内容
1．采样率和采样点数的确定
在本实验中要用到正弦波、矩形波和正弦调制波
正弦波：sin(20πt)；
矩形波：频率为50Hz、占空比为1的矩形波；
正弦波调制波：sin(20πt)×cos(100πt)
根据上述波形确定采样频率。假定所有波形的频率分辨率均为0.5Hz，确定频域采样点数。
2．信号的频谱分析
① 正弦波进行快速傅里叶变换；
② 矩形波进行快速傅里叶变换；
③ 正弦调制波进行快速傅里叶变换；
3．分析各信号的频谱与时域波形之间的关系
四、实验步骤
1．复习并理解离散傅里叶变换的定义和物理意义；
2．编写Matlab程序对信号进行频谱分析(参看例题中的程序)；
3．调试程序，排除程序中的错误；
4．分析程序运行结果，检验是否与理论一致；
5．如结果不理想，调整有关参数，得到较理想的结果。
五、实验报告要求
1．阐明实验的目的、原理和内容；
2．打印主要程序并粘贴在实验报告中；
3．打印实验结果并粘贴在实验报告中；
4．针对实验结果加以分析和总结。
六、思考题
（1）频谱的幅度有没有物理意义？如没有，怎样处理才能有物理意义？
（2）为什么所得信号的频谱均是关于中心点对称的？
（3）要让所得频谱近似为理想的冲激，该如何调整参数？
附例题
例1 试对信号x(t)＝2sin(30πt)－cos(32πt)+ sin(60πt)进行频谱分析。
解：信号中包含了3种频率：15Hz、16Hz和30Hz，最高频率为30Hz，所以采样率不
傅里叶变换对 用快速傅里叶变换对信号进行频谱分析
能低于60Hz，这里取100Hz。（www.61k.com]没有明确告诉频率分辨率，但是有两个频率仅相差1Hz，因此，频率分辨率不能低于1Hz，取0.1Hz。当然采样率越高、频率分辨率越高，则计算量就越大。程序如下：
deltf=0.1;%频率分辨率
Fs=100;%采样率
N=Fs/deltf;%采样点数
n=0:N-1;%采样点
x=2*sin(30*pi*n/Fs)-cos(32*pi*n/Fs)+sin(60*pi*n/Fs);%采样
y=fft(x);%快速傅里叶变换
ye=y.*conj(y);%计算能量
subplot(2,2,1);plot(n*Fs/N,real(y),'k');
xlabel('频率/Hz');ylabel('幅度');text(45,100,'实部');
subplot(2,2,2);plot(n*Fs/N,imag(y),'k');
xlabel('频率/Hz');ylabel('幅度');
axis([0 100 -1500 1500]);text(45,1200,'虚部');
subplot(2,2,3);plot(n*Fs/N,ye,'k');
xlabel('频率/Hz');ylabel('幅度');
axis([0 100 0 12e5]);text(45,10e5,'能量');
subplot(2,2,4);plot(n*Fs/N,ye/N^2,'k');
xlabel('频率/Hz');ylabel('幅度');
axis([0 100 0 1.5]);text(45,1.2,'功率');
二 : 快速傅里叶变换是什么快速傅里叶变化公式是什么，有什么用处？
快速傅里叶变换是什么
快速傅里叶变化公式是什么，有什么用处？
计算离散傅里叶变换的一种快速算法，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。
三 : BZOJ 3257 ZJOI2014快速傅里叶变换
题目大意：给定n个点，第i个点和第j个点之间的库仑力为(qi*qj)/(i-j)^2，定义左侧为正方向，求每个点受的合力与电荷量的比值
#include#include#include#include#include#define M 263000#define PI 3.1415926535897932384626433832795028841971using namespace std;struct Complex{long double a,b;Complex() {}Complex(long double _,long double __):a(_),b(__) {}Complex operator + (const Complex &x) const{return Complex(a+x.a,b+x.b);}Complex operator - (const Complex &x) const{return Complex(a-x.a,b-x.b);}Complex operator * (const Complex &x) const{return Complex(a*x.a-b*x.b,a*x.b+b*x.a);}void operator *= (const Complex &x){*this=(*this)*x;}}a[M],b[M],c[M];int n;long double q[M],ans[M];void FFT(Complex x[],int n,int type){static Complex temp[M];if(n==1) return ;int i;for(i=0;i>1]=x[i],temp[i+n>>1]=x[i+1];memcpy(x,temp,sizeof(Complex)*n);Complex *l=x,*r=x+(n>>1);FFT(l,n>>1,type);FFT(r,n>>1,type);Complex root(cos(type*2*PI/n),sin(type*2*PI/n)),w(1,0);for(i=0;i>1;i++,w*=root)temp[i]=l[i]+w*r[i],temp[i+(n>>1)]=l[i]-w*r[i];memcpy(x,temp,sizeof(Complex)*n);}int main(){int i,digit;double x;cin>>n;for(digit=1;digit<=n<<1;digit<<=1);for(i=0;i
扩展：bzoj3257 / 快速傅里叶变换 / 快速傅里叶变换原理
四 : 05-第五章 快速傅里叶变换(蝶形运算)
第五章 快速傅里叶变换
本章目录
直接计算DFT的问题及改进的途径 的问题及改进的途径 直接计算 按时间抽取的基2-FFT算法 时间抽取的 FFT算法 按频率抽取的基2-FFT算法 频率抽取的基 FFT算 的基2 快速傅里叶逆变换 快速傅里叶逆变换(IFFT)算法 傅里叶逆变换(IFFT)算法 Matlab实现 Matlab实现
2
5.1 引言
DFT在实际应用中很重要 DFT在实际应用中很重要: 可以计算信号的频 在实际应用中很重要: 功率谱和线性卷积等。 谱、功率谱和线性卷积等。 直接按DFT变换进行计算，当序列长度N 直接按DFT变换进行计算，当序列长度N很 变换进行计算 大时，计算量非常大，所需时间会很长。 大时，计算量非常大，所需时间会很长。 FFT并不是一种与 FFT并不是一种与DFT不同的变换，而是 并不是一种与DFT不同的变换， 不同的变换 DFT的一种快速计算的算法。 DFT的一种快速计算的算法。 的一种快速计算的算法
3
5.2 直接计算DFT的问题及改进的途径 直接计算DFT的问题及改进的途径
DFT的运算量 DFT的运算量 设复序列x(n) 长度为 点，其DFT为 长度为N点 设复序列 为
X (k ) = ∑ x( n)WNnk
n= n=0 N ?1
k=0，， ，N-1 ，，…， ，，
（1）计算一个X(k) 值的运算量 ）计算一个
复数乘法次数： 复数乘法次数： N 复数加法次数： 复数加法次数： N－1
4
5.2.1 DFT的运算量 DFT的运算量
（2）计算全部 个X(k) 值的运算量 ）计算全部N个
复数乘法次数： 复数乘法次数： N2 复数加法次数： 复数加法次数： N(N－1)
（3）对应的实数运算量 ）
X ( k ) = ∑ x( n)W
n =0
N ?1 n =0
N ?1
nk N
= ∑ [Re x( n) + j Im x( n)][Re WNnk + j Im WNnk ]
n =0
N ?1
= ∑ {[Re x( n) ? Re WNnk ? Im x(n) ? Im WNnk ]
+ j[Re x ( n) ? Im WNnk + Im x ( n) ? Re WNnk ]}
5
一次复数乘法： 次实数乘法 一次复数乘法： 4次实数乘法 一个X(k) ： 一个
＋ 2次实数加法 次实数加法
4N次实数乘法 ＋ 次实数乘法 2N+2(N-1)= 2(2N-1)次实数加法 次实数加法
所以 整个N点DFT运算共需要： 整个 点 运算共需要： 运算共需要 实数乘法次数： 实数乘法次数： 4 N2 实数加法次数： 实数加法次数： N×2(2N-1)= 2N(2N-1)
6
DFT运算量的结论 DFT运算量的结论
N点DFT的复数乘法次数举例 点 的复数乘法次数举例 N 2 4 8 16 32 N2 4 16 64 256 1028 N 64 128 256 512 1024 N2 4049 16384 65 536 262 144 1 048 576
结论： 很大时， 结论：当N很大时，其运算量很大，对实时性很强的信号 很大时 其运算量很大， 处理来说，要求计算速度快，因此需要改进DFT的计算 处理来说，要求计算速度快，因此需要改进 的计算 方法，以大大减少运算次数。 方法，以大大减少运算次数。
7
5.2.2 减少运算工作量的途径
nk 的以下特性对DFT进
行分解： 进行分解： 主要原理是利用系数 W N 的以下特性对 进行分解
（1）对称性 ）
? (WNnk )? = WN nk = WNk ( N ? n )
（2）周期性 ） （3）可约性 ） 另外， 另外，
( WN n + N ) k = WNn ( k + N ) = WNnk
mnk WmN = WNnk
/ WNnk = WNnk mm /
( k WN k + N / 2 ) = ?WN
WNN / 2 = ?1
8
5.3 按时间抽取的基2-FFT算法 按时间抽取的基2 FFT算法
算法原理 按时间抽取基-2FFT算法与直接计算 按时间抽取基-2FFT算法与直接计算 DFT运算量的比较 DFT运算量的比较 按时间抽取的FFT算法的特点 按时间抽取的FFT算法的特点 按时间抽取FFT算法的其它形式流程图 按时间抽取FFT算法的其它形式流程图
9
5.3.1 算法原理
设N＝2L，将x(n)按 n 的奇偶分为两组： ＝ 按 的奇偶分为两组：
x(2r ) = x1 (r )
x(2r + 1) = x2 (r )
r =0，1，…， ? 1 ， ， ，
N 2
则
nk X (k) = DFT[x(n)] = ∑x(n)WN n=0
N ?1 N ?1
N?1
=
n=0 n为偶数
∑ x(n)W
nk N
+
n=0 n为奇数
nk x (n)WN ∑
10
=
2 ( = ∑x(2r)WN rk + ∑x(2r +1)WN2r+1)k r =0 r =0
N ?1 2 r=0 N ?1 2 r =0
n =0 n为偶数 N ?1 2
∑ x(n)W
N ?1
nk N
+
n =0 n为奇数 N ?1 2
nk x(n)WN ∑
N ?1
rk k rk = ∑x1 (r)WN +WN ∑x2 (r)WN 2 2
k = X1 (k) +WN X 2 (k)
式中， 分别是x 式中，X1(k)和X2(k)分别是 1(n)和x2(n)的N/2的DFT。 和 分别是 和 的 的 。 另外，式中 的取值范围是 的取值范围是： ， ， ， 另外，式中k的取值范围是：0，1， …，N/2－1 。 －
11
k X 因此， 只能计算出X(k)的前一半值。 的前一半值。 因此， (k) = X1(k) +WN X2 (k) 只能计算出 的前一半值
后一半X(k) 值， N/2 ， N/2 ＋1， …，N ？ 后一半 ， ， 利用 可得到
N X1 ( + k ) = 2
rk r N WN (2 2+ k ) = WN 2
N 2 ?1 r =0
∑ x (r )W
1
r ( N 2+ k ) N 2
rk = ∑ x1 (r )WN 2 = X 1 (k ) r =0
N 2 ?1
同理可得
X2(
N + k ) = X 2 (k ) 2
12
考虑到 及前半部分X(k) 及前半部分
( WN N 2 + k ) = WNN 2 ? WNk = ?WNk
k X (k) = X1(k) +WN X2 (k)
因此可得后半部分X(k) 因此可得后半部分
X (k +
k=0，1， …，N/2－1 ， ， ， －
N N N k ) = X 1 (k + ) + WN + N 2 X 2 (k + ) 2 2 2
k = X 1 ( k ) ? WN X 2 ( k )
k=0，1， …，N/2－1 ， ， ， －
13
蝶形运算
k X (k) = X1(k) +WN X2 (k)
蝶形运算式
X ( k ) = X 1 ( k ) ? WNk X 2 ( k )
蝶形运算信 号流图符号
因此，只要求出 个 点的DFT，即X1(k)和X2(k)，再 因此，只要求出2个N/2点的 点的 ， 和 ， 经过蝶形运算就可求出全部X(k)的值，运算量大大减少。 的值， 经过蝶形运算就可求出全部 的值 运算量大大减少。
14
以8点为例第一次按奇偶分解
为例， 以N=8为例， 为例 分解为2个 点 分解为 个4点 的DFT，然后 ， 做8/2=4次蝶形 次蝶形 运算即可求出 所有8点 所有 点X(k)的 的 值。
0 WN
1 WN
2 WN
3 WN
15
蝶形运算量比较
N点D
FT的运算量 DFT的运算量
复数乘法次数： N2 复数加法次数： N(N－1)
分解一次后所需的运算量＝ N/2的DFT＋ 分解一次后所需的运算量＝2个N/2的DFT＋ N/2蝶形 N/2蝶形： 蝶形：
复数乘法次数： 2*(N/2)2+N/2=N2/2+N/2 复数加法次数： 2*(N/2)(N/2－1)+2*N/2=N2/2
因此通过一次分解后，运算工作量减少了差 因此通过一次分解后， 不多一半。 不多一半。
16
进一步按奇偶分解
由于N＝2L，因而N/2仍是偶数 ，可以进一步把每个N/2点 子序列再按其奇偶部分分解为两个N/4点的子序列。 以N/2点序列x1(r)为例 则有
X 1 (k ) =
N 2 ?1 r =0
x1 (2l ) = x3 (l ) ? N l = 0,1,…, ? 1 ? x1 (2l + 1) = x4 (l ) ? 4
rk N 2
∑ x1 (r )W
=
N 4 ?1 l =0
∑ x1 (2l )W
N 4 ?1 l =0 k N 2 4
2 lk N 2
+
N 4 ?1 l =0
( l x1 (2l + 1)WN22+1) k ∑
=
N 4 ?1 l =0
∑ x (l )W
3
lk N 4
+W
∑ x (l )W
lk N 4
k = X 3 ( k ) + WN / 2 X 4 ( k )
k=0,1,…,
N ?1 4
17
且
?N ? X 1 ? + k ? = X 3 ( k ) ? WNk / 2 X 4 ( k ) ?4 ?
k=0,1,…,
N ?1 4
由此可见，一个N/2点DFT可分解成两个 点DFT。 可分解成两个N/4点 由此可见，一个 点 可分解成两个 。 同理，也可对 进行同样的分解， 同理，也可对x2(n)进行同样的分解，求出 2(k)。 进行同样的分解 求出X 。
18
以8点为例第二次按奇偶分解
19
算法原理
对此例N=8，最后剩下的是4个N/4= 2点的 ，最后剩下的是 个 点的DFT，2点 对此例 点的 ， 点 DFT也可以由蝶形运算来完成。以X3(k)为例。 也可以由蝶形运算来完成。 为例。 也可以由蝶形运算来完成 为例
X 3 (k ) =
即
N / 4 ?1
∑
l =0
x3 (l )W
lk N /4
=
lk x3 (l )WN / 4 ∑ l =0
1
k=0, 1
0 0 X 3 (0) = x3 (0) + W20 x3 (1) = x(0) + W2 x(4) = x (0) + WN x(4)
X 3 (1) = x3 (0) + W21 x3 (1) = x(0) + W21 x(4) = x(0) ? WN0 x(4)
这说明， 可全部由蝶形运算来完成。 这说明，N=2M的DFT可全部由蝶形运算来完成。 可全部由蝶形运算来完成
20
以8点为例第三次按奇偶分解
N=8按时间抽取法 按时间抽取法FFT信号流图 按时间抽取法 信号流图
21
5.3.2 按时间抽取基2-FFT算法与直接计算DFT运算量的比较 按时间抽取基2 FFT算法与直接计算DFT运算量的比较
由按时间抽取法FFT的信号流图可知，当N=2L时，共有 L 级 蝶形运算；每级都由 N/2 个蝶形运算组成，而每个蝶形有 1 次复乘、 2 次复加，因此每级运算都需 N/2 次复乘和 N 次复加。
22
级运算总共需要： 这样 L 级运算总共需要：
N N 复数乘法： ? L = log 2 N 2 2
复数加法：N ? L = N log2 N 直接DFT算法运算量 算法运算量 直接 复数乘法： N2 复数加法： N(N－1) 直接计算DFT与FFT算法的计算量之比为 与 算法的计算量之比为M 直接计算 算法的计算量之比为 N2 2N M= = N log 2 N log 2 N 2
23
FFT算法与直
接DFT算法运算量的比较 FFT算法与直接DFT算法运算量的比较
N 2 4 8 16 32 64 N2 4 16 64 256 1028 4049
N log 2 N 2
计算量 之比M 4.0 4.0 5.4 8.0 12.8 21.4
N 128 256 512 1024 2048
N2 16 384 65 536 262 144 1 048 576 4 194 304
N log 2 N 2
计算量 之比M 36.6 64.0 113.8 204.8 372.4
1 4 12 32 80 192
448 1 024 2 304 5 120 11 264
24
5.3.3 按时间抽取的FFT算法的特点 按时间抽取的FFT算法的特点
序列的逆序排列 同址运算（原位运算） 蝶形运算两节点间的距离
r WN
的确定
25
序列的逆序排列
序列的逆序排列
由于 x(n) 被反复地按奇、偶分组，所以流图输入端的 被反复地按奇、偶分组，所以流图输入端的 排列不再是顺序的，但仍有规律可循： 排列不再是顺序的，但仍有规律可循： 因为 N=2M ， 对于任意 n（0≤n ≤N-1)，可以用 个 （ ，可以用M个 二进制码表示为： 二进制码表示为：
n(DEC) = (nM ?1nM ?2 ?n2n1n0 )(BIN)
?0 nM?1, nM?2 ,?, n2 , n1, n0 = ? ?1
n 反复按奇、偶分解时，即按二进制码的“0” “1” 分解。 反复按奇、偶分解时，即按二进制码的“ 分解。
26
倒位序的树状图（N=8） 倒位序的树状图（N=8）
27
码位的倒位序(N=8) 码位的倒位序(N=8)
自然顺序 n 0 1 2 3 4 5 6 7 二进制数 000 001 010 011 100 101 110 111 倒位序二进制数 000 100 010 110 001 101 011 111
? 倒位序顺序数 n
0 4 2 6 1 5 3 7
28
倒位序的变址处理（N=8） 倒位序的变址处理（N=8）
29
同址运算（原位运算） 同址运算（原位运算）
同址运算（原位运算）
某一列任何两个节点k 某一列任何两个节点 和j 的节点变量进行蝶形运算 得到结果为下一列k、 两节点的节点变量 两节点的节点变量， 后，得到结果为下一列 、j两节点的节点变量，而和其他 节点变量无关。这种原位运算结构可以节省存储单元， 节点变量无关。这种原位运算结构可以节省存储单元， 降低设备成本。 降低设备成本。 例 运算前
A(k) ← X1(k) X (k)
运算后
→ A(k)
A(k +
N ) ← X2 (k) 2
k WN
N → A(k + N ) X (k + ) 2 2
30
观察原位运算规律
31
蝶形运算两节点间的距离
蝶形运算两节点间的距离
以N=8为例： 为例： 为例 第一级蝶形，距离为： 第一级蝶形，距离为： 第二级蝶形，距离为： 第二级蝶形，距离为： 第三级蝶形，距离为： 第三级蝶形，距离为： 1 2 4
规律：对于共 级的蝶形而言 级的蝶形而言， 规律：对于共L级的蝶形而言，其m级蝶形运算的节 级蝶形运算的节 点间的距离为 2 m ?1
32
r WN
的确定
r WN 的确定
以N=8为例：
r m = 1时，WN = WNj / 4 = W2jm = W20 , j = 0
r m = 2时，WN = WNj / 2 = W2jm = W4j , j = 0,1
r m = 3时，WN = WNj = W2jm = W8 j , j = 0,1,2,3
N = 2 M , 第L级：
r WN = W2jL , j = 0,1,2,? ,2 L ?1 ? 1
∵
2 L = 2 M × 2 L?M = N × 2 L?M
∴W = W
r N j N ?2 L ? M
=e
?j
2π
N ?2 L ? M
?j
=e
?j
2π ? j ?2 M ? L N
=W
j ?2 M ? L N
33
5.4 按频率抽取的基2-FFT算法 按频率抽取的基2 FFT算法
算法原理
先把输入按n的顺序分成前后两半 先把输入按 的顺序分成前后两半 再把输出X(k)按k的奇偶分组 按 的奇偶分组 再把输出 设序列长度为N=2L，L为整数 设序列长度为 为整数 前半子序列x(n) 前半子序列 后半子序列 x(n +
N ) 2
0≤n≤ 0≤n≤
N ?1 2 N ?1 2
34
5.4.1 算法原理
由DFT定义得 定义得
X (k ) = ∑ x(n)WNnk
n =0 N ?1
=
N / 2 ?1
∑
n =0
x(n)W
nk N
+
n= N / 2
N / 2 ?1
∑
N ?1
x(n)WNnk
( n+ ) k N x(n + )W N 2 2 N
=
N / 2 ?1 n =0
∑ x(n)W
n=0
nk N
+
∑
n =0
=
N / 2 ?1
∑
N k? ? N nk 2 ? x(n) + x(n + )WN ? ? WN 2 ? ?
k=0，1， …，N
35
X (k ) =
N / 2 ?1
∑
n =0
N k? ? N 2 x(n) + x(n + )WN ? ?WNnk ? 2 ? ?
由于 所以
W
N 2 N
=e
?j
2π N ? N 2
= e ? jπ = ? 1
W
N k 2 N
= (?1) k
N / 2 ?1
则
X (k ) =
∑
n=0
N ? nk ? k ? x (n) + (?1) x(n + 2 )?WN ? ?
k=0，1， …，N
36
然后按k的奇偶可将 然后按 的奇偶可将X(k)分为两部分 的奇偶可将 分为两部分
? k = 2r ? ? k = 2r + 1
则式
X (k ) =
N / 2 ?1
N r=0，1， …， ? 1 ， ， ， 2
N ? nk ? x(n) + (?1) k x(n + )?WN ? 2 ? ?
∑
n=0
可转化为
X ( 2r ) =
N / 2 ?1
∑
n=0
N ? 2 nr N / 2 ?1 ? N ? ? x(n) + x( n + ) ?WN = ∑ ? x(n) + x(n + )?WNnr/ 2 ? 2 ? 2 ? n=0 ? ?
[ x ( n) ? x ( n +
N / 2 ?1 N n N )] ? WN ( 2 r +1) = {[ x (n) ? x(n + )]WNn } ? WNnr2 ∑ 2 2 n=0
X (2r + 1) =
N / 2 ?1
∑
n=0
37
令 ? N x1 (n) = x(n) + x(n + ) ?
2 ? ? n ? x2 (n) = ? x(n) ? x(n + N )?WN ? ? 2 ? ? ? ?
X (2r ) =
N / 2 ?1
N n=0，1， …， ? 1 ， ， ， 2
代入
∑
n =0
N ? nr ? x(n) + x(n + ) ? WN / 2 ? 2 ? ?
{[ x(n) ? x(n +
nr N 2
X (2r + 1) =
N / 2 ?1
∑
n =0
可得
Χ ( 2r ) =
N n nr )]WN } ? WN 2 2
N 2 ?1 n =0
∑ x (n)W
1 N 2 ?1 n =0 2
Χ (2r + 1) =
∑ x (n)W
N r=0，1， …， ? 1 ， ， ， 2
nr N 2
点的DFT，合起来正好是 点X(k)的值。 的值。 为2个N/2点的 个 点的 ，合起来正好是N点 的值
38
蝶形运算
将
N ? x1 (n) = x(n) + x(n + ) ? 2 ? ? n ? x2 (n) = ? x(n) ? x(n + N )?WN ? ? 2 ? ? ? ?
称为蝶形运算
与时间抽选基2FFT算法中的蝶形运算符号略有不同。 算法中的蝶形运算符号略有不同。 与时间抽选基 算法中的蝶形运算符号略有不同
39
例 按频率抽取(N=8) 按频率抽取(N=8)
按频率抽取， 分解为两个N/2点 的组合(N=8) 例 按频率抽取，将N点DFT分解为两个 点DFT的组合 点 分解为两个 的组合
40
与时间抽取法的推导过程一样
， 与时间抽取法的推导过程一样，由于 N=2L，N/2仍然是 仍然是 一个偶数，因而可以将每个N/2点DFT的输出再分解为偶数组 一个偶数，因而可以将每个 点 的输出再分解为偶数组 与奇数组，这就将N/2点DFT进一步分解为两个 点DFT。 进一步分解为两个N/4点 与奇数组，这就将 点 进一步分解为两个 。
N=8
41
5.4.2 频率抽取法与时间抽取法的异同
频率抽取法输入是自然顺序，输出是倒位序 的；时间抽取法正好相反。 频率抽取法的基本蝶形与时间抽取法的基本 蝶形有所不同。 频率抽取法运算量与时间抽取法相同。 频率抽取法与时间抽取法的基本蝶形是互为 转置的。
42
5.5 快速傅里叶逆变换(IFFT)算法 快速傅里叶逆变换(IFFT)算法
IDFT公式 公式
(n ) = IDFT [X (k )] = 1 x N
DFT公式 公式
N ?1 n=0
? X (k )WN nk ∑ k =0
N ?1
nk X (k ) = DFT [x(n)] = ∑ x (n)WN
比较可以看出， 比较可以看出， W Nnk IDFT多出 多出
? W N nk
N = 2M
1 ? 1? =? ? N ? 2?
M
M个1/2可分解到 级蝶形运算中。 个 可分解到 级蝶形运算中。 可分解到M级蝶形运算中
43
例 频率抽取IFFT流图(N=8) 频率抽取IFFT流图(N=8)
44
快速傅里叶逆变换另一种算法
1 N?1 ? IDFT[ X (k)] = ∑X (k)WN nk N k=0
1? 1? ? ?nk ? ? nk ? ? = ?∑X (k) WN ? = N ?∑X (k)WN ? N ?k=0 ? ?k=0 ?
∴
? 1 ? IFFT[ X (k)] = {FFT[ X (k)]} N
N?1
(
)
?
N?1
?
求 FFT X ?(k) ??? FFT[ X ?(k)] X (k) ? ?? → →
求 轭 共 除 N 以 ? → {FFT[X (k)]} ??? x(n) ?? → 求 轭 共
? ?
45
5.8 Matlab实现 Matlab实现
进行谱分析的Matlab实现 用FFT进行谱分析的 进行谱分析的 实现 进行谱分析的Matlab实现 用CZT进行谱分析的 进行谱分析的 实现 中使用的线性调频z变换函数为 在Matlab中使用的线性调频 变换函数为 ， 中使用的线性调频 变换函数为czt， 其调用格式为
>>X= czt(x, M, W, A) 其中， 是待变换的时域信号 是待变换的时域信号x(n)，其长度为 ， 其中，x是待变换的时域信号 ，其长度为N， M是变换的长度，W确定变换的步长，A确定变换 是变换的长度， 确定变换的步长 确定变换的步长， 确定变换 是变换的长度 的起点。 变成DFT。 的起点。若M= N，A= 1，则CZT变成 ， ， 变成 。
46
5.8.1 用FFT进行谱分析的Matlab实现 FFT进行谱分析的Matlab实现
例5.1 设模拟信号 x (t ) = 2 sin(4π t ) + 5 cos(8π t ) ，以 t= 0.01n (n=0: N-1) 进行取样，试用 函数对其做频谱分析。N分别 进行取样，试用fft函数对其做频谱分析 函数对其做频谱分析。 分别 为：(1) N=45；(2) N=50；(3) N=55；(2) N=60。 ； ； ； 。 程序清单如下 %计算 计算N=45的FFT并绘出其幅频曲线 计算 的 并绘出其幅频曲线 N=45;n=0:N-1;t=
本文标题：
快速傅里叶变换-用快速傅里叶变换对信号进行频谱分析 本文地址：
http://www.61k.com/1213065.html 