短时傅立叶变换试验
为了克服傅立叶变换的时频局部化方面的不足，也是为了对时域信号作局部分析，
D.Gabor 于1946年提出了窗口傅立叶变换（简记为WFT ）。
WFT 的公式形式
()(,)()()j w t R G f w b f t w t b e d t -=-?()(,)()()j w t R G f w b f t w t b e d t
-=-? 其中，实函数w(t)为是时窗函数，窗函数w(t)具有较强的衰减性，所以要精心选择窗函数。
下面是一个短时傅立叶变换的代码程序
function timefreq(x,Nw,window)
% 待分析信号，行向量，Nw 时窗宽度
subplot(2,2,1);
plot(real(x));%描绘待分析信号
X=fft(x);%快速傅里叶变换
X=fftshift(X);%调整0频位置
subplot(2,2,2);
plot(abs(X));%描绘幅度谱
Lap=Nw/2;%重叠宽度
Tn=(length(x)-Lap)/(Nw-Lap);%计算分段数目
nfft=2^ceil(log2(Nw));%做fft 的点数
TF=zeros(Tn,nfft);%时频矩阵
for i=1:Tn
if(strcmp(window,'rec'))
Xw=x((i-1)*10+1:i*10+10);%加窗矩形处理
elseif(strcmp(window,'Hamming'))
Xw=x((i-1)*10+1:i*10+10).*Hamming(Nw)';%加hamming 处理 elseif(strcmp(window,'Blackman'))
Xw=x((i-1)*10+1:i*10+10).*Blackman(Nw)';%加black 处理 elseif(strcmp(window,'Gauss'))
Xw=x((i-1)*10+1:i*10+10).*Gauss(Nw)';%加Gauss 处理 else
return;
end
temp=fft(Xw,nfft);%求fft
temp=fftshift(temp);%调整0频位置
TF(i,:)=temp;%保存分段fft 结果
end
%绘制时频分析结果
subplot(2,2,3);
fnew=((1:nfft)-nfft/2)/nfft;
tnew=(1:Tn)*Lap;
[F,T]=meshgrid(fnew,tnew);
mesh(T,F,abs(TF));
xlabel('n');ylabel('w');zlabel('Gf'); subplot(2,2,4);
contour(T,F,abs(TF));
xlabel('n');ylabel('w');
例子：clc ;clear;
N=400;
x=zeros(1,N);
T=0:N-1;
x=exp(j*4*pi*(T/80).^2);
figure(1);
timefreq(x,20,’rec’);
figure(2);
timefreq(x,20,’Blackman’);
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,:));
短时傅立叶变换试验 为了克服傅立叶变换的时频局部化方面的不足，也是为了对时域信号作局部分析，D.Gabor 于1946年提出了窗口傅立叶变换（简记为WFT ）。 WFT 的公式形式 ()(,)()()j w t R G f w b f t w t b e d t -=-?()(,)()()j w t R G f w b f t w t b e d t -=-? 其中，实函数w(t)为是时窗函数，窗函数w(t)具有较强的衰减性，所以要精心选择窗函数。 下面是一个短时傅立叶变换的代码程序 function timefreq(x,Nw,window) % 待分析信号，行向量，Nw 时窗宽度 subplot(2,2,1); plot(real(x));%描绘待分析信号 X=fft(x);%快速傅里叶变换 X=fftshift(X);%调整0频位置 subplot(2,2,2); plot(abs(X));%描绘幅度谱 Lap=Nw/2;%重叠宽度 Tn=(length(x)-Lap)/(Nw-Lap);%计算分段数目 nfft=2^ceil(log2(Nw));%做fft 的点数 TF=zeros(Tn,nfft);%时频矩阵 for i=1:Tn if(strcmp(window,'rec')) Xw=x((i-1)*10+1:i*10+10);%加窗矩形处理 elseif(strcmp(window,'Hamming')) Xw=x((i-1)*10+1:i*10+10).*Hamming(Nw)';%加hamming 处理 elseif(strcmp(window,'Blackman')) Xw=x((i-1)*10+1:i*10+10).*Blackman(Nw)';%加black 处理 elseif(strcmp(window,'Gauss')) Xw=x((i-1)*10+1:i*10+10).*Gauss(Nw)';%加Gauss 处理 else return; end temp=fft(Xw,nfft);%求fft temp=fftshift(temp);%调整0频位置 TF(i,:)=temp;%保存分段fft 结果 end %绘制时频分析结果 subplot(2,2,3); fnew=((1:nfft)-nfft/2)/nfft; tnew=(1:Tn)*Lap; [F,T]=meshgrid(fnew,tnew); mesh(T,F,abs(TF)); xlabel('n');ylabel('w');zlabel('Gf');
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
%傅里叶变换 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('幅值频谱')
%%%%%%%%%% 线性调频信号%%%%%%%%%%%%%%%%%%%%%%%%%% clear all k=50; fs=250; n=0:1000-1; x=exp(j*pi*k*(n/fs).^2); subplot(221); plot(n/fs,real(x)); title('线性调频信号'); xlabel('时间'); ylabel('振幅'); %%%%%%%%%%%%%% 线性调频信号频谱%%%%%%%%%%%%%%%%%%%%%%%%% P=fft(x); subplot(222); plot(n/length(n)*fs,abs(P)); title('线性调频信号频谱'); xlabel('频率'); ylabel('振幅'); %%%%%%%%%%%%%%% 短时傅立叶变换%%%%%%%%%%%%%%%%%%%%%%% N=40; NFFT=512; N1=round((1000-N)*2/N+1); w=gausswin(N); %w=hanning(N); %w=rectwin(N); %w=hamming(N); %w=blackman(N); w=rot90(w); STFT=zeros(N1,NFFT); for i=1:N1; xt=x((i-1)*N/2+1:(i-1)*N/2+N); xt=xt.*w; %加高斯窗 STFT(i,:)=fft(xt,NFFT); end fn=(1:NFFT)*fs/NFFT; tn=(1:N1)*N/2/fs; [F,T]=meshgrid(fn,tn); subplot(223); mesh(F,T,abs(STFT)); title('STFT'); xlabel('频率'); ylabel('时间');
M A t l a b傅里叶变换实 验报告 文档编制序号：[KKIDT-LLE0828-LLETD298-POI08]
班级信工142 学号 22 姓名何岩实验组别 实验日期室温报告日期成绩 报告内容：(目的和要求，原理，步骤，数据，计算，小结等) 1.求信号的离散时间傅立叶变换并分析其周期性和对称性； 给定正弦信号x(t)=2*cos(2*pi*10*t),fs=100HZ,求其DTFT。(a)代码： f=10;T=1/f;w=-10::10; t1=0::1;t2=0::1; n1=-2;n2=8;n0=0;n=n1::n2; x5=[n>=]; 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 *min(x2) *max(x2)]); xlabel('x(n)');ylabel('x(n)'); title('原信号x1'); xlabel('t');ylabel('x1'); subplot(2,2,3);stem(t2,x2); axis([0 1 *min(x2) *max(x2)]); title('原信号采样结果x2'); xlabel('t');ylabel('x2'); subplot(2,2,2);stem(n,x5); axis([0 1 *min(x5) *max(x5)]); xlabel('n');ylabel('x2'); title('采样函数x2'); subplot(2,2,4);stem(t2,x4); axis([0 1 +*min(x4) *max(x4)]); xlabel('t');ylabel('x4'); title('DTFT结果x4'); （b）结果： 2.用以下两个有限长序列来验证DTFT的线性、卷积和共轭特性；
课程设计任务书 学生姓名:专业班级: 指导教师:工作单位: 题目: 连续时间信号傅里叶级数分析及MATLAB实现 初始条件： MATLAB 6.5 要求完成的主要任务： 深入研究连续时间信号傅里叶级数分析的理论知识，利用MA TLAB强大的图形处理功能，符号运算功能以及数值计算功能，实现连续时间周期信号频域分析的仿真波形。 1.用MATLAB实现周期信号的傅里叶级数分解与综合。 2.用MATLAB实现周期信号的单边频谱及双边频谱。 3.用MATLAB实现典型周期信号的频谱。 4.撰写《MATLAB应用实践》课程设计说明书。 时间安排： 学习MATLAB语言的概况第1天 学习MATLAB语言的基本知识第2、3天 学习MATLAB语言的应用环境，调试命令，绘图能力第4、5天 课程设计第6-9天 答辩第10天 指导教师签名：年月日 系主任（或责任教师）签名：年月日
目录 摘要................................................................................................................................................ I Abstract .......................................................................................................................................... II 绪论. (1) 1 MATLAB简介 (2) 1.1 MATLAB语言功能 (2) 1.2 MATLAB语言特点 (2) 2 傅里叶级数基本原理概要 (4) 2.1 周期信号的傅里叶分解 (4) 2.2 三角形式和指数形式傅里叶级数及各系数间的关系 (4) 2.3 周期信号的频谱 (5) 3 用MATLAB实现周期信号的傅立叶级数分解与综合 (6) 3.1 合成波形与原波形之间的关系 (6) 3.2 吉布斯现象 (6) 4 用MATLAB实现周期信号的单边频谱及双边频谱。 (8) 4.1 单边，双边（幅度，相位）频谱及其关系 (8) 4.1.1单边，双边（幅度，相位） (8) 4.1.2 单边，双边频谱关系 (9) 4.2以单边幅度频谱为例，研究脉冲宽度与频谱的关系 (10) 4.3以单边幅度频谱为例，研究脉冲周期与频谱的关系 (11) 5用MATLAB实现典型周期信号的频谱 (13) 5.1 周期方波脉冲频谱的MATLAB实现 (13) 5.2 周期三角波脉冲频谱的MATLAB 实现 (14) 6 小结及心得体会 (17) 参考文献 (18) 附录： (19)
傅里叶变换m a t l a b代 码 Document serial number【UU89WT-UU98YT-UU8CB-UUUT-UUT108】
%傅里叶变换 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
短时傅里叶变换短时傅里叶变换（STFT，short-time Fourier transform，或short-term Fourier transform)）是和傅里叶变换相关的一种数学变换，用以确定时变信号其局部区域正弦波的频率与相位。它的思想是：选择一个时频局部化的窗函数，假定分析窗函数g(t)在一个短时间间隔内是平稳（伪平稳）的，移动窗函数，使f(t)g(t)在不同的有限时间宽度内是平稳信号，从而计算出各个不同时刻的功率谱。短时傅里叶变换使用一个固定的窗函数，窗函数一旦确定了以后，其形状就不再发生改变，短时傅里叶变换的分辨率也就确定了。如果要改变分辨率，则需要重新选择窗函数。短时傅里叶变换用来分析分段平稳信号或者近似平稳信号犹可，但是对于非平稳信号，当信号变化剧烈时，要求窗函数有较高的时间分辨率；而波形变化比较平缓的时刻，主要是低频信号，则要求窗函数有较高的频率分辨率。短时傅里叶变换不能兼顾频率与时间分辨率的需求。短时傅里叶变换窗函数受到W.Heisenberg不确定准则的限制，时频窗的面积不小于2。这也就从另一个侧面说明了短时傅里叶变换窗函数的时间与频率分辨率不能同时达到最优。 短时距傅里叶变换 维基百科，自由的百科全书 汉漢▼
[编辑]与傅里叶转换在概念上的区别 将信号做傅里叶变换后得到的结果，并不能给予关于信号频率随时间改变的任何信息。以下的例子作为说明： 傅里叶变换后的频谱和短时距傅里叶转换后的结果如下： 傅里叶转换后, 横轴为频率(赫兹) 短时距傅里叶转换, 横轴为时间( 秒 ) ， 纵轴 为频率(赫兹) 由上图可发现，傅里叶转换只提供了有哪些频率成份的信息，却没有提供时间信息；而短时傅里叶转换则清楚的提供这两种信息。这种时频分析的方法有利于频率会随着时间改变的信号(例如:音乐信号、语音信号等)分析。 [编辑]定义 [编辑]数学定义 简单来说，在连续时间的例子，一个函数可以先乘上仅在一段时间不为零的窗函数再进行一维的傅里叶变换。再将这个窗函数沿着时间轴挪移，所得到一系列的傅里叶变换结果排开则成为二维表象。数学上，这样的操作可写为： 另外也可用角频率来表示：
应用MATLAB对信号进行频谱分析及滤波 设计目的 要求学生会用MATLAB语言进行编程，绘出所求波形，并且运用FFT求对连续信号进行分析。 一、设计要求 1、用Matlab产生正弦波，矩形波，并显示各自的时域波形图； 2、进行FFT变换，显示各自频谱图，其中采样率、频率、数据长度自选，要求注明； 3、绘制三种信号的均方根图谱； 4、用IFFT回复信号，并显示恢复的正弦信号时域波形图。 二、系统原理 用FFT对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行频谱分析的信号是模拟信号和时域离散信号。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现频率分辨率是2π/N。 x(n)是一个长度为M的有限长序列，则x(n)的N点离散傅立叶变换为： N?1?2?kn)(nx j?W W NN e?0?n N X(k)=DFT[x(n)]=,k=0,1,...,N-1N?11?kn?)(WXk N N0?n x(n) =IDFT[X(k)]= 逆变换：,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;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f*t); figure(1); subplot(211); plot(t,x);%作正弦信号的时域波形 axis([0,0.1,-1,1]); title('正弦信号时域波形'); z=square(50*t); subplot(212) plot(t,z) axis([0,1,-2,2]); title('方波信号时域波形');grid;
用Matlab实现快速傅立叶变换 FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 现在就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此啰嗦了。 采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。 假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。 假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。 下面以一个实际的信号来做说明。假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)。式中cos 参数为弧度，所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？我们来看看FFT的结果的模值如图所示。
目录 用Matlab对信号进行傅里叶变换 (2) Matlab的傅里叶变换实例 (5) Matlab方波傅立叶变换画出频谱图 (7)
用Matlab对信号进行傅里叶变换 1.离散序列的傅里叶变换DTFT(Discrete Time Fourier Transform) 代码： 1 N=8; %原离散信号有8点 2 n=[0:1:N-1] %原信号是1行8列的矩阵 3 xn=0.5.^n; %构建原始信号，为指数信号 4 5 w=[-800:1:800]*4*pi/800; %频域共-800----+800 的长度（本应是无穷，高频分量很少，故省去） 6 X=xn*exp(-j*(n'*w)); %求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.离散傅里叶变换DFT(Discrete Fourier Transform)
与1中DTFT不一样的是，DTFT的求和区间是整个频域，这对 结果图：
分析：DFT只是DTFT的现实版本，因为DTFT要求求和区间无穷，而DFT只在有限点内求和。 3.快速傅里叶变换FFT（Fast Fourier Transform） 虽然DFT相比DTFT缩减了很大的复杂度，但是任然有相当大的计算量，不利于信息的实时有效处理，1965年发现的DFT解决了这一问题。 实现代码： 1 N=64; %原离散信号有8点 2 n=[0:1:N-1] %原信号是1行8列的矩阵 3 xn=0.5.^n; %构建原始信号，为指数信号 4 Xk=fft(xn,N); 5 subplot(221); 6 stem(n,xn); 7 title('原信号'); 8 subplot(212); 9 stem(n,abs(Xk)); 10 title('FFT变换') 效果图： 分析：由图可见，fft变换的频率中心不在0点，这是fft算法造成的，把fft改为fftshift可以将频率中心移到0点。
西南交通大学峨眉校区 （作业小论文） 工程测试技术课程设计 短时傅里叶变换及其谱图分 姓名：xxxx 学号：2wwwww 班级：wwww 专业：工程机械 2013.03.20
短时傅里叶变换及其谱图分析 摘要：本文讨论了有噪信号的短时傅里叶变换STFT及其谱围．分析和仿真结果表明，受 白噪声污染的信号的STFT可以无偏估计原信号的STFT，而其谱图可以对愿信号的谱图作有偏 估计，估计方差是有限的，且是时间和频率的函数．在短窗的情况下，求得了该方差上限的 近似表示． 关键词：短时傅里叶变换谱图噪声污染信号估计 1.前言 信号的短时傅里叶变换STFT是最早提出的一种时。频二维表示方法，它采用加窗的复正 弦作为基函数，也称为加窗傅里叶变换。由于它采用单一的分析窗处理所有频率分量，在时- 频平面内所有点上的分辨率是相同的，因而适合于准平稳信号的处理。STFT 简单易实现，许 多联合时．频分析的应用都是由它开始的。尽管STFT按定义属于线性变换，但在各种实际应 用中常采用它的能谱分布表示方法，这就是基于短时傅里叶变换的谱图Spectrogram)表示。 谱图定义为STFT的模平方，它是二次型时．频分布，尽管不满足时一频边缘条件，但可以认为 是信号能量在时．频平面上的分布。谱图已经在信号检测，语音处理等方面得到了广泛应用 [1Ⅱ2】。 谱图具有非线性性质，对于多分量信号将产生类似于Wigner分布中的交叉项干扰，从而 引入了模糊，影响信号分析结果。在利用谱图对信号的谱估计中，加性噪声的影响使信号具有 了多分量特性．可能使得估计产生较大偏差。本文就确定性信号受自噪声污染后的STFT及其 谱图的最小方差估计问题进行了分析。文中第二部分做了理论推导，求得了有噪信号的sTFT 及其谱图的均值和方差，第三部分对短窗的情况作了近似分析，最后给出了一例简单的仿真结 果。 2.傅立叶变换的提出
关于傅立叶变换，无论是书本还是在网上可以很容易找到关于傅立叶变换的描述，但是大都是些故弄玄虚的文章，太过抽象，尽是一些让人看了就望而生畏的公式的罗列，让人很难能够从感性上得到理解，最近，我偶尔从网上看到一个关于数字信号处理的电子书籍，是一个叫Steven W. Smith, Ph.D.外国人写的，写得非常浅显，里面有七章由浅入深地专门讲述关于离散信号的傅立叶变换，虽然是英文文档，我还是硬着头皮看完了有关傅立叶变换的有关内容，看了有茅塞顿开的感觉，在此把我从中得到的理解拿出来跟大家分享，希望很多被傅立叶变换迷惑的朋友能够得到一点启发，这电子书籍是免费的，有兴趣的朋友也可以从网上下载下来看一下，URL地址是： https://www.doczj.com/doc/bd10367085.html,/pdfbook.htm 要理解傅立叶变换，确实需要一定的耐心，别一下子想着傅立叶变换是怎么变换的，当然，也需要一定的高等数学基础，最基本的是级数变换，其中傅立叶级数变换是傅立叶变换的基础公式。 二、傅立叶变换的提出 让我们先看看为什么会有傅立叶变换？傅立叶是一位法国数学家和物理学家的 名字，英语原名是Jean Baptiste Joseph Fourier(1768-1830), Fourier对热传递很感兴趣，于1807年在法国科学学会上发表了一篇论文，运用正弦曲线来描述温度分布，论文里有个在当时具有争议性的决断：任何连续周期信号可以由一组适当的正弦曲线组合而成。当时审查这个论文的人，其中有两位是历史上著名的数学家拉格朗日(Joseph Louis Lagrange, 1736-1813)和拉普拉斯(Pierre Simon de Laplace, 1749-1827)，当拉普拉斯和其它审查者投票通过并要发表这个论文时，拉格朗日坚决反对，在近50年的时间里，拉格朗日坚持认为傅立叶的方法无法表示带有棱角的信号，如在方波中出现非连续变化斜率。法国科学学会屈服于拉格朗日的威望，拒绝了傅立叶的工作，幸运的是，傅立叶还有其它事情可忙，他参加了政治运动，随拿破仑远征埃及，法国大革命后因会被推上断头台而一直在逃避。直到拉格朗日死后15年这个论文才被发表出来。 谁是对的呢？拉格朗日是对的：正弦曲线无法组合成一个带有棱角的信号。但是，我们可以用正弦曲线来非常逼近地表示它，逼近到两种表示方法不存在能量差别，基于此，傅立叶是对的。 为什么我们要用正弦曲线来代替原来的曲线呢？如我们也还可以用方波或三角波来代替呀，分解信号的方法是无穷的，但分解信号的目的是为了更加简单地处理原来的信号。用正余弦来表示原信号会更加简单，因为正余弦拥有原信号所不具有的性质：正弦曲线保真度。一个正弦曲线信号输入后，输出的仍是正弦曲线，只有幅度和相位可能发生变化，但是频率和波的形状仍是一样的。且只有正弦曲线才拥有这样的性质，正因如此我们才不用方波或三角波来表示。
傅里叶变换的本质 傅里叶变换的公式为 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 下面从公式解释下傅里叶变换的意义 因为傅里叶变换的本质是内积，所以f(t)和t j e ω求内积的时候，只有f(t)中频率为ω的分量才会有内积的结果，其余分量的内积为0。可以理解为f(t)在t j e ω上的投影，积分值是时间从负无穷到正无穷的积分，就是把信号每个时间在ω的分量叠加起来，可以理解为f(t)在t j e ω上的投影的叠加，叠加的结果就是频率为ω的分量，也就形成了频谱。 傅里叶逆变换的公式为 ωωπωd e F t f t j ?+∞∞-=)(21 )( 下面从公式分析下傅里叶逆变换的意义 傅里叶逆变换就是傅里叶变换的逆过程，在)(ωF 和t j e ω-求内积的时候，)(ωF 只有t 时刻的分量内积才会有结果，其余时间分量内积结果为0，同样积分值是频率从负无穷到正无穷的积分，就是把信号在每个频率在t 时刻上的分量叠加起来，叠加的结果就是f(t)在t 时刻的值，这就回到了我们观察信号最初的时域。 对一个信号做傅里叶变换，然后直接做逆变换，这样做是没有意义的，在傅里叶变换和傅里叶逆变换之间有一个滤波的过程。将不要的频率分量给滤除掉，然后再做逆变换，就得到了想要的信号。比如信号中掺杂着噪声信号，可以通过滤波器将噪声信号的频率给去除，再做傅里叶逆变换，就得到了没有噪声的信号。 优点：频率的定位很好，通过对信号的频率分辨率很好，可以清晰的得到信号所包含的频率成分，也就是频谱。 缺点：因为频谱是时间从负无穷到正无穷的叠加，所以，知道某一频率，不能判断，该频率的时间定位。不能判断某一时间段的频率成分。 例子： 平稳信号：x(t)=cos(2*pi*5*t)+cos(2*pi*10*t)+cos(2*pi*20*t)+cos(2*pi*50*t)
— 短时傅立叶变换试验 为了克服傅立叶变换的时频局部化方面的不足，也是为了对时域信号作局部分析，于1946年提出了窗口傅立叶变换（简记为WFT ）。 WFT 的公式形式 ()(,)()()jwt R Gf w b f t w t b e dt -=-?()(,)()()jwt R Gf w b f t w t b e dt -=-? 其中，实函数w(t)为是时窗函数，窗函数w(t)具有较强的衰减性，所以要精心选择窗函数。 下面是一个短时傅立叶变换的代码程序 function timefreq(x,Nw,window) % 待分析信号，行向量，Nw 时窗宽度 、 subplot(2,2,1); plot(real(x));%描绘待分析信号 X=fft(x);%快速傅里叶变换 X=fftshift(X);%调整0频位置 subplot(2,2,2); plot(abs(X));%描绘幅度谱 Lap=Nw/2;%重叠宽度 Tn=(length(x)-Lap)/(Nw-Lap);%计算分段数目 ' nfft=2^ceil(log2(Nw));%做fft 的点数 TF=zeros(Tn,nfft);%时频矩阵 for i=1:Tn if(strcmp(window,'rec')) Xw=x((i-1)*10+1:i*10+10);%加窗矩形处理 elseif(strcmp(window,'Hamming')) Xw=x((i-1)*10+1:i*10+10).*Hamming(Nw)';%加hamming 处理 elseif(strcmp(window,'Blackman')) ( Xw=x((i-1)*10+1:i*10+10).*Blackman(Nw)';%加black 处理 elseif(strcmp(window,'Gauss')) Xw=x((i-1)*10+1:i*10+10).*Gauss(Nw)';%加Gauss 处理 else return; end temp=fft(Xw,nfft);%求fft temp=fftshift(temp);%调整0频位置 ； TF(i,:)=temp;%保存分段fft 结果 end
一、傅立叶变化的原理； （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
短时傅里叶变换和小波变换 吴桐 （西南交通大学峨眉校区机械工程系铁道车辆一班学号20116432）摘要：短时傅里叶变换（STFT，short-time Fourier transform，或 short-term Fourier transform)）是和傅里叶变换相关的一种数学变换，用以确定时变信号其局部区域正弦波的频率与相位。小波变换与Fourier变换相比，是一个时间和频域的局域变换因而能有效地从信号中提取信息，通过伸缩和平移等运算功能对函数或信号进行多尺度细化分析（Multiscale Analysis），解决了Fourier变换不能解决的许多困难问题。 关键词：傅里叶变换；实质；缺陷；小波变换 0引言 傅立叶变换是信号分析技术的基础，它在分析平稳信号时起着至关重要的作用。傅立叶变换是一种全局的变换，只能获得信号的整个频谱，并不能反映某一局部时间内信号的频谱特性。在许多科学领域的实验和工程测量中，普遍存在着非平稳信号。针对传统的傅立叶变换无法表达信号的时频局域性质，人们提出了一系列新的信号分析理论，其中以短时傅立叶和小波变换应用最为广泛。 1傅里叶变换的实质与缺陷 傅里叶变换的实质是将一个信号分离为无穷多多正弦/复指数信号的加成，也就是说，把信号变成正弦信号相加的形式——既然是无穷多个信号相加，那对于非周期信号来说，每个信号的加权应该都是零——但有密度上的差别。所以，傅里叶变换之后，横坐标即为分离出的正弦信号的频率，纵坐标对应的是加权密度。对于周期信号来说，因为确实可以提取出某些频率的正弦波成分，所以其加权不为零——在幅度谱上，表现为无限大——但这些无限大显然是有区别的，所以我们用冲激函数表示已经说过，傅里叶变换是把各种形式的信号用正弦信号表示，因此非正弦信号进行傅里叶变换，会得到与原信号频率不同的成分——都是原信号频率的整数倍。这些高频信号是用来修饰频率与原信号相同的正弦信号，使之趋近于原信号的。所以说，频谱上频率最低的一个峰（往往是幅度上最高的），就是原信号频率。傅里叶变换把信号由时域转为频域，因此把不同频率的信号在时域上拼接起来进行傅里叶变换是没有意义的——实际情况下，隔一段时间采集一次信号进行变换，才能体现出信号在频域上随时间的变化。 基于傅里叶(Fourier)变换的信号频域表示,揭示了时间函数和频谱函数之间的内在联系，在传统的平稳信号分析和处理中发挥了极其重要的作用，很多理论研究和应用研究都把傅里叶变换当作最基本的经典工具来使用。但是傅里叶变换存在着严重的缺点：用傅里叶变换的方法提取信号频谱时，需要利用信号的全部时域信息，这是一种整体变换，缺少时域定位功能，因此必须对其加以改进。 信号的瞬时频率，表示了信号的谱峰在时间-频率平面上的位置及其随时间的变化情况，一般平稳信号的瞬时频率为常数，而非平稳信号的瞬时频率是时间t的函数。从傅里叶变换