运用Matlab 实现光学中的几个傅立叶变换
摘要：光学中的傅立叶变换具有难度高，抽象性大的特点，而Matlab 却具有强大的信号处理功能，结合光学中傅立叶变换、傅立叶级数、卷积定理的内容，通过Matlab 程序来体现 光学中几个傅里叶调动的直观演示。 关键词：傅立叶变换 Matlab 程序 傅立叶光学 1. 引言：
傅立叶变换的原理由正交级数的展开来完成，是将一个在时域收敛的函数展开成一系列不同频率谐波的叠加，从而达到解决周期性函数问题的目的。在此根底上实行变动，对非周期函数进行时频变更。
跟着科学技术的不停成长与发展，信号处理在人们生产生活中应允了越来越普遍的应用，其中，采用通信频域的方法比起经典的方法体现出越来越多的优点，当今，光学中傅立叶变换和信号处理成为信号分析与处理的一种非常重要的手段和工具。
Matlab 是数学计算过程中可以用到的一种数学工具，其中，电子计算机为应用数学管理实际问题创造了必不可少的物质条件，如今，在电子技术高度发达的数学领域，技术科学中最有研究价值的是数学研究领域，而数学领域中数值分析和数学建模成为了重点中的重点。
1. 光学傅立叶变换
光学傅立叶变动是光学中信息处理的根底，其中所运用的道理便是通过光学中透镜来实现傅立叶的变换效应,相干光的照射会使透镜L 后焦面上光场的复振幅分布即是位于透镜物体前焦面的复振幅函数F(x,y)的立叶变换F(u,v)。将两个傅立叶变换串联就构成了典型的4厂(厂为透镜焦距)光学处理系统：Ll 的前焦面称为输入平面或物平面,Ll 的后焦面与L2的前焦面重合,称为频谱面或傅立叶面,L2的后焦面为输出平面或像平面。待处理的图像信号可以是胶片或干板,也可以是空间光学调制器。相干光(激光)经准直后形成宽束平行光照射位于输入平面上振幅透过率为厂(x,y)的胶片(相当于灰度为厂(x,y)的图像),由于Ll 的傅立叶变换效应,在频谱面上形成图像的傅立叶空间频谱,经L2的再次变换(傅立叶逆变换)在像平面上则形成原图的反像。若在频谱面上设置适当的滤波器或相关器,就能实现对图像的各种处理,如设置高通、低通、带通或方向滤波器作高通、低通、带通或方向滤波；设置逆滤波器进行图像复原；以匹配滤波器(滤波器的复振幅透过率为特定信号空间频谱的复数共轭)或光学相关器实现图像的相关、匹配和目标检测。设计更复杂的光学处理系统还可以完成广义傅立叶变换、维纳变换、小波变换及神经网络等较复杂的处理运算
2. 傅里叶变换分析
由(,)g x y 表示物体的波动散布，相干的前提之下用(,g x y ）表现xy 平面上的复式振幅，那么其所对应的模便是每一个点对应点的振幅，辐角则代表每一个点的初相位，由g(,)x y 作傅里叶变换。
,)(,)e x p {-2()}x y x y
f f
g x y j f
x f y d x d y
π-
=+??G( 作傅立叶变换，可以由G(,)x y f f 将(,)g x y 代替出：
(,=(,)exp{2()}x y x y x y g x y G f f j f x f y df df π+??）来，由式子可以看出，物函数(,)g x y 可
被分解为无穷多个不同频率的(,x y f f ）和无穷多个不同方位取向的(tan θ=
)y x
f f ，外加不同
权重{(,)x y x y G f f df df }指数基元，由物平面上复振幅公式：0(,)exp{(cos cos )U x y U jkx y αβ=+)，指数基元exp{2()}x y j f x f y π+。相当于方向是cos α=λx f ,cos β=y f λ为一个单位振幅的平面波，故此，该傅立叶变换式由可以理解为：物函数(,)g x y 可以被看作无数振幅差异的（|(,)x y x y G f f df df |），方向不同（cos α=x f λ
,cos β=y f λ）平面波相干叠加的结果，或
者是， ()g x y ，可以被分解为无数个振幅与方向均有不同方向的无数平面波。由方向和振幅的不同，体现其相干性。 同理，由非相干条件可知，()g x y ，是xy 上强度的分布，它是一种实函数，由逆向傅立叶变换可写成物函数：
(,)2|(,)|c o s {2()(}x y x y x y x y
g x y G f f f x f y f f d f d f
πψ∞
-∞
=++??，
）， 也就是在非相干照射的条件下，光强分布(,)g x y 可以被分解成不同幅值，不同空间取向，不同频率的余弦形式的强度分布，换言之，它也可以被解析成无穷个幅值相同，且相互之间
方向对称的平面波。
对于上述理解，在傅立叶光学中相当重要，我们在理解和研究傅立叶光学时，考察光学系统对物函数(,)g x y 的作用情况，就可以考察它对每个基元的作用情况，然后再加以作用分析，或者说可以分别考察系统的不同方面平行光的作用，然后再加以综合分析。
3. 傅里叶光学变换的方法
持续的傅里叶变更可以将平方的可积函数()f t 表示成为复指数函数的积分或者也
可表示为级数样式。
-(){()}=()iwt
F w f f t f t e
dt
∞
-∞
=?
这是将频域的函数()F W 表示时间域的函数()f t 的积分形式。连续傅立叶逆变换为：
1
1
()[()]()2iwt
f t f F W F w e
dw
π
∞
--∞
==
?
这是把时间域的函数表示成频率域的函数。一般可认为，
()f t 为原函数，而()F W 为傅
立叶的像函数，原函数和像函数之间构成傅立叶变换的一个变换对。
4. 光学傅立叶变换的本质
光学傅立叶变换的本质即公式的变换，即：
()()jwt
F W f t e dt
∞
-
-∞
=
?
也可以把傅立叶变换成另外一种形式
1()(),2jwt
F W f t e π
=??
由此可得，傅立叶函数变更的本质就是内积的变换和内积的运算，如果给定的三角函数是一个完备的正交函数集，那么在频率不相等的条件之下，给定三角函数间的内积为0，只有在
频率相等的条件之下，给此三角函数做内积运算时，其得到内积不为0。
1212()12,2()
j t
j j t
t
e
e
e
dt πδΩ-ΩΩΩ??==Ω-Ω?
5 .MATLAB在傅立叶光学中的应用
（1）傅立叶光学告诉我们夫琅禾费衍射场强度分布就得等于屏函数的功率谱，因此我们可以直接在光屏上进行傅立叶变更，之后通过处理就可得到衍射的图样。物体中图像的生成可以直接运用矩阵运算来完成，也可以直接利用WINDOWS中下载的图画工具，生成一幅黑
白图像，使用命令函数imread()对其进行输入命令，输入的图像为一个二维矩阵，此二维矩阵必须要求很大，之后再利用MATLAB函数中fft2（）对此矩阵进行二位离散傅里叶变更，以得到图像屏谱，该屏谱为一个复矩阵，然后运用取模函数abs()对该复数取模，得其屏谱，之后再利用函数FFTSHIFT对需要取模的矩阵进行屏谱位移，这是由于变换后的二维矩阵其直流分量是位于图像的周边角，况且图像交换矩阵的1、3象限和2、4象限，对直流分量移到谱中间，以便来使FFT屏谱可视效果与所得到的实际图像相吻合，最后利用imread（）将所需函数图像显示出来。
(2) 假设任意输入一个函数，然后，输出函数的傅里叶变更函数，然后得振幅频率。
x=sin(2*pi*t); %任意输入一个函数。
y=fft(x); %傅里叶变换函数。
plot(abs(y)); %振幅频率。
函数（function）表示一种对应关系，由每个输入值对应唯一输出值。这是一种一一对应的关系，由一个集合里的任意一个元素对应于另一个（可能相同的）集合里的唯一对应的元素。设一个函数中所对应输入值的输出值x的标准符号为f(x)。定义域是函数的三要素之一，对应法则的作用对象，所有属于函数输入值的集合为定义域，包含所有的输出值的集合为值域。定义在非空数集之间的映射称为函数。许多程序设计语言中，可以将一段经常需要使用的代码封装起来，在需要使用时可以直接调用，这就是程序中的函数。
在一定前提之下，傅里叶变更能满足某个函数表示成三角函数（正弦和/或余弦函数）或者满足其积分的线性组合的需求。在不同的研究领域和不同前提之下，傅里叶变更具多种不同的形式变换，如连续傅里叶变换和离散傅里叶变换。
以正弦波为例来进行说明:
1．用Matlab数学工具生成正弦波,之后作用形成并显示出时域波形图。
2. 对波形进行FFT傅里叶变换，显示频率波谱图形，其中所涉及到的采样率，频率、数据长度自选。
3．做出上述信号的均方根图谱,功率图谱,以及对数均方根图谱。4．用IFFT傅立叶反变换恢复信号，并显示恢复的正弦信号时域波形图。
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);
xlabel('频率(Hz)');
ylabel('对数谱');
title('正弦信号y=2*pi*10t对数谱');
grid;
%用IFFT做傅里叶反变换恢复原始信号
xifft=ifft(y);
magx=real(xifft);
ti=[0:length(xifft)-1]/fs;
figure(1);
subplot(236);
plot(ti,magx);
xlabel('t');
ylabel('y');
title('通过IFFT转换的正弦信号波形');
grid
6. 结论：
通过光学傅里叶变换和matlab 数学工具的有效结合，通过matlab强大的信号处理功能光学中的傅立叶变换具有难度高，抽象性大的特点，对数据处理和波形图像的处理，而matlab 也具有强大的信号处理功能，结合光学中傅立叶变换、傅立叶级数、卷积定理的内容，通过matlab程序来体现光学中几个傅里叶调动的直观演示。
7、参考文献：
[1]谢家宁，陈伟成，赵建林，陈国杰，张潞英.Matlab在光学信息处理仿真实验中的应用[J].
物理实验，2004,24（6）：23-25
[2]欧攀，戴一堂，王爱明，柳强，李立京，伊小素.高等光学仿真（matlab版）——光导波，
激光[M].北京：北京航空航天大学出版社，2001
[3]任玉杰.数值分析及其MATLAB实现[M].北京：高等教育出版社，2007
[4]张志勇.精通Matlab[M].北京：北京航空航天大学出版社，2000
[5]熊万杰，陆建龙.MATLAB软件支持下得理论物理教学改革[J].内蒙古师范大学学报，
2003,16（5）.
[6]蔡履中，王成彦，周玉芳.光学（修订版）[M].济南：山东大学出版社，2002.
[7]陈怀琛.MATLAB及其在理工课程中的应用指南[M].西安：西安电子科技大学出版社，
2000.
[8]潘柏根. 信息光学仿真实验软件的研究[M].安徽：合肥工业大学出版社，2010.
[9]徐蓁. 计算机生成全息图的快速硬件算法与光电重建[M].南京：东南大学出版社. 2008.
[10]赵盾. 光学实验计算机仿真平台的构建[M].武汉：武汉大学出版社，2010.
[11]何福. 数字信号处理方法在地震工程中的应用研究[M].北京：中国地震局工程力学研
究所, 2007.
[12] 常永. 基于谐波小波分析的矿山机械故障诊断研究[M]. 河南：河南理工大学, 2012.
[13] 尹珅. 基于光学原理的图像存储和恢复[M].南京：中南大学出版社，2014.
[14] 吕宪魁. 大孔径单色傅里叶变换透镜的设计[M].重庆：重庆大学出版社，2005.
[15] 刘茂先. 高一学生对函数概念理解的调查研究[M].贵州：贵州大学出版社，2014.
Matlab implementation of optical fourier transform
Abstract:in the optical fourier transform,with the high difficulty and abstract characteristics
致谢：
历时将近两个月的时间终于将这篇论文写完，在论文的写作过程中遇到了无数的困难和障碍，都在同学和老师的帮助下度过了。尤其要强烈感谢我的论文老师，她对我进行了无私的指导和帮助，不厌其烦的帮助进行论文的修改和改进。另外，在校图书馆查找资料的时候，
图书馆的老师也给我很多方面的支持和帮助。在此，向帮助和指导过我的各位老师表示最衷心的感谢！
班级信工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];
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/9911847952.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/9911847952.html,/view/fea1e985b9d528ea81c779ee.html %% 输入参数个数检查
function [Spec,Freq]=STFT(Sig,nLevel,WinLen,SampFreq) %计算离散信号的短时傅里叶变换； % Sig 待分析信号； % nLevel 频率轴长度划分（默认值512）； % WinLen 汉宁窗长度（默认值64）； % SampFreq 信号的采样频率（默认值1）； if (nargin <1), error('At least one parameter required!'); end; Sig=real(Sig); SigLen=length(Sig); if (nargin <4), SampFreq=1; end if (nargin <3), WinLen=64; end if (nargin <2), nLevel=513; end nLevel=ceil(nLevel/2)*2+1; WinLen=ceil(WinLen/2)*2+1; WinFun=exp(-6*linspace(-1,1,WinLen).^2); WinFun=WinFun/norm(WinFun); Lh=(WinLen-1)/2; Ln=(nLevel-1)/2; Spec=zeros(nLevel,SigLen); wait=waitbar(0,'Under calculation,please wait...'); for iLoop=1:SigLen, waitbar(iLoop/SigLen,wait); iLeft=min([iLoop-1,Lh,Ln]); iRight=min([SigLen-iLoop,Lh,Ln]); iIndex=-iLeft:iRight; iIndex1=iIndex+iLoop; iIndex2=iIndex+Lh+1; Index=iIndex+Ln+1; Spec(Index,iLoop)=Sig(iIndex1).*conj(WinFun(iIndex2)); end; close(wait); Spec=fft(Spec); Spec=abs(Spec(1:(end-1)/2,:));
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
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;
%傅里叶变换 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('幅值频谱')
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）原理 正交级数的展开是其理论基础！将一个在时域收敛的函数展开成一系列不同频率谐波的叠加，从而达到解决周期函数问题的目的。在此基础上进行推广，从而可以对一个非周期函数进行时频变换。 从分析的角度看，他是用简单的函数去逼近（或代替）复杂函数，从几何的角度看，它是以一族正交函数为基向量，将函数空间进行正交分解，相应的系数即为坐标。从变幻的角度的看，他建立了周期函数与序列之间的对应关系；而从物理意义上看，他将信号分解为一些列的简谐波的复合，从而建立了频谱理论。 当然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
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对信号进行频谱分析及滤波 设计目的 要求学生会用MATLAB语言进行编程，绘出所求波形，并且运用FFT求对连续信号进行分析。 一、设计要求 1、用Matlab产生正弦波，矩形波，并显示各自的时域波形图； 2、进行FFT变换，显示各自频谱图，其中采样率、频率、数据长度自选，要求注明； 3、绘制三种信号的均方根图谱； 4、用IFFT回复信号，并显示恢复的正弦信号时域波形图。 二、系统原理 用FFT对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行频谱分析的信号是模拟信号和时域离散信号。频谱分辨率直接和FFT的变换区间N 有关，因为FFT能够实现频率分辨率是2π/N。 x(n)是一个长度为M的有限长序列，则x(n)的N点离散傅立叶变换为： X(k)=DFT[x(n)]= kn N W N n n x ∑ - = 1 ) ( ,k=0,1,...,N-1 N j e N W π 2 - = 逆变换：x(n) =IDFT[X(k)]= kn N W k X N n N - ∑ - = 1 ) ( 1 ,k=0,1,...,N-1 但FFT是一种比DFT更加快速的一种算法，提高了DFT的运算速率,为数字信号处理技术应用于各种信号处理创造了条件，大大提高了数字信号处理技术的发展。本实验就是采用FFT，IFFT对信号进行谱分析。 三、程序设计 fs=input('please input the fs:');%设定采样频率 N=input('please input the N:');%设定数据长度 t=0:0.001:1; f=100;%设定正弦信号频率
目录 用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
陕西科技大学实验报告 班级信工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）结果：
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 计算傅里叶变换
maｔlab实现信号的傅立叶变换 一、设计目的 ?1．熟悉和掌握mａtlab的基本使用方法，能够熟练运用matlab。 2．巩固信号与系统中的傅立叶变换内容,加深对这部分内容的理解。 二、设计任务 1.掌握matlａb的基本操作。 2.利用matlaｂ实现典型非周期信号的傅立叶变换，画出信号的时域图和频域图。 ?3.利用maｔlab实现傅立叶变换的基本性质。 三、设计原理 １．ｍatｌａb简介 MＡTLAB是ＭａthＷoｒks公司推出的一套高性能的数值计算和可视化软件，经过多年大量的、坚持不懈的改进，现在ＭATLAB已经更新至7．x版。MAＴLAB集数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个方便的、界面友好的用户环境。在这个环境下，对所要求解的问题,用户只需简单地列出数学表达式，其结果便以人们十分熟悉的数值或图形方式显示出来。 MATLAＢ可用来解决实际的工程和数学问题，其典型应用有:通用的数值计算,算法设计，各种学科（如自动控制、数字信号处理、统计信号处理)等领域的专门问题求解。ＭATＬＡＢ语言易学易用，不要求用户有高深的数学和程序语言知识，不需要用户深刻了解算法及编程技巧。MATＬAB既是一种编程环境,又是一种程序设计语言。这种语言与C、FＯRTRＡＮ等语言一样，有其内定的规则，但MAＴＬAB的规则更接近数学表示。使用更为简便,可使用户大大节约设计时间,提高设计质量。 2．ｍatlab2013ｂ基本界面介绍 maｔｌaｂ2013ｂ主界面窗口基本分为五个部分: 1）主菜单界面 在此界面我们只需要用到新建命令文件和对程序进行间断调试的功能 2）文件查看窗口,双击可快速打开文件
离散信号的变换（MATLAB 实验） 一、实验目的 掌握用Z 变换判断离散系统的稳定与否的方法，掌握离散傅立叶变换及其基本性质和特点，了解快速傅立叶变换。 二、实验内容 1、已经系统函数为 5147.13418.217.098.2250 5)(2342-++--+=z z z z z z Z H (1) 画出零极点分布图，判断系统是否稳定； (2)检查系统是否稳定； (3) 如果系统稳定，求出系统对于u(n)的稳态输出和稳定时间b=[0,0,1,5,-50];a=[2,-2.98,0.17,2.3418,-1.5147]; subplot(2,1,1);zplane(b,a);title('零极点分布图'); z=roots(a); magz=abs(z) magz = 0.9000 0.9220 0.9220 0.9900 n=[0:1000]; x=stepseq(0,0,1000); s=filter(b,a,x); subplot(2,1,2);stem(n,s);title('稳态输出'); （1）因为极点都在单位园内，所以系统是稳定的。 （2）因为根的幅值（magz ）都小于1，所以这个系统是稳定的。 （3）稳定时间为570。 2、综合运用上述命令，完成下列任务。 (1) 已知)(n x 是一个6点序列： ???≤≤=其它,050,1)(n n x
计算该序列的离散时间傅立叶变换，并绘出它们的幅度和相位。 要求：离散时间傅立叶变换在[-2π,2π]之间的两个周期内取401个等分频率上进行数值求值。 n=0:5;x=ones(1,6); k=-200:200;w=(pi/100)*k; X=x*(exp(-j*pi/100)).^(n'*k); magX=abs(X);angX=angle(X); subplot(2,1,1);plot(w/pi,magX);grid;title('幅度'); subplot(2,1,2);plot(w/pi,angX);grid;title('相位'); (2) 已知下列序列： a. ,1000),52.0cos()48.0cos()(≤≤+=n n n n x ππ； b ．)4sin()(πn n x =是一个N ＝32的有限序列； 试绘制)(n x 及它的离散傅立叶变换 )(k X 的图像。 a ． n=[0:1:100];x=cos(0.48*pi*n)+cos(0.52*pi*n); subplot(2,1,1);plot(n,x);title('x(n)的图像'); X=dft(x,101); magX=abs(X); subplot(2,1,2);plot(n,magX);title('丨X(k)丨的图像');
一、实验目的 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;
MATLAB 在离散傅立叶变换(DFT)中的应用 一、序列的移位和周期延拓运算。 已知)()8.0()(8n R n x n =，利用MATLAB 生成并图示序列),(),(m n x n x -和)())((8n R n x N ),())((8n R m n x N -其中为周期的延拓。以表示8)())((,0,248n x n x N m N <<= 解：MATLAB 程序清单如下： N=24; M=8; m=3;% 设移位值为3 n=0:N-1; xn=0.8.^n.*(n>=0 & n
实验13 离散傅里叶变换的性质 (完美格式版，本人自己完成，所有语句正确，不排除极个别错误，特别适用于山大，勿用冰点等工具下载，否则下载之后的word格式会让很多部分格式错误，谢谢) XXXX学号姓名处XXXX 一、实验目的 1 加深对离散傅里叶变换(DFT)基本性质的理解。 2 了解有限长序列傅里叶变换(DFT)性质的研究方法。 3 掌握用MATLAB语言进行离散傅里叶变换性质分析时程序编写的方法。 二、实验内容 1 线性性质。 2 循环移位性质。 3 循环折叠性质。 4 时域和频域循环卷积特性。 5 循环对称性。 三、实验环境 MA TLAB7.0 四、实验原理 1 线性性质 如果两个有限长序列分别为x1(n)和x2(n)，长度分别为N1和N2，且 y(n)＝ax1(n)＋bx2(n) (a、b均为常数) 则该y(n)的N点DFT为 Y(k)＝DFT［y(n)］＝aX1(k)＋bX2(k) 0≤k≤N－1 其中：N＝max［N1，N2］，X1(k)和X2(k)分别为x1(n)和x2(n)的N点DFT。 例13-1已知x1(n)＝［0，1，2，4］，x2(n)＝［1，0，1，0，1］，求： (1)y(n)＝2x1(n)＋3x2(n)，再由y(n)的N点DFT获得Y(k)； (2)由x1(n)、x2(n)求X1(k)、X2(k)，再求Y(k)＝2X1(k)＋3X2(k)。 用图形分别表示以上结果，将两种方法求得的Y(k)进行比较，由此验证有限长序列傅里叶变换(DFT)的线性性质。 解MA TLAB程序如下： >> xn1=[0,1,2,4]; %建立xn1序列 >> xn2=[1,0,1,0,1]; %建立xn2序列 >> N1=length(xn1);N2=length(xn2); >> N=max(N1,N2); %确定N >> if N1>N2 xn2=[xn2,zeros(1,N1-N2)]; %对长度短的序列补0 >> elseif N2>N1 xn1=[xn1,zeros(1,N2-N1)]; >> end >> yn=2*xn1+3*xn2; %计算yn >> n=0:N-1;k=0:N-1;
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：矩形函数傅里叶变换的理论结果