用Matlab对信号进行傅里叶变换实例

目录
用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
虽然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 点。

Matlab 的傅里叶变换实例
1. 傅里叶变换的定义傅里叶变换从数学上的定义，为
F(w)=int(x*exp(-jwt),t=-inf...inf) 其中，int表示积分，t是时间，x是时域信号，inf 表示无穷，exp 表示指数运算。

其含义说的是给一个
无限长的时域信号和一个频点w，可以唯一确定一个复数F。

于是， F 和w 就有了这种对应关系，考虑到 F 是个复数。

F 的绝对值和w 关系叫幅频， F 的幅角和w 关系叫相频。

help fft 可以知道这个和数学上的傅里叶不一样，因为计算机是离散的！因为计算机的时域信号存储量是有限的！比如等时采样得到的信号，高频分量是不可能获得的，对于比较大的w 将无法计算。

于是，fft 这样计算傅里叶变换：把时域信号进行周期延拓，取一组w，就是时域信号的周期及该周期的二分之一，三分之一，直到n 分之一，其中n是一个周期内的采样点。

这样做的结果，就是对一段有限长的时域信号，将其长度作为基频率，分析基频和高频含量。

当然，能分析到的最高频为n 次谐波，再高次谐波由于香农定理而无法体现。

3. 写一个数学定义傅里叶变换的程序将有限长时域信号不延拓，时域信号外的时间内，认为信号为零。

于是获得无限长时域信号，取频点若干，分析其傅里叶变换。

考虑到matlab 对于由描点法定义的函数，数值积分时常用的方法有：矩形法，
梯形法。

一下代码采用梯形法，算例如下：
clear
clc
%% 输入信号
t=0:1e-3:20;% 时域信号的时间范围x=sin(t)+sin(1.5*t+1)+5*cos(0.5*t)+2*randn(size(t));% 时域信号x w=[0:1e-2:2];% 想要观察的频率范围
%% 预定义
y=w;
a=w;
j=sqrt(-1);% 先定义变量维度，提高运算速度
%% 计算频点
for i=1:length(w)
f=trapz(t,x.*exp(-j*w(i)*t));
y(i)=abs(f);
a(i)=angle(f);
end
%% 输出
subplot(3,1,1),
plot(t,x)
subplot(3,1,2),
plot(w,y)
subplot(3,1,3),
plot(w,a)
算例中的时域信号，里有三个正余弦分量，一个干扰分量。

等时采样，并认为采样频率足够高，即得到的信号是连续信号。

Matlab 方波傅立叶变换画出频谱图
代码如下： clc;clear;close all; fs=30; % 采样频率 T=1/fs;
t=0:T:2*pi;
A=2;P=4; y=A*square(P*t); subplot(2,1,1),plot(t,y) title(' 方波信号 ')
Fy=abs(fft(y,512));
f2=fs*(0:256)/512;
subplot(2,1,2),plot(f2,Fy(1:257)) title('频谱图 '); set(gcf,'unit','normalized','position',[0 0 1 1]); set(gca,'xtick',0:0.6:8); axis([0,8,0 300]);
方波大概 0.6 左右。

角频率 W=4 ，所以频率 f=4/（2*pi ）=2/pi 没错。

求信号的 512点 FFT ，点数越大 （2的整数幂 ），频谱越精细。

把 512 改成 1024 后，频率变到 1.2 左右了。

代码：
clc;clear;close all; fs=30; % 采样频率
N=input （
' 输入 FFT 点数： '）; T=1/fs; t=0:T:2*pi;
A=2;P=4; y=A*square(P*t); subplot(2,1,1),plot(t,y) title(' 方波信号 ') Fy=abs(fft(y,N)); f2=fs*(0:N/2)/N; subplot(2,1,2),plot(f2,Fy(1:(N/2+1))) title('频谱图 ');
set(gcf,'unit','normalized','position',[0 0 1 1]); set(gca,'xtick',0:0.6:8); axis([0,8,0 300]);
输入N 分别为128 ，256 ，512，1024 看看，是不是点越多，曲线越平滑。