实验三快速傅里叶变换及其应用
一:实验目的
（1）加深对FFT的理解，熟悉matlab中的有关函数。
（2）应用FFT对典型信号进行频谱分析。
（3）了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT.
（4）应用FFT实现序列的线性卷积和相关。
二:实验原理:
在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，
它的DFT定义为：反变换为：
有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。
FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次分解，
使其成为若干小点数的组合，从而减少运算量。常用的FFT是以2为基数的，其长度。它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。
(一)在运用DFT进行频谱分析的过程中可能的产生三种误差
(1) 混叠
序列的频谱是被采样信号的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。
避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。
(2) 泄漏
实际中我们往往用截短的序列来近似很长的甚至是无限长的序列，这样可以使用较短的DFT来对信号进行频谱分析，这种截短等价于给原信号序列乘以一个矩形窗函数，也相当于在频域将信号的频谱和矩形窗函数的频谱卷积，所得的频谱是原序列频谱的扩展。
泄漏不能与混叠完全分开，因为泄漏导致频谱的扩展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减至最小。
(3) 栅栏效应
DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数，就一定意义上看，用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住，不能别我们观察到。
减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值，从而变动DFT的点数，这一方法实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或谷点暴露出来。
(二)用FFT计算线性卷积
用FFT可以实现两个序列的圆周卷积。在一定的条件下，可以使圆周卷积等于线性卷积。一般情况，设两个序列的
长度分别为N1和N2，要使圆周卷积等于线性卷积的充要条件是FFT 的长度N ≥N1＋N2 对于长度不足N 的两个序列，分别将他们补零延长到N 。
当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。有两种方法：
（1）重叠相加法。将长序列分成与短序列相仿的片段，分别用FFT 对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。
（2） 重叠保留法。这种方法在长序列分段时，段与段之间保留有互相重叠的部分，在构成总的卷积输出时只需将各段线性卷积部分直接连接起来，省掉了输出段的直接相加。
（三）用FFT 计算相关函数
两个长为N 的实离散时间序列)(n x 与)(n y 的互相关函数定义为 ：
)
()()()()()()(1
1
m y m x n y m n x m n y n x m r N n N m xy *-=
-=+=
∑∑-=-= )
(n r xy 的离散傅里叶变换为：
1
0),
()()(-≤≤=*
N k k Y k X k R xy
当)()(n y n x =时，得到)(n x 的自相关函数为：
km
N
j
N n N k xx e
k X N
m n x n x m r π22
1
1
)(1)()()(∑∑
-=-==
+=
利用FFT 求两个有限长序列线性相关的步骤（设)(n x 长1N ，)(n y 长2N ）：
（1）为了使两个有限长序列的线性相关可用其圆周相关代替而不产生混淆，选择周期L
N 2=≥121-+N N ，以便使
用FFT ，将)(n x ，)(n y 补零至长为N 。 （2）用FFT 计算1
0),(),(-≤≤N k k Y k X （3）
1
0),
()()(-≤≤=*
N k k Y k X k R xy
（4）对
)
(k R xy 作IFFT ；取后1-N 项，得11)(-≤≤+-m N m r xy ；取前N 项，得1
0)(-≤≤N m m r xy 。
三、实验内容及步骤 实验中用到的信号序列：
a)高斯（Gaussian ）序列 b)衰减正弦序列
??
?
??≤≤=--其它
0150)(2
)
(n e
n x q p n a
??
?≤≤=-其它
0150)
2sin()(n fn e n x an b π
c)三角波序列 d)反三角波序列
??
?
??≤≤-≤≤=其它07483
0)(n n
n n n x c
???
??≤≤-≤≤-=其它
0744
304)(n n n n n x d
上机实验内容：
（1）观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8，改变q 的值，使q 分别等于2，4，8，观察它们的时域和幅频特性，了解当q 取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p ，使p 分别等于8，13，14，观察参数p 变化对信号序列的时域及幅频特性的影响，观察p 等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。 参数p=8不变 主要代码如下: i=1:15; p=8;q=2;
subplot(3,2,1); x(i)=exp((-(i-p).^2)/q); stem(x);xlabel('n'); subplot(3,2,2); G=fft(x);
plot(abs(G(1:15)));xlabel('k'); 运行结果如下：
对比三张图可知：
P 不变，随着q 值的增大，时域信号幅值变化缓慢时域幅度对应变大,波形变胖，低频分量变多，频域信号频谱泄露程
参数q=8不变
代码和之前相比只要改变相应几个参数而已，故不列出；
运行结果如下:
可见，当q不变，随着p的增大，时域信号幅值不变，会在时间轴移位，对应右移,可见p决定了波形位置，
当实验中q=8,p=13时，明显出现泄漏。
(2)观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。
其中一段代码如下（其他的都是改变参数而已）：
a=0.1;f=0.0625;
x(i)=exp(-a*(i-1))*sin(2*pi*f*(i-1));
end
for i=17:100
x(i)=0;
end
n=0:15;
subplot(3,2,1);
plot(n(1:16),x(1:16));xlabel('n');
subplot(3,2,2);
G=fft(x,16);
plot(n(1:16),abs(G(1:16)));xlabel('k');
结果：
(3)观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列xc(n)和xd(n)的幅频特性，观察两者的序
列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。
主要代码如下：
for i=1:4
x(i)=i-1;
end
for i=5:8
x(i)=9-i;
end
n=0:7;
subplot(2,2,1);
plot(n,x(1:8));xlabel('n');
subplot(2,2,2);
G=fft(x,8);
plot(n(1:8),abs(G(1:8)));xlabel('k');
结果：
在xc(n)和xd(n)末尾补零，用N=32点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFT
频谱还有相同之处吗？这些变化说明了什么？
代码只需改动几处数据即可，其中一段如下所示：
for i=1:4
x(i)=i-1;
end
for i=5:8
x(i)=9-i;
end
for i=9:32
x(i)=0;
end
n=0:31;
subplot(2,1,1);
G=fft(x,32);
plot(n(1:32),abs(G(1:32)));xlabel('k');
运行结果：
变化：反三角波的低频分量增多，对信号末尾补零加长整数个周期可以对原信号达到细化频谱的作用。
(4) 一个连续信号含两个频率分量，经采样得x(n)=sin*2π*0.125n++cos*2π*(0.125+Δf)n+n=0,1……,N-1
已知N=16，Δf分别为1/16和1/64，观察其频谱；当N=128时，Δf不变，其结果有何不同，为什么？
其中一段代码如下：实现方法和之前一样
N=16;
f=1/16;
for n=1:N
x(n)=sin(2*pi*0.125*(n-1))+cos(2*pi*(0.125+f)*(n-1));
end
n=0:15;
subplot(2,2,1);
G=fft(x,16);
plot(n(1:16),abs(G(1:16)));xlabel('k');
运行结果：
（5）用FFT分别实现xa(n)（p＝8，q＝2）和 xb(n)（a＝0.1，f＝0.0625）的16点循环卷积和线性卷积。for i=1:16
x(i)=exp((-(i-1-8).^2)/2);
y(i)=exp(-0.1*(i-1))*sin(2*pi*0.0625*(i-1));
end
%for i=17:31
% x(i)=0;
% y(i)=0;
%end
n=0:30
G1=fft(x,31);
G2=fft(y,31);
z=ifft(G1.*G2,31);
subplot(2,1,1);
plot(n(1:31),z(1:31));
subplot(2,1,2);
G1=fft(x,16);
G2=fft(y,16);
z=ifft(G1.*G2,16);
n=0:15;
plot(n(1:16),z(1:16));
运行结果：
（6）产生一512点的随机序列xe(n),并用xc(n)和xe(n)作线性卷积，观察卷积前后xe(n)频谱的变化。要求将xe(n)分成8段，分别采用重叠相加法和重叠保留法。
xe=rand(1,512);
n1=0:1:3;
xc1=n1;
n2=4:7;
xc2=8-n2;
xc=[xc1,xc2];
yn=zeros(1,519);
for j=0:7
xj=xe(64*j+1:64*(j+1));
xak=fft(xj,71);
xck=fft(xc,71);
yn1=ifft(xak.*xck);
temp=zeros(1,519);
temp(64*j+1:64*j+71)=yn1;
yn=yn+temp;
end;
n=0:518;
figure(1)
subplot(2,1,1);
plot(n,yn);
xlabel('n');ylabel('y(n)');
subplot(2,1,2);
plot(n,abs(fft(yn)));
xlabel('k');ylabel('Y(k)');
axis([0,600,0,300]);
xe=rand(1,512);
k=1:7;
xe1=k-k;
xe_1=[xe1,xe];
yn_1=zeros(1,519);
for j=0:7
xj_1=xe_1(64*j+1:64*j+71);
xak_1=fft(xj_1);
xck_1=fft(xc,71);
yn1_1=ifft(xak_1.*xck_1);
temp_1=zeros(1,519);
temp_1(64*j+1:64*j+64)=yn1_1(8:71); yn_1=yn_1+temp_1;
end;
n=0:518;
figure(2)
subplot(2,1,1);
plot(n,yn_1);
xlabel('n');ylabel('y(n)'); subplot(2,1,2);
plot(n,abs(fft(yn_1))); xlabel('k');ylabel('Y(k)'); axis([0,600,0,300]);
运行结果：
（7）用FFT分别计算xa(n)(p=8,q=2)和 xb(n)（a＝0.1，f＝0.0625）的16点循环相关和线性相关，问一共有多少种结果，它们之间有何共同点？
n=0:1:15;
xan=exp(-(n-8).^2/2);
xbn=exp(-0.1*n).*sin(2*pi*0.0625*n);
k=length(xbn);
xan1=[xan zeros(1,k-1)];
xbn1=[xbn zeros(1,k-1)];
xak=fft(xan1);
xbk=fft(xbn1);
rm=real(ifft(conj(xak).*xbk));
rm1=[rm(k+1:2*k-1) rm(1:k)];
m=(-k+1):(k-1);
subplot(2,1,1);
stem(m,rm1);xlabel('n');
subplot(2,1,2);
xak=fft(xan);
xbk=fft(xbn);
rm=real(ifft(conj(xak).*xbk));
stem(n,rm);xlabel('n');
运行结果：
（8）用FFT分别计算xa(n)(p=8,q=2)和 xb(n)（a＝0.1，f＝0.0625）的自相关函数。n=0:1:15;
xan=exp(-(n-8).^2/2);
k=length(xan);
xak=fft(xan,2*k);
rm=real(ifft(conj(xak).*xak));
rm=[rm(k+2:2*k) rm(1:k)];
m=(-k+1):(k-1);
subplot(2,1,1);
stem(m,rm);xlabel('m');
运行结果：
四，思考题
（1）实验中的信号序列xc(n)和xd(n)，在单位圆上的Z变换频谱|Xc(jω)|和|Xd(jω)|会相同吗？如果不同，你能说出哪一个低频分量更多一些吗？为什么？
答：不同。反三角波的低频分量更多一些。
（2）对一个有限长序列进行DFT等价于将该序列周期延拓后进行DFS展开，因为DFS也只是取其中一个周期来计算，所以FFT在一定条件下也可以用以分析周期信号序列。如果实正弦信号sin(2πfn),f=0.1用16点FFT来做DFS运算，得到的频谱时信号本身的真实谱吗？
function x=MyIFFT_FB(y) %MyIFFT_TB:My Inverse Fast Fourier Transform Time Based %按频率抽取基2-傅里叶逆变换算法 %input: % y -- 傅里叶正变换结果,1*N的向量 %output: % x -- 逆变换结果，1*N的向量 %参考文献： % https://www.doczj.com/doc/9213568079.html,/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 -- 标量，用于临时储存 %参考文献： % https://www.doczj.com/doc/9213568079.html,/view/fea1e985b9d528ea81c779ee.html %% 输入参数个数检查
一、实验目的 1在理论学习的基础上，通过本实验加深对快速傅立叶变换的理解； 2熟悉并掌握按时间抽取FFT算法的程序； 3了解应用FFT进行信号频谱分析过程中可能出现的问题，例如混淆、泄漏、栅栏效应等，以便在实际中正确应用FFT。 二、实验内容 1仔细分析教材第六章‘时间抽取法FFT ’的算法结构，编制出相应的用FFT 进行信号分析的C语言（或MATLAB 语言）程序； 用MATLAB语言编写的FFT源程序如下： %% 输入数据f、N、T及是否补零 clc; clear; f=input('输入信号频率f：'); N=input('输入采样点数N：'); T=input('输入采样间隔T：'); C=input('信号是否补零（补零输入1，不补零输入0）：'); %补零则输入1，不补则输入0 if(C==0) t=0:T:(N-1)*T; x=sin(2*pi*f*t); b=0; e lse b=input('输入补零的个数：'); while(log2(N+b)~=fix(log2(N+b))) b=input('输入错误，请重新输入补零的个数：'); end t=0:T:(N+b-1)*T; x=sin(2*pi*f*t).*(t<=(N-1)*T); end %% fft算法的实现 A=bitrevorder(x); % 将序列按二进制倒序 N=N+b; M=log2(N); % M为蝶形算法的层数 W=exp(-j*2*pi/N); for L=1:1:M %第L层蝶形算法 B=2^L/2; % B为每层蝶形算法进行加减运算的两个数的间隔 K=N/(2^L); % K为每层蝶形算法中独立模块的个数 for k=0:1:K-1 for J=0:1:B-1
班级信工142 学号 22 姓名何岩实验组别实验日期室温报告日期成绩报告内容：(目的和要求，原理，步骤，数据，计算，小结等) 1.求信号的离散时间傅立叶变换并分析其周期性和对称性； 给定正弦信号x(t)=2*cos(2*pi*10*t),fs=100HZ,求其DTFT。 (a)代码： f=10;T=1/f;w=-10:0.2:10; t1=0:0.0001:1;t2=0:0.01:1; n1=-2;n2=8;n0=0;n=n1:0.01:n2; x5=[n>=0.01]; x1=2*cos(2*f*pi*t1); x2=2*cos(2*f*pi*t2); x3=(exp(-j).^(t2'*w)); x4=x2*x3; subplot(2,2,1);plot(t1,x1); axis([0 1 1.1*min(x2) 1.1*max(x2)]); xlabel('x(n)');ylabel('x(n)'); title('原信号x1'); xlabel('t');ylabel('x1'); subplot(2,2,3);stem(t2,x2); axis([0 1 1.1*min(x2) 1.1*max(x2)]); title('原信号采样结果x2'); xlabel('t');ylabel('x2'); subplot(2,2,2);stem(n,x5); axis([0 1 1.1*min(x5) 1.1*max(x5)]); xlabel('n');ylabel('x2'); title('采样函数x2'); subplot(2,2,4);stem(t2,x4); axis([0 1 -0.2+1.1*min(x4) 1.1*max(x4)]); xlabel('t');ylabel('x4'); title('DTFT结果x4'); （b）结果： 2.用以下两个有限长序列来验证DTFT的线性、卷积和共轭特性； (n) x1(n)=[1 2 3 4 5 6 7 8 9 10 11 12];x2(n)=R 10 (1)线性:（a）代码： w=linspace(-8,8,10000); nx1=[0:11]; nx2=[0:9]; x1=[1 2 3 4 5 6 7 8 9 10 11 12];
Matlab傅里叶变换傅里叶逆变换 %% 信号经过傅里叶变换然后进行傅里叶逆变换后信号的变化 clear all;clc; %------Author&Date------ %Author: %Date: 2013/07/31 %========================================================================== Fs=8e3; %采样率 t=0:1/Fs:1; %采样点 len=length(t); %采样长度 f1=10; %频率1 f2=100; %频率2 f3=1000; %频率3 A1=1; %幅度1 A2=0.8; %幅度2 A3=0.3; %幅度3 MaxS=A1+A2+A3; %信号幅度的最大值 signal=A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t)+A3*sin(2*pi*f3*t); X=fft(signal,len); %傅里叶变换 magX=abs(X); %信号的幅度 angX=angle(X); %信号的相位 Y=magX.*exp(1i*angX); %信号的频域表示 y=ifft(Y,len); %信号进行傅里叶逆变换 y=real(y); er=signal-y; %原始信号和还原信号的误差 subplot(311);plot(t,signal);axis([0 1 -MaxS MaxS]);xlabel('时间');ylabel('振幅');title('原始信号'); subplot(312);plot(t,y);axis([0 1 -MaxS MaxS]);xlabel('时间');ylabel('振幅');title('还原信号'); subplot(313);plot(t,er);xlabel('时间');ylabel('振幅');title('误差'); % End Script
陕西科技大学实验报告 班级信工142 学号22 姓名何岩实验组别实验日期__________ 室温_____________ 报告日期________________ 成绩报告内容：(目的和要求，原理，步骤，数据，计算，小结等) 1.求信号的离散时间傅立叶变换并分析其周期性和对称性； 给定正弦信号x(t)=2*cos(2*pi*10*t),fs=100HZ, 求其DTFT (a)代码： f=10；T=1/f；w=-10:0.2:10; t1=0:0.0001:1;t2=0:0.01:1; n1=-2; n2=8; n0=0; n=n 1:0.01: n2; x5=[ n>=0.01]; x1=2*cos(2*f*pi*t1); x2=2*cos(2*f*pi*t2); x3=(exp(-j)4(t2'*w)); x4=x2*x3; subplot(2,2,1);plot(t1,x1); axis([0 1 1.1*mi n(x2) 1.1*max(x2)]); xlabel('x( n)');ylabel('x( n)'); title('原信号x1'); xlabel('t');ylabel('x1'); subplot(2,2,3);stem(t2,x2); axis([0 1 1.1*mi n(x2) 1.1*max(x2)]); title(' 原信号采样结果x2'); xlabel('t');ylabel('x2'); subplot(2,2,2);stem( n, x5); axis([0 1 1.1*mi n(x5) 1.1*max(x5)]); xlabel(' n');ylabel('x2'); title(' 采样函数x2'); subplot(2,2,4);stem(t2,x4); axis([0 1 -0.2+1.1*mi n(x4) 1.1*max(x4)]); xlabel('t');ylabel('x4'); title('DTFT 结果x4'); (b)结果：
%傅里叶变换 clc;clear all;close all; tic Fs=128;%采样频率，频谱图的最大频率 T=1/Fs;%采样时间，原始信号的时间间隔 L=256;%原始信号的长度，即原始离散信号的点数 t=(0:L-1)*T;%原始信号的时间取值范围 x=7*cos(2*pi*15*t-pi)+3*cos(2*pi*40*t-90*pi/180)+3*cos(2*pi*30*t-90*pi/ 180); z=7*cos(2*pi*15*t-pi)+3*cos(2*pi*40*t-90*pi/180); z1=6*cos(2*pi*30*t-90*pi/180); z1(1:L/2)=0; z=z+z1; y=x;%+randn(size(t)); figure; plot(t,y) title('含噪信号') xlabel('时间（s）') hold on plot(t,z,'r--') N=2^nextpow2(L);%N为使2^N>=L的最小幂 Y=fft(y,N)/N*2; Z=fft(z,N)/N*2;%快速傅里叶变换之后每个点的幅值是直流信号以外的原始信号幅值的N/2倍（是直流信号的N倍） f=Fs/N*(0:N-1);%频谱图的频率取值范围 A=abs(Y);%幅值 A1=abs(Z); B=A; %让很小的数置零. B1=A1; A(A<10^-10)=0; % A1(A1<10^-10)=0; P=angle(Y).*A./B; P1=angle(Z).*A1./B1; P=unwrap(P,pi);%初相位值，以除去了振幅为零时的相位值 P1=unwrap(P1,pi); figure subplot(211) plot(f(1:N/2),A(1:N/2))%函数ffs返回值的数据结构具有对称性，因此只取前一半 hold on plot(f(1:N/2),A1(1:N/2),'r--') title('幅值频谱')
1.请用MATLAB编写程序，实现任意两个有限长度序列的卷积和。要求用图 形显示两个序列及卷积结果。 解：y(n)=∑x(i)h(n-i) 假设x(n)={1,2,3,4,5}; h(n)={3,6,7,2,1,6}; y(n)=x(n)*h(n) 验证：y[n]=[1,12,28,46,65,72,58,32,29,30] 【程序】 N=5 M=6 L=N+M-1 x=[1,2,3,4,5] h=[3,6,7,2,1,6] y=conv(x,h) nx=0:N-1 nh=0:M-1 ny=0:L-1 subplot(131);stem(nx,x,'*b');xlabel('n');ylabel('x(n)');grid on subplot(132);stem(nh,h,'*b');xlabel('n');ylabel('h(h)');grid on subplot(133);stem(ny,y,'*r');xlabel('n');ylabel('y(h)');grid on 【运行结果】
2.已知两个序列x[n]=cos(n*pi/2), y[n]=e j*pi*n/4x[n],请编写程序绘制 X(e jw)和Y(e jw)和幅度和相角,说明它们的频移关系。 –提示：用abs函数求幅度，用angle求相角。 【程序】 n=0:15; x=cos(n*pi/2); y=exp(j*pi*n/4).*x; X=fft(x); Y=fft(y); magX=abs(X); angX=angle(X); magY=abs(Y); angY=angle(Y); subplot(221);stem(n,magX,'*r');xlabel('频率');ylabel('幅度');grid on; subplot(222);stem(n,angX,'*b');xlabel('频率');ylabel('相位');grid on; subplot(223);stem(n,magY,'*r');xlabel('频率');ylabel('幅度');grid on; subplot(224);stem(n,angY,'*b');xlabel('频率');ylabel('相位');grid on;
MATLAB实验傅里叶分析
实验七 傅里叶变换 一、实验目的 傅里叶变换是通信系统、图像处理、数字信号处理以及物理学等领域内的一种重要的数学分析工具。通过傅里叶变换技术可以将时域上的波形分 布变换为频域上的分布，从而获得信号的频谱特性。MATLAB 提供了专门的函数fft 、ifft 、fft2(即2维快速傅里叶变换)、ifft2以及fftshift 用于实现对信号的傅里叶变换。本次实验的目的就是练习使用fft 、ifft 以及fftshift 函数，对一些简单的信号处理问题能够获取其频谱特性（包括幅频和相频特性）。 二、实验预备知识 1. 离散傅里叶变换(DFT)以及快速傅里叶变换(FFT)简介 设x (t )是给定的时域上的一个波形，则其傅里叶变换为 2()() (1)j ft X f x t e dt π∞--∞=? 显然X ( f )代表频域上的一种分布(波形)，一般来说X ( f )是复数。而傅里叶逆变换定义为： 2()() (2)j ft x t X f e df π∞-∞ =?
因此傅里叶变换将时域上的波形变换为频域上的波形，反之，傅里叶逆变换则将频域上的波形变换为时域上的波形。 由于傅里叶变换的广泛应用，人们自然希望能够使用计算机实现傅里叶变换，这就需要对傅里叶变换(即(1)式)做离散化处理，使 之符合电脑计算的特征。另外，当 把傅里叶变换应用于实验数据的分 析和处理时，由于处理的对象具有 离散性，因此也需要对傅里叶变换 进行离散化处理。而要想将傅里叶 变换离散化，首先要对时域上的波 形x (t )进行离散化处理。采用一个 时域上的采样脉冲序列: δ (t －nT ), n = 0, 1, 2, …, N －1； 可以实现上述目的，如图所示。其中N 为采样点数，T 为采样周期；f s = 1/T 是采样频率。注意采样时，采样频率f s 必须大于两倍的信号频率（实际是截止频率），才能避免混迭效应。 接下来对离散后的时域波形()()()(x t x t t n T x n T δ= -=的傅里叶变换()X f 进行离散处理。与上述做法类 似，采用频域上的δ脉冲序列： x (t δ x (t )δ t t t
快速傅里叶变换实验报告 班级： 姓名： 学号:
快速傅里叶变换 一．实验目的 1.在理论学习的基础上，通过本实验加深对快速傅立叶变换的理解； 2.熟悉并掌握按时间抽取FFT 算法的程序； 3.了解应用FFT 进行信号频谱分析过程中可能出现的问题，例如混淆、泄漏、栅栏效应等，以便在实际中正确应用FFT 。 二．实验内容 1.仔细分析教材第六章‘时间抽取法FFT ’的算法结构，编制出相应的用FFT 进行信号分析的C 语言（或MATLAB 语言）程序； 2.用FFT 程序分析正弦信号 ()sin(2)[()(*)],(0)1y t f t u t u t N T t u π=---∞<<+∞=设 分别在以下情况进行分析并讨论所得的结果： a ） 信号频率f ＝50Hz ，采样点数N=32,采样间隔T=0.000625s b ） 信号频率f ＝50Hz ，采样点数N=32,采样间隔T=0.005s c ） 信号频率f ＝50Hz ，采样点数N=32,采样间隔T=0.0046875s d ） 信号频率f ＝50Hz ，采样点数N=32,采样间隔T=0.004s e ） 信号频率 f ＝50Hz ，采样点数N=64,采样间隔T=0.000625s f ） 信号频率f ＝250Hz ，采样点数N=32,采样间隔T=0.005s g ） 将c ） 信号后补32个0，做64点FFT 三．实验要求 1.记录下实验内容中各种情况下的X (k)值，做出频谱图并深入讨论结果，说明参数的变化对信号频谱产生哪些影响。频谱只做模特性，模的最大值＝1，全部归一化；
2.打印出用C 语言（或MATLAB 语言）编写的FFT 源程序，并且在每一小段处加上详细的注释说明； 3.用C 语言（或MATLAB 语言）编写FFT 程序时，要求采用人机界面形式： N , T , f 变量均由键盘输入，补零或不补零要求设置一开关来选择。 四．实验分析 对于本实验进行快速傅里叶变换，依次需要对信号进行采样，补零（要求补零时），码位倒置，蝶形运算，归一化处理并作图。 此外，本实验要求采用人机界面形式，N,T,F 变量由键盘输入，补零或不补零设置一开关来选择。 1.采样 本实验进行FFT 运算，给出的是正弦信号，需要先对信号进行采样，得到有限 长序列()n x ， N n ...... 2,1,0= Matlab 实现： t=0:T:T*(N-1); x=sin(2*pi*f*t); 2.补零 根据实验要求确定补零与否，可以用if 语句做判断，若为1，再输入补零个数， 并将补的零放到采样得到的序列的后面组成新的序列，此时新的序列的元素个数等于原采样点个数加上补零个数，并将新的序列个数赋值给N 。 Matlab 实现： a=input('是否增加零点? 是请输入1 否请输入0\n'); if (a) ZeroNum=input('请输入增加零点的个数:\n'); else ZeroNum=0; end if (a) x=[x zeros(1, ZeroNum)];%%指令zeros(a,b)生成a 行b 列全0矩阵，在单行矩阵x 后补充0 end N=N+ZeroNum; 3.码位倒置 本实验做FFT 变换的级数为M ，N M 2log =
Matlab数字图像处理实验指导 实验目的： 通过实验，深入理解和掌握图像处理的基本技术，提高动手实践能力。 实验环境： Matlab变成 实验一图像的几何变换 实验内容：设计一个程序，能够实现图像的各种几何变换。 实验要求：读入图像，打开图像，实现图像的平移变换、比例缩放、转置变换、镜像变换、旋转变换等操作。 实验原理： 图像几何变换又称为图像空间变换，它将一幅图像中的坐标位置映射到另一幅图像中的新坐标位置。学习几何变换的关键就是要确定这种空间映射关系，以及映射过程中的变化参数。 几何变换不改变图像的像素值，只是在图像平面上进行像素的重新安排。一个几何变换需要两部分运算：首先是空间变换所需的运算，如平移、镜像和旋转等，需要用它来表示输出图像与输入图像之间的（像素）映射关系；此外，还需要使用灰度插值算法，因为按照这种变换关系进行计算，输出图像的像素可能被映射到输入图像的非整数坐标上。 设原图像f(x0,y0)经过几何变换产生的目标图像为g(x1,y1),则该空间变换（映射）关系可表示为： x1=s(x0,y0) y1=t(x0,y0) 其中，s(x0,y0)和t(x0,y0)为由f(x0,y0)到g(x1,y1)的坐标换变换函数。 一、图像平移 图像平移就是将图像中所有的点按照指定的平移量水平或者垂直移动。
二、图像镜像 镜像变换又分为水平镜像和垂直镜像。水平镜像即将图像左半部分和右半部分以图像竖直中轴线为中心轴进行对换；而竖直镜像则是将图像上半部分和下半部分以图像水平中轴线为中心轴进行对换。 三、图像转置 图像转置是将图像像素的x坐标和y坐标呼唤。图像的大小会随之改变——高度和宽度将呼唤。
实验二傅里叶分析及应用 姓名学号班级 一、实验目的 （一）掌握使用Matlab进行周期信号傅里叶级数展开和频谱分析 1、学会使用Matlab分析傅里叶级数展开，深入理解傅里叶级数的物理含义 2、学会使用Matlab分析周期信号的频谱特性 （二）掌握使用Matlab求解信号的傅里叶变换并分析傅里叶变换的性质 1、学会运用Matlab求连续时间信号的傅里叶变换 2、学会运用Matlab求连续时间信号的频谱图 3、学会运用Matlab分析连续时间信号的傅里叶变换的性质 （三）掌握使用Matlab完成信号抽样并验证抽样定理 1、学会运用MATLAB完成信号抽样以及对抽样信号的频谱进行分析 2、学会运用MATLAB改变抽样时间间隔，观察抽样后信号的频谱变化 3、学会运用MATLAB对抽样后的信号进行重建 二、实验条件 需要一台PC机和一定的matlab编程能力 三、实验内容 2、分别利用Matlab符号运算求解法和数值计算法求下图所示信号的FT，并画出其频谱图（包括幅度谱和相位谱）[注：图中时间单位为：毫秒(ms)]。
符号运算法： Ft= sym('t*(Heaviside(t+2)-Heaviside(t+1))+Heaviside(t+1)-Heaviside(t-1)+(-t)*(Heavi side(t-1)-Heaviside(t-2))'); Fw = fourier(Ft); ezplot(abs(Fw)),grid on; phase = atan(imag(Fw)/real(Fw)); ezplot(phase);grid on; title('|F|'); title('phase'); 3、试用Matlab 命令求ω ωωj 54 -j 310)F(j ++= 的傅里叶反变换，并绘出其时域信号图。
信息工程学院实验报告 课程名称： 实验项目名称：实验3 傅里叶变换及其性质 实验时间：2015/11/17 班级：通信141 姓名： 学号： 一、实 验 目 的: 学会运用MATLAB 求连续时间信号的傅里叶（Fourier ）变换；学会运用MATLAB 求连续时间信号的频谱图；学会运用MATLAB 分析连续时间信号的傅里叶变换的性质。 二、实 验 设 备 与 器 件 软件：Matlab 2008 三、实 验 原 理 3.1傅里叶变换的实现 信号()f t 的傅里叶变换定义为： ()[()]()j t F F f t f t e dt ωω∞ --∞ ==? ， 傅里叶反变换定义为：1 1()[()]()2j t f t F F f e d ωωωωπ ∞ --∞ == ? 。 信号的傅里叶变换主要包括MATLAB 符号运算和MATLAB 数值分析两种方法，下面分别加以探讨。同时，学习连续时间信号的频谱图。 3.1.1 MATLAB 符号运算求解法 MATLAB 符号数学工具箱提供了直接求解傅里叶变换与傅里叶反变换的函数fourier( )和ifourier( )。Fourier 变换的语句格式分为三种。 （1）F=fourier(f)：它是符号函数f 的Fourier 变换，默认返回是关于ω的函数。 （2）F=fourier(f,v)：它返回函数F 是关于符号对象v 的函数，而不是默认的 ω，即 ()()j v t F v f t e d t ∞ --∞ =? 。 （3）F=fourier(f,u,v)：是对关于u 的函数f 进行变换，返回函数F 是关于v 的函数，即 ()()jvu F v f t e du ∞ --∞ =?。 傅里叶反变换的语句格式也分为三种。 （1）f=ifourier(F)：它是符号函数F 的Fourier 反变换，独立变量默认为ω，默认返回是关于x 的函数。 （2）f=ifourier(F,u)：它返回函数f 是u 的函数，而不是默认的x 。 （3）f=ifourier(F,u,v)：是对关于v 的函数F 进行反变换，返回关于u 的函数f 。
实验二傅里叶分析及应用 、实验目的 （一）掌握使用Matlab 进行周期信号傅里叶级数展开和频谱分析 1、学会使用Matlab 分析傅里叶级数展开，深入理解傅里叶级数的物理含义 2、学会使用Matlab 分析周期信号的频谱特性 二）掌握使用Matlab 求解信号的傅里叶变换并分析傅里叶变换的性质 1、学会运用Matlab 求连续时间信号的傅里叶变换 2、学会运用Matlab 求连续时间信号的频谱图 3、学会运用Matlab 分析连续时间信号的傅里叶变换的性质 三）掌握使用Matlab 完成信号抽样并验证抽样定理 1、学会运用MATLAB 完成信号抽样以及对抽样信号的频谱进行分析 2、学会运用MATLAB 改变抽样时间间隔，观察抽样后信号的频谱变化 3、学会运用MATLAB 对抽样后的信号进行重建 、实验条件 Win7系统，MATLAB R2015a 三、实验内容 1、分别利用Matlab 符号运算求解法和数值计算法求下图所示信号的FT，并画出其频谱图（包括幅度谱和相位谱）［注：图中时间单位为：毫秒（ms）］。
Code: ft = sym( ' (t+2)*(heaviside(t+2)-heavisi de(t+1))+(heaviside(t+1)-heav iside(t- 1))+(2-t)*(heaviside( t-1)-heaviside(t- 2))' ); fw = simplify(fourier(ft)); subplot(2, 1, 1); ezplot(abs(fw)); grid on; title( 'amp spectrum' ); phi = atan(imag(fw) / real(fw)); subplot(2, 1, 2); ezplot(phi); grid on ; title( 'phase spectrum' ); 符号运算法 Code: dt = 0.01; t = -2: dt: 2; ft (t+2).*(uCT(t+2)- uCT(t+1))+(u CT(t+1)-uCT(t- 1))+(2-t).*(uCT (t-1)- uCT(t-2)); N = 2000; k = -N: N; w = pi * k / (N*dt); fw = dt*ft*exp(-i*t'*w); fw = abs(fw); plot(w, fw), grid on; axis([-2*pi 2*pi -1 3.5]); 数值运算法
信息工程学院实验报告 课程名称：信号与系统 实验项目名称：实验3 傅里叶变换及其性质实验时间：2013-11-29 班级： 姓名： 学号： 一、实验目的: 1、学会运用MATLAB 求连续时间信号的傅里叶（Fourier ）变换； 2、学会运用MATLAB 求连续时间信号的频谱图； 3、学会运用MATLAB 分析连续时间信号的傅里叶变换的性质。 二、实验环境: 1、硬件：在windows 7 操作环境下； 2、软件：Matlab 版本7.1 三、实验原理： 3.1傅里叶变换的实现 信号()f t 的傅里叶变换定义为： ()[()]()j t F F f t f t e dt ωω∞ --∞ == ? ， 傅里叶反变换定义为：1 1 ()[()]()2j t f t F F f e d ωωωωπ ∞ --∞ == ? 。 信号的傅里叶变换主要包括MATLAB 符号运算和MATLAB 数值分析两种方法，下面分别加以探讨。同时，学习连续时间信号的频谱图。 3.1.1 MATLAB 符号运算求解法 MATLAB 符号数学工具箱提供了直接求解傅里叶变换与傅里叶反变换的函数fourier( )和ifourier( )。Fourier 变换的语句格式分为三种。 （1）F=fourier(f)：它是符号函数f 的Fourier 变换，默认返回是关于ω的函数。 （2）F=fourier(f,v)：它返回函数F 是关于符号对象v 的函数，而不是默认的ω，即()()jvt F v f t e dt ∞ --∞ = ? 。 （3）F=fourier(f,u,v)：是对关于u 的函数f 进行变换，返回函数F 是关于v 的函数，即 ()()jvu F v f t e du ∞ --∞ =? 。 傅里叶反变换的语句格式也分为三种。 （1）f=ifourier(F)：它是符号函数F 的Fourier 反变换，独立变量默认为ω，默认返回是关于x 的函数。 （2）f=ifourier(F,u)：它返回函数f 是u 的函数，而不是默认的x 。 （3）f=ifourier(F,u,v)：是对关于v 的函数 F 进行反变换，返回关于u 的函数f 。
目录 用Matlab 对信号进行傅里叶变换 (2) Matlab 的傅里叶变换实例 (5) Matlab 方波傅立叶变换画出频谱图 (7)
用 Matlab 对信号进行傅里叶变换 1. 离散序列的傅里叶变换 DTFT(Discrete Time Fourier Transform) 代码： %原离散信号有 8 点 %原信号是 1行 8列的矩阵 %构建原始信号，为指数信号 %频域共-800 +800 的长度（本应是无穷， 高 %求 dtft 变换，采用原始定义的方法，对复指 7 subplot(311) 8 stem(n,xn); 9 title('原始信号(指数信号 )'); 10 subplot(312); 11 plot(w/pi,abs(X)); 12 title('DTFT 变换 ') 结果： 分析：可见，离散序列的 dtft 变换是周期的，这也符合 Nyquist 采样 定理的描述， 连续时间信号经周期采样之后， 所得的离散信号的频谱 是原连续信号频谱的周期延拓。 2. 离散傅里叶变换 1 N=8; 2 n=[0:1:N-1] 3 xn=0.5.^n; 4 5 w=[-800:1:800]*4*pi/800; 频分量很少，故省去) 6 X=xn*exp(-j*(n'*w)); 数分 量求和而得
与 1 中 DTFT 不一样的是， DTFT 的求和区间是整个频域，这对 N=8; % 原离散信号有 8 点 n=[0:1:N-1] %原信号是 1行 8列的矩阵 xn=0.5.^n; %构建原始信号，为指数信号 w=[-8:1:8]*4*pi/8; %频域共 -800 +800 的长度(本应是无穷， 高频分量很少， 故省去) X=xn*exp(-j*(n'*w)); %求 dtft 变换，采用原始定义的方法，对复指数分量求和而得 subplot(311) stem(n,xn); w1=[-4:1:4]*4*pi/4; X1=xn*exp(-j*(n'*w1)); title(' 原始信号 (指数信号 )'); subplot(312); stem(w/pi,abs(X)); title(' 原信号的 16 点 DFT 变换 ') subplot(313) stem(w1/pi,abs(X1)); title(' 原信号的 8 点 DFT 变换 ') 计算机的计算来说是不可以实现的， DFT 就是序列的有限傅里叶变换。 实际上， 1 中代码也只是对频域的 -800 +800 中间的 1601 结果图： 分析： DFT 只是 DTFT 的现实版本，因为 DTFT 要求求和区间无穷， 而 DFT 只在有限点内求和。 3. 快速傅里叶变换 FFT ( Fast Fourier Transform ) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
一、傅立叶变化的原理； （1）原理 正交级数的展开是其理论基础！将一个在时域收敛的函数展开成一系列不同频率谐波的叠加，从而达到解决周期函数问题的目的。在此基础上进行推广，从而可以对一个非周期函数进行时频变换。 从分析的角度看，他是用简单的函数去逼近（或代替）复杂函数，从几何的角度看，它是以一族正交函数为基向量，将函数空间进行正交分解，相应的系数即为坐标。从变幻的角度的看，他建立了周期函数与序列之间的对应关系；而从物理意义上看，他将信号分解为一些列的简谐波的复合，从而建立了频谱理论。 当然Fourier积分建立在傅氏积分基础上，一个函数除了要满足狄氏条件外， 一般来说还要在积分域上绝对可积，才有古典意义下的傅氏变换。引入衰减因子e^(-st)，从而有了Laplace变换。（好像走远了）。 （2）计算方法 连续傅里叶变换将平方可积的函数f（t）表示成复指数函数的积分或级数形式。 这是将频率域的函数F(ω)表示为时间域的函数f（t）的积分形式。 为 连续傅里叶变换的逆变换 (inverse Fourier transform) 即将时间域的函数f（t）表示为频率域的函数F(ω)的积分。 一般可称函数f（t）为原函数，而称函数F(ω)为傅里叶变换的像函数，原函数和像函数构成一个傅里叶变换对（transform pair）。 二、傅立叶变换的应用； DFT在诸多多领域中有着重要应用，下面仅是颉取的几个例子。需要指出 的是，所有DFT的实际应用都依赖于计算离散傅里叶变换及其逆变换的快速算
法，即快速傅里叶变换（快速傅里叶变换（即FFT ）是计算离散傅里叶变换及其逆变换的快速算法。）。（1）、频谱分析DFT 是连续傅里叶变换的近似。因此可以对连续信号x(t)均匀采样并截断以得到有限长的离散序列，对这一序列作离散傅里叶变换，可以分析连续信号x(t)频谱的性质。前面还提到DFT 应用于频谱分析需要注意的两个问题：即采样可能导致信号混叠和截断信号引起的频谱泄漏。可以通过选择适当的采样频率（见奈奎斯特频率）消减混叠。选择适当的序列长度并加窗可以抑制频谱泄漏。（2）、数据压缩由于人类感官的分辨能力存在极限，因此很多有损压缩算法利用这一点将语音、音频、图像、视频等信号的高频部分除去。高频信号对应于信号的细节，滤除高频信号可以在人类感官可以接受的范围内获得很高的压缩比。这一去除高频分量的处理就是通过离散傅里叶变换完成的。将时域或空域的信号转换到频域，仅储存或传输较低频率上的系数，在解压缩端采用逆变换即可重建信号。（3）、OFDM OFDM （正交频分复用）在宽带无线通信中有重要的应用。这种技术将带宽为N 个等间隔的子载波，可以证明这些子载波相互正交。尤其重要的是，OFDM 调制可以由IDFT 实现，而解调可以由DFT 实现。OFDM 还利用DFT 的移位性质，在每个帧头部加上循环前缀（Cyclic Prefix ），使得只要信道延时小于循环前缀的长度，就能消除信道延时对传输的影响。三、傅里叶变换的本质； 傅里叶变换的公式为dt e t f F t j ?+∞∞--=ωω)()(可以把傅里叶变换也成另外一种形式： t j e t f F ωπ ω),(21)(=可以看出，傅里叶变换的本质是内积，三角函数是完备的正交函数集，不同频率的三 角函数的之间的内积为0，只有频率相等的三角函数做内积时，才不为0。)(2,21)(2121Ω-Ω==?Ω-ΩΩΩπδdt e e e t j t j t j
快速傅里叶变换实验报告
————————————————————————————————作者：————————————————————————————————日期: ?
快速傅里叶变换实验报告 机械34班 刘攀 2０1３0105５8 一、 基本信号（函数)的FF Ｔ变换 1. 000()sin()sin 2cos36x t t t t π ωωω=+++ 1) 采样频率08s f f =,截断长度N ＝１６； 取02ωπ=ｒad/ｓ，则0f =1Hz ，s f =8Ｈz ，频率分辨率 f ?=s f f N ?= =0.5Hz 。 最高频率c f ＝30f =３Hz ，s f ＞2c f ,故满足采样定理,不会发生混叠现象。 截断长度02T T =,整周期截取,不会发生栅栏效应。理论上有一定的泄漏,但在整周期 截取的情况下，旁瓣上的采样都约为 0，泄漏现象没有体现出来。 频谱图如下:
幅值误差0A ?=,相位误差0??=。 2) 采样频率08s f f =，截断长度N=3２; 取02ωπ=rad/s ，则0f =１Hz,s f ＝８Hz ，频率分辨率f ?=s f f N ?==０.2５Hz 。 最高频率c f =30f =3H ｚ,s f >２c f ,故满足采样定理,不会发生混叠现象。 截断长度04T T =,整周期截取,不会发生栅栏效应。理论上有一定的泄漏，但在整周期 截取的情况下,旁瓣上的采样都约为 0，泄漏现象没有体现出来。 频谱图如下：
幅值误差0A ?=，相位误差0??=。 2. 00()sin()sin116x t t t π ωω=++ 1) 采样频率08s f f =,截断长度N=16； 取02ωπ=ra ｄ/s,则0f =１Hz ，s f =８Hz,频率分辨率f ?=s f f N ?==0．5H ｚ。 最高频率c f =110f =11H ｚ,s f <２c f ，故不满足采样定理，会发生混叠现象。 截断长度02T T =,整周期截取，不会发生栅栏效应。理论上有一定的泄漏,但在整周期 截取的情况下，旁瓣上的采样都约为 0，泄漏现象没有体现出来。 频谱图：
1 利用FFT 计算连续时间信号的傅里叶变换 设()x t 是连续时间信号，并假设0t <时()0x t =，则其傅里叶变换由下式给出 0()()i t X x t e dt ωω∞ -=? 令Γ是一个固定的正实数，N 是一个固定的正整数。当,0,1,2,,1k k N ω=Γ=-L 时，利用FFT 算法可计算()X ω。 已知一个固定的时间间隔T ，选择T 足够小，使得每一个T 秒的间隔(1)nT t n T ≤<+内，()x t 的变化很小，则式中积分可近似为 (1)0 ()()()n T iwt nT n X e dt x nT ω∞+-==∑? (1)01[ ]()i t t n T t nT n e x nT i ωω ∞-=+==-=∑ 0 1()i T i nT n e e x nT i ωωω-∞-=-=∑ （27） 假设N 足够大，对于所有n N ≥的整数，幅值()x nT 很小，则式（27）变为 1 01()()i T N i nT n e X e x nT i ωωωω---=-=∑ （28） 当2/k NT ωπ=时，式（28）两边的值为 2/2/12/0211()()[]2/2/i k N i k N N i nk N n k e e X e x nT X k NT i k NT i k NT ππππππ----=--==∑ （29） 其中[]X k 代表抽样信号[]()x n x nT =的N 点DFT 。最后令2/NT πΓ=，则上式变为 2/1()[]0,1,2,,12/i k N e X k X k k N i k NT ππ--Γ==-L （30） 首先用FFT 算法求出[]X k ，然后可用上式求出0,1,2,,1k N =-L 时的()X k Γ。 应该强调的是，式（28）只是一个近似表示，计算得到的()X ω只是一个近似值。通过取更小的抽样间隔T ，或者增加点数N ，可以得到更精确的值。如果B ω>时，幅度谱()X ω很小，对应于奈奎斯特抽样频率2s B ω=，抽样间隔T 选择/B π比较合适。如果已知信号只在时间区间10t t ≤≤内存在，可以通过对1nT t >时的抽样信号[]()x n x nT =补零，使N 足够大。 例1 利用FFT 计算傅里叶变换
实验三周期信号的傅里叶级数分析及MATLAB实现 一、实验目的： 1.利用MATLAB实现周期信号的分解与合成，并图示仿真结果； 2.用MATLAB实现周期信号的频谱，画图观察和分析周期信号的频谱； 3.通过MATLAB对周期信号频谱的仿真，进一步加深对周期信号频谱理论知识的理解。 二、实验内容 9.1（a）：程序： display('Please input the value of m(傅里叶级数展开项数)'); m=input('m='); t=-3*pi:0.01:3*pi; n=round(length(t)/4); f=cos(t).*(heaviside(t+2.5*pi)-heaviside(t+1.5*pi)+heaviside(t+0.5*pi)-heaviside(t-0.5 *pi)+heaviside(t-1.5*pi)-heaviside(t-2.5*pi)); y=zeros(m+1,max(size(t))); y(m+1,:)=f'; figure(1); plot(t/pi,y(m+1,:)); grid; axis([-3 3 -1 1.5]); title('半波余弦'); xlabel('单位:pi','Fontsize',8); x=zeros(size(t)); kk='1'; syms tx n T=2*pi; fx=sym('cos(tx)'); Nn=30; An=zeros(m+1,1); Bn=zeros(m+1,1); a0=2*int(fx,tx,-T/4,T/4)/T an=2*int(fx*cos(2*pi*(n+eps/2)*tx/T),tx,-T/4,T/4)/T bn=2*int(fx*sin(2*pi*(n+eps/2)*tx/T),tx,-T/4,T/4)/T An(1)=double(vpa(a0,Nn)); An(2)=0.5; for k=2:m An(k+1)=double(vpa(subs(an,n,k),Nn)); Bn(k+1)=double(vpa(subs(bn,n,k),Nn));