实验七 傅里叶变换
一、实验目的
傅里叶变换是通信系统、图像处理、数字信号处理以及物理学等领域内的一种重要的数学分析工具。通过傅里叶变换技术可以将时域上的波形分 布变换为频域上的分布,从而获得信号的频谱特性。MA TLAB 提供了专门的函数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 nT x nT δ=-=的傅里叶变换()X f 进行离散处理。与上述做法类似,采用频域上的δ脉冲序列:
δ ( f -n/T 0), n = 0, 1, 2, …, N -1;T 0= NT 为总采样时间
可以实现傅里叶变换()X f 的离散化,如下图示。不难看出,离散后的傅里叶变换其频率间隔(频率轴上离散点的间隔,即频域分辨率)
x (t )
δ 脉冲序列
x (t )δ (t -nT )
t
t
t
011
(3)s f f T NT N
?===
因此要增加分辨率须增加采样点数目N 。频域上每个离散点对应的频率为:
0; 0,1,2,...,-1 (4)s n f n n
f n n N T NT N
=
===显然n = 0的点对应于直流成分。
经过以上离散化处理之后,连续积分的傅里叶变换(1)式转变为如下离散形式:
1
2/0
()(), 0,1,2,..., 1 (5)
N j nk N n k k X f x t e n N π--===-∑其中t k = kT (k =0,1,2,…,N-1)代表采样点时刻。X ( f n )一般是复数,因此离散傅里叶变换(DFT)后变成一个N 点(采样点数)的复数序列。X ( f n )绝对值代表振幅,其幅角代表相位,因此由(5)式可以给出DFT 的振幅频谱和相位频谱。(5)式通常又简写成如下形式:
1
0()(), 0,1,2,..., 1 (6)N nk
N k X n x k W n N -===-∑
其中 2/ j N N W e π-=,x 是采样点数据,它是一个N 个点的向量,DFT 的结果X 是N 个点的复数向量。(5)式或(6)式就是对傅里叶变换进行数值计算的基础。
一般采样点数N 越大,DFT 的结果越接近真实的情况,但是当N 较大时,(6)式的计算量很大,因为使用计算机求解(6)式时,总共要执行N 2次复数乘法和N×(N-1)次复数加法。所以直接用DFT 算法(即(5)式)进行谱分析和信号的实时处理是不切实际的。为了减轻计算的压力,人们提出了一种所谓快速傅里叶变换(FFT )的思想:
取N =2m ,首先将N 个点的采样数据011[,,...,]N x x x x -=分成两个N /2点的序列:
1022[,,...,]N x x x x -= (偶数序列) 2131[,,...,]N x x x x -= (奇数序列)
这样处理的好处是可以把(6)式分解为两个N/2点的DFT ,使计算量降下来。接下来再将N/2点的序列x 1仿照上述做法进一步分裂成2个N/4点的序列x 3和x 4,另一序列x 2亦做如此处理,分裂成2个N/4点的序列x 5和x 6。这样两个N/2点的序列分成了更短的4个N/4点的序列,依次类推,最后的结果是将一个N 点的序列x 裂成了N 个点的单点序列:x 0, x 1, x 2, …, x N-1。这样做可以将DFT 的运算效率提高1-2个数量级,为数字信号处理技术应用于各种信号的实时处理创造了条件,从而推动数字处理技术的发展。由此可见FFT 的思想实质是
f
f
f
()X f
δ (f -n/T 0)
X ( f
)
混迭
不断地把长序列的DFT计算分解成若干短序列的DFT,并利用旋转因子(即W N )的周期性和对称性来减少DFT的运算次数。所以FFT就是DFT的快速算法。
有关FFT算法的详细介绍和理论推导参见有关的书籍,这里不做进一步介绍。
2. FFT的MATLAB实现
为了实现快速傅里叶变换,MA TLAB提供了fft、ifft、fft2、ifft2以及fftshift函数,分别用于一维和二维离散傅里叶变换(DFT)及其逆变换。借助这些函数可以完成很多信号处理任务。考虑到信号处理包含的领域很广泛,这里只介绍一维傅里叶变换及其逆变换函数。
(1) fft函数
该函数使用了快速算法来实现时域信号的离散傅里叶变换。常用的格式:
Y= fft (x)
Y= fft (x, m)
Y返回值(复数),返回m点的DFT序列,即(6)式左边的X;
m计算时使用的数据点数(样本数);
x时域信号x(t)在采样点t k处的值,即(6)式右边的x;若实际采样点数目为N(m和N都须是2的幂次),则x为N个元素(即长度N)的向量;若向量x的长度小于m,那么计算时将自动在x序列的后面补0;若x的长度大于m,则x自动截断,使之长度为m。对信号进行频谱分析时,数据样本应有足够的长度,一般FFT程序中所用数据点数(m)最好与原信号含有的数据点数(即输入的样本数N)相同,这样的频谱图具有较高的质量,可减小因补零或截断而产生的影响。
两点说明:
①关于FFT振幅频谱和相位频谱的计算
由于傅里叶变换的结果一般是复数,所以
●对fft的结果取绝对值abs()可以得到振幅,即
Amplitude = abs(Y)
需要注意的是这样得到的幅值实际并非真正的信号振幅,因其值与FFT使用的数据点数N 有关,但不影响分析结果,在IFFT(逆变换)时已经做了处理。要得到真实的振幅值的大小,只要将上述结果除以N/2即可。
●对fft的结果使用函数angle()可以得到相位的结果。但是使用angle函数计算复数的相角时,系统规定一、二象限的角为0~π;三、四象限的角为-π~0。因此若一个角度本来应该从0变到2π,但计算得到的结果却是0~π,再由-π~0,在π处发生跳变,跳变幅度为2π,这就叫相位的卷绕。这种相位的卷绕会使得相频图不连续,呈现锯齿状,为了平滑相频图,通常要再使用unwrap()函数进行相位的解卷绕。因此FFT的相位频谱图应该如下实现Phase = unwrap(angle(Y))
②FFT的振幅频谱具有对称性
如下图所示。
因此用FFT 对信号做谱分析,只需考察0~Nyquist 频率范围内(共N/2+1个频率点)的幅频特性。 (2) fftshift 函数
其作用是将零频点移到频谱的中间(即Nyquist 频率处),使用格式: Y=fftshift(X)
X 是向量,该命令将零频点移动到频谱X 的中间,并交换频谱X 的左右两半。将零频点放到频谱的中间对于观察傅立叶变换是有用的。
例1:对时域信号()0.5sin(215)2sin(240)x t t t ππ=?+?进行频谱分析。
fs=100; % 采样频率>2倍的信号频率 N=256; % 采样点数目(=2的幂次) n=0:N-1; % 构造采样点序列
t=n/fs; % 得到采样时间序列,t=nT=n/fs
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); % 产生时域信号的样本值,向量 Y=fft(x,N); % N 点的DFT 计算 mag=abs(Y); % FFT 的振幅
phase=unwrap(angle(Y)); % FFT 的相位 % 1. 以下绘制物理频谱图(即正频部分)
fn=(0:N/2)*fs/N; % 频率轴上的离散频率点,起始于0频(对应直流成分),终
%于Nyquist 频率fs/2,共N/2+1个频率点
subplot(2,2,1)% 将图形窗口分割为2×2的子窗口,并指定第1个子窗口为绘图区 plot(fn,mag(1:N/2+1)) % 取出前N/2+1个振幅作图,即正频率分量
…
…
f n : f 轴
1T 0
2T 0/22
s N f T =
1N T -
s N
f T = 对称轴
(Nyquist 频率)
负频部分→
← 正频部分 FFT :
注意:频率轴上最后一个点的频率等于采样频率, 实际上不存在。 0, (0,1,2,...,1)n s n n n f f n N T NT N
=
===- 为频率轴上的频率点。 特别: 221
2
N s s N f f f N =
= 为Nyquist 频率,整个FFT 频谱关于此频率对称, 右边实际上是负频
xlabel('频率/Hz');ylabel('振幅');
title('图1:物理(正频)幅频图');grid on % 加网格线
% 2. 以下绘制全频率的幅频图
fn1=(0:N-1)*fs/N;
subplot(2,2,2) %指定第2个子窗口为绘图区
plot(fn1,mag);
xlabel('频率/Hz');ylabel('振幅');
title('图2:全频率的幅频图'); grid on
%3. 以下绘制正频部分的相频图
subplot(2,2,3) %指定第3个子窗口为绘图区
plot(fn,phase(1:N/2+1));
xlabel('频率/Hz');ylabel('相位');title('图3:相频图');grid
%4. 以下移动零频点
Y1=fftshift(Y); % fftshift移动频率零点,并将Y的左右两部分交换
mag1=abs(Y1); % 重新计算振幅
fn2=fn1-fs/2; % 零点移动到fs/2处,故需重新标记频率轴
subplot(2,2,4); %指定第4个子窗口为绘图区,最终4幅图绘制在一张图上了plot(fn2,mag1);
xlabel('频率/Hz');ylabel('振幅');
title('图4: fftshift后的幅频图');grid
运行结果如下:
图
说明:1是物理谱图(正频部分),从中看到,该信号包含两个频率15Hz 和40Hz 。由于使用的采样频率fs=100Hz ,所以Nyquist 频率为50Hz ,在图2中明显能看到整个频谱图关于Nyquist 频率对称,不过Nyquist 频率右边的谱图实际上是负频部分,没有意义。图4是fftshift 之后的幅频图,由于它是图2结果的左右交换,因此图2右边变成了负频。另外,图中的振幅不是真实的信号振幅,从信号x(t)的表达式我们知道15Hz 和40Hz 这两种频率成分的振幅分别是0.5和2。要得到真实的振幅,只需要将程序中的mag 除以N/2即可。 (3) ifft 函数
执行离散傅里叶变换的逆变换,格式 x = ifft (Y) 或者 x = ifft (Y , m)
Y 是FFT 的输出结果,返回值x 是时域上的结果,m 仍然是计算使用的数据点数。在上例中若程序末尾使用: xx=ifft(Y , N),则得到采样时刻点上,信号x(t)的样本值。
三、实验内容及要求
实验项目: 给定采样频率51.2Hz 及采样点数N=512,计算矩形函数 1 01()0 1
t x t t <≤?=?
>?
的振幅频谱,并与理论计算结果对比。
A. 显然该信号x (t )是无限长的非周期信号,因此做FFT 计算时必须先将信号x(t)截断为有限长度。令采样频率为f s ,采样点数目N ,则截断长度是:
T 0 = N / f s (即总的采样时间)
因此截断长度和采样点数目N 成正比。对于无限长的非周期信号,截断长度应尽可能的大,以接近实际信号,避免结果失真;如果是周期信号,则要求截断长度为信号周期的整数倍,以免出现频谱的“泄漏”。
若给定采样点数N=512,则时间采样序列可用向量t 表示:t=(0:N-1)/fs ,矩形函数x(t)的样本值可以使用MA TLAB 提供的符号函数sign (请使用help sign 命令查询sign 函数的定义)来表示:x=0.5-0.5*sign(t-1),然后使用fft 命令即可获得DFT 计算结果。 B. 根据(1)式不难算出上述信号x(t)的傅里叶变换的理论结果(精确值):
()[sin 2(cos21)] (7)X f f j f f πππ=+-
其振幅为:
() (8)X f f π=
然后再根据(8)式即可绘出x(t)的振幅频谱,将其与fft 计算结果比较:改变N 以改变截断长度观察FFT 结果与(8)式结果的差异。 要求: ⑴ 认真阅读例1提供的程序;
⑵ 参考例1中的程序,编写出本实验项目的fft 计算程序,绘制出零点移动到
Nyquist频率处的振幅频谱1;对图形进行标注(如参考图所示)2;使用(8)式绘制矩形函数x(t)精确的振幅频谱图3;使用subplot命令将2幅图在一张图上显示4;频率区间[-5,5]观察比较得到的两个振幅频谱5;改变采样数N,再比较两幅图的差异,分析采样点数对fft 计算结果精确度的影响6。
附实验结果参考图:
主瓣
边瓣
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.sodocs.net/doc/1211580736.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.sodocs.net/doc/1211580736.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
傅里叶分析之掐死教程(完整版)更新于2014.06.06 Heinrich · 6 个月前 作者:韩昊知乎:Heinrich 微博:@花生油工人知乎专栏:与时间无关的故事 谨以此文献给大连海事大学的吴楠老师,柳晓鸣老师,王新年老师以及张晶泊老师。 转载的同学请保留上面这句话,谢谢。如果还能保留文章来源就更感激不尽了。 我保证这篇文章和你以前看过的所有文章都不同,这是12年还在果壳的时候写的,但是当时没有来得及写完就出国了……于是拖了两年,嗯,我是拖延症患者…… 这篇文章的核心思想就是: 要让读者在不看任何数学公式的情况下理解傅里叶分析。 傅里叶分析不仅仅是一个数学工具,更是一种可以彻底颠覆一个人以前世界观的思维模式。但不幸的是,傅里叶分析的公式看起来太复杂了,所以很多大一新生
上来就懵圈并从此对它深恶痛绝。老实说,这么有意思的东西居然成了大学里的杀手课程,不得不归咎于编教材的人实在是太严肃了。(您把教材写得好玩一点会死吗?会死吗?)所以我一直想写一个有意思的文章来解释傅里叶分析,有可能的话高中生都能看懂的那种。所以,不管读到这里的您从事何种工作,我保证您都能看懂,并且一定将体会到通过傅里叶分析看到世界另一个样子时的快感。至于对于已经有一定基础的朋友,也希望不要看到会的地方就急忙往后翻,仔细读一定会有新的发现。 ——————————————以上是定场诗—————————————— 下面进入正题: 抱歉,还是要啰嗦一句:其实学习本来就不是易事,我写这篇文章的初衷也是希望大家学习起来更加轻松,充满乐趣。但是千万!千万不要把这篇文章收藏起来,或是存下地址,心里想着:以后有时间再看。这样的例子太多了,也许几年后你都没有再打开这个页面。无论如何,耐下心,读下去。这篇文章要比读课本要轻松、开心得多…… p.s.本文无论是cos还是sin,都统一用“正弦波”(Sine Wave)一词来代表简谐波。 一、什么是频域 从我们出生,我们看到的世界都以时间贯穿,股票的走势、人的身高、汽车的轨迹都会随着时间发生改变。这种以时间作为参照来观察动态世界的方法我们称其为时域分析。而我们也想当然的认为,世间万物都在随着时间不停的改变,并且永远不会静止下来。但如果我告诉你,用另一种方法来观察世界的话,你会发现世界是永恒不变的,你会不会觉得我疯了?我没有疯,这个静止的世界就叫做频域。 先举一个公式上并非很恰当,但意义上再贴切不过的例子: 在你的理解中,一段音乐是什么呢?
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
实验报告 课程名称信号与线性系统分析 实验名称周期信号的傅里叶级数和频谱分析实验类型验证(验证、综合、设计、创新) 3日实验四、周期信号的傅里叶级数和频谱分析1实验目的 1)学会利用MATLAB分析傅里叶级数展开,并理解傅里叶级数的物理含义; 2)学会利用MATLAB分析周期信号的频谱特性。 2实验原理及实例分析 周期信号可以再函数的区间里展成在完备正交信号空间中的无穷级数。如果完备的正交函数集是三角函数集或指数函数集,那么,周期信号所展开的无穷级数就分别成为“三角型傅里叶级数”或“指数型傅里叶级数”,统称为傅里叶级数。
2.1周期信号的傅里叶级数 (基本原理请参阅教材第四章的4.1节和4.2节。) 例1:周期方波信号)(t f 如图1所示,试求出该信号的傅里叶级数,利用MATLAB 编程实现其各次谐波的叠加,并验证Gibbs 现象。 图1 周期方波信号)(t f 的波形图 解:从理论分析可知,周期方波信号)(t f 的傅里叶级数展开式为 )9sin 9 17sin 715sin 513sin 31(sin 4)(00000 +++++=t t t t t t f ωωωωωπ 其中,ππω220== T 。则可分别求出1、3、5、9、19、39、79、159项傅里叶级数求和的结果,其MATLAB 程序如下,产生的图形如图2所示。 close all;clear all; clc t = -2:0.0001:2; omega = 2 * pi; y = square(2 * pi * t,50); n_max = [1 3 5 9 19 39 79 159]; N = length(n_max); for k = 1:N fk = zeros(1,length(t)); for n = 1:2:n_max(k) bn = 4 / (pi * n); fk = fk + bn * sin(n * omega * t); end figure;plot(t,y,t,fk,'Linewidth',2); xlabel('t(sec)');ylabel('部分和的波形'); f(t) t(sec)
%傅里叶变换 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
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 ++= 的傅里叶反变换,并绘出其时域信号图。
实验二傅里叶分析及应用 、实验目的 (一)掌握使用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]); 数值运算法
目录 用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
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));
一、实验目的 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 。 function [y,z]=disbutterfly(x) N=length(x); n=0:N/2-1; w=exp(-2*1i*pi/N).^n; x1=x(n+1); x2=x(n+1+N/2); y=x1+x2; z=(x1-x2).*w; 2.定义函数rader ,纠正输出序列的输出顺序。 function y=rader(x,N) n=[0:N-1]; bn=dec2bin(n); rbn=fliplr(bn); rn=bin2dec(rbn); y=x(rn+1); 3.定义函数myfft ,程序中套了两个循环。 function X=myfft(x) N=length(x); h=log2(N); %h=3 for i=1:h %第一次i=1;第二次i=2 s=[]; for j=1:2^(i-1);%i=1时,j=1;i=2时,j=1:2 M=2^(h-i+1);%M:M=8;M=4 xj=x([1:M]+(j-1)*M);%xj=x([1:8]+(1-1)*8)=x(1)+x(2)...+x(8); %j=1:xj=x([1:4]);j=2:xj=x([1:4]+4) [y,z]=disbutterfly(xj); s=[s,y,z]; end x=s;
实验七 傅里叶变换 一、实验目的 傅里叶变换是通信系统、图像处理、数字信号处理以及物理学等领域内的一种重要的数学分析工具。通过傅里叶变换技术可以将时域上的波形分 布变换为频域上的分布,从而获得信号的频谱特性。MA TLAB 提供了专门的函数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 nT x nT δ=-=的傅里叶变换()X f 进行离散处理。与上述做法类似,采用频域上的δ脉冲序列: δ ( f -n/T 0), n = 0, 1, 2, …, N -1;T 0= NT 为总采样时间 可以实现傅里叶变换()X f 的离散化,如下图示。不难看出,离散后的傅里叶变换其频率间隔(频率轴上离散点的间隔,即频域分辨率) x (t ) δ 脉冲序列 x (t )δ (t -nT ) t t t
实验二傅里叶分析及应用 -、实验目的 (一)掌握使用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,并画出其频
谱图(包括幅度谱和相位谱)
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)); gridon ; 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]); 数值运算法amp spectrum -6-4 -2 0 2 4 6 w x10 phase spectrurri -6 -4 -2 0 2 4 6 w
fs=51.2; N=1024; n=0:N-1; t=n/fs; x=0.5-0.5*sign(t-1); Y=fft(x,N); mag=abs(Y); Y1=fftshift(Y); mag1=abs(Y1); fn2=(-N/10.24:N/10.24)*fs/N; subplot(2,1,1) plot(fn2,mag1((N/2-N/10.24+1):(N/2+N/10.24+1))); set(gca,'XTick',(-5:0.5:5)); set(gca,'YTick',(0:10:60)); xlabel('频率/Hz');ylabel('振幅'); title('图1:矩形函数的FFT结果,N=512,fs=51.2Hz'); grid on; f=linspace(-5,5,1000); y=sqrt(2-2*cos(2*pi*f))./abs((2*pi*f)); subplot(2,1,2) plot(f,y); set(gca,'XTick',(-5:0.5:5)); set(gca,'YTick',(0:0.2:1)); xlabel('频率/Hz');ylabel('振幅'); title('图2:矩形函数傅里叶变换的理论结果'); grid on
-5-4.5-4-3.5-3-2.5-2-1.5-1-0.500.51 1.52 2.53 3.54 4.55010203040 50 60 频率/Hz 振幅 图1:矩形函数的FFT 结果,N=512,fs=51.2Hz -5-4.5-4-3.5-3-2.5-2-1.5-1-0.500.51 1.52 2.53 3.54 4.55 00.2 0.4 0.60.8 1频率/Hz 振幅图2:矩形函数傅里叶变换的理论结果