前言
 
做OFDM通信少不了频谱分析，基带信号DA后的频谱，以及基带数字上变频后的DA信号都要频谱分析。我觉得其实做任何工程都是这样，先规定实施方案，然后仿真成功，再实际开发，不过也可以一边开发，一边仿真，开发结果要与仿真预期结果一致。
 所以分析与仿真工具MATLAB就很重要了，既可以仿真，又可以通过示波器或其他方法把实际信号采下来分析。
 
matlab使用FFT函数分析信号频谱
 
一般我使用的FFT分析频谱流程如下：
 
%% 两个频率分别为15HZ 和 20HZ 的正弦信号[1]
Fs=50;%采样率
f1=15;
f2=20;
t = 0:1/Fs:10-1/Fs; %500个点
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
figure;
plot(t,x);
y = fft(x); 
%将横坐标转化，显示为频率f= n*(fs/N)
f = (0:length(y)-1)*Fs/length(y);
figure;
plot(f,abs(y));
title('Magnitude');
%该变换还会生成尖峰的镜像副本，该副本对应于信号的负频率。
%为了更好地以可视化方式呈现周期性，可以使用 fftshift 函数对变换执行以零为中心的循环平移。
n = length(x);                         
fshift = (-n/2:n/2-1)*(Fs/n);
yshift = fftshift(y);
figure;
plot(fshift,abs(yshift));
 
 
 

  
  
 
 
其中有3个注意的点：
1.FFT的结果看的是频谱，所以怎么把横坐标的值从原来的FFT点数0:N-1转换为频率值呢？
 首先要引出频谱分辨率的概念，即分辨两个不同频率信号的最小间隔，FFT结果相邻点间的间隔。因为N点FFT对应采样率为fs的序列，其频率分辨率为,其中Ts为采样周期，T为整个序列的时间长度。有关频率分辨率的就不多说了。所以我们横坐标转换为：f = (0:length(y)-1)*Fs/length(y);
 
2.直接FFT的结果里怎么又多余的信号频率（镜像频率）图2？
 DFT具有对称性，因为其是周期序列DFS在一个周期内的点，时域序列是有限长实序列，DFT的结果的实部周期偶对称，虚部周期奇对称，也就是模值周期偶对称，相位周期奇对称。其实从奈奎斯特定律也可以看出，fs>=2f，fs的采样率最多也就显示fs/2的真实频率（感性理解哈哈）。
 所以程序处理方式就是周期延拓后取-N/2:N/2-1.用到函数fftshift()，结果如图3.如注释所述：
 %该变换还会生成尖峰的镜像副本，该副本对应于信号的负频率。
 %为了更好地以可视化方式呈现周期性，可以使用 fftshift 函数对变换执行以零为中心的循环平移。
 
其实这和设计数字滤波器IIR与FIR也一样，采样率为fs的信号，设计的滤波器的通带阻代也限制在0-fs/2内。
 
3.程序中的信号幅度值都是1，500点的FFT画出来的幅度值怎么变成了250，应该是1吧？
 是的，应该是1。所以怎么变换为1呢，注意到FFT的结果是偶对称的，且其反应的真实频谱是0-fs/2。所以需要的操作是直接取0-N/2的FFT结果，乘以2，然后除以N。即2*abs(y(1:N/2+1))/N，在上面的程序下接着写：
 注意到要除以N，也就是FFT的长度，为什么除以N，这个有很多说法，我自己理性理解的也不透彻，所以这里就不解释了自行百度（我只能结合本例子感性理解哈哈）。
 
%FFT的结果所要展现的真实的频谱幅值[2]
realy=2*abs(y(1:n/2+1))/n;
realf=(0:n/2)*(Fs/n);
figure;
plot(realf,realy);
 
 结果如图4，可见横坐标频率是0-25，纵坐标信号的幅度为1.
 
然而，有不少人（包括我）平时分析信号都是直接FFT画频谱：图2，且看幅值都是看相对的大小，或者有没有频率分量，就很少做图3与图4的变换。但是我不知道其中的缘由，今天终于明白了。
 
 
 

  
  
 
 
 
 
参考资料
 
[1]傅里叶变换 - MATLAB & Simulink - MathWorks 中国
 
[2]快速傅里叶变换- MATLAB fft - MathWorks 中国

实验代码仅供参考!
 
一、 对无限长时域离散信号进行频谱分析
 
 试用DFT对其进行频谱分析。 使用矩形窗进行截取，并分别
 选择三种不同的窗长，即N=16, 60, and 120。
 实验要求：
 1、用Matlab编程实现上述频谱分析；
 2、画出不同窗长截取后的信号的频谱；
 3、窗固定，分析不同窗长对频谱的影响；
 4、窗长固定，不同窗对频谱的影响。
 
close all;
clear all;
 
N1 = 16;
N2 = 60;
N3 = 120;
n = 1:200;
x = cos(pi*n/10) + sin(pi*n/6) + cos(2*pi*n/5) ;
 
%图1，时域x[n]------------------------------------
figure
stem(x)
title('时域')
xlabel('n')
ylabel('x[n]')
 
%图2，DFT,矩形窗----------------------------------
% N = 16
y1 = fft(x,N1);
Y1 = abs(fftshift(y1));
ly1 = length(y1);
k1 = (-ly1/2:ly1/2-1)/N1;
 
figure
subplot(311)
stem(k1,Y1)
title('频域（N=16,矩形窗）')
xlabel('k')
ylabel('X(k)')
 
%N = 60
y2 = fft(x,N2);
Y2 = abs(fftshift(y2));
ly2 = length(y2);
k2 = (-ly2/2:ly2/2-1)/N2;
 
subplot(312)
stem(k2,Y2)
title('频域（N=60，矩形窗）')
xlabel('k')
ylabel('X(k)')
 
%N = 120
y3 = fft(x,N3);
Y3 = abs(fftshift(y3));
ly3 = length(y3);
k3= (-ly3/2:ly3/2-1)/N3;
 
subplot(313)
stem(k3,Y3)
title('频域（N=120，矩形窗）')
xlabel('k')
ylabel('X(k)')
 
%图3，DFT,hamming窗----------------------------------
 
x1= x(1:N1).*hamming(N1)';
x2= x(1:N2).*hamming(N2)';
x3= x(1:N3).*hamming(N3)';
 
% N = 16
y11 = fft(x1,N1);
Y11= abs(fftshift(y11));
ly11 = length(y11);
k11 = (-ly11/2:ly11/2-1)/N1;
 
figure
subplot(311)
stem(k11,Y11)
title('频域（N=16，hamming窗）')
xlabel('f')
ylabel('X(f)')
 
%N = 60
y22 = fft(x2,N2);
Y22 = abs(fftshift(y22));
ly22 = length(y22);
k22 = (-ly22/2:ly22/2-1)/N2;
 
subplot(312)
stem(k22,Y22)
title('频域（N=60，hamming窗）')
xlabel('f')
ylabel('X(f)')
 
%N = 120
y33 = fft(x3,N3);
Y33 = abs(fftshift(y33));
ly33 = length(y33);
k33= (-ly33/2:ly33/2-1)/N3;
 
subplot(313)
stem(k33,Y33)
title('频域（N=120，hamming窗）')
xlabel('f')
ylabel('X(f)')
 

 
 
 结果分析:
 由图1-2可以看出，N=60和N=120时的频谱一样且正确,而N=16时的不正确,说明了:当矩形窗截取长度大于或等于信号周期时，可以正确显示频谱图;当矩形窗截取长度小于信号周期时，无法正确显示频谱图。
 此外,对比图1-2和图1-3,可以看出加矩形窗和Hamming窗的频谱有区别,但在对应的频率点都有最大值,说明Hamming窗影响并不是太大.
 
二、对连续信号进行谱分析
 试用DFT进行频谱分析，要求频率分辨率为1Hz。
 注：频率分辨率是指频谱分析中能够分辨的两个相邻频率点谱线的最小距离，即频率域的采样间隔
 实验要求：
 1、确定对x(t)采样的最大采样间隔和频域采样的最少采样点数
 2、编程实现上述频谱分析；
 3、画出该信号的幅频特性；（横坐标用信号的模拟频率！！！）
 4、分析该信号的幅频特性，考虑频域二、对连续信号进行谱分析
 试用DFT进行频谱分析，要求频率分辨率为1Hz。
 注：频率分辨率是指频谱分析中能够分辨的两个相邻频率点谱线的最小距离，即频率域的采样间隔
 实验要求：
 1、确定对x(t)采样的最大采样间隔和频域采样的最少采样点数
 2、编程实现上述频谱分析；
 3、画出该信号的幅频特性；（横坐标用信号的模拟频率！！！）
 4、分析该信号的幅频特性，考虑频域采样的采样点数对频谱的影响。
 
close all;
clear all;
t = 1;   %time =1s
fs = 300; %采样频率 fs>=2*fm,fm = 100Hz
N1 = 50;
T = 1/fs;
t1 = 0:T:t;
x = cos(200*pi*t1) + sin(100*pi*t1) + cos(50*pi*t1);
y = fft(x,fs);
Y = abs(fftshift(y));
 
figure
subplot(211)
stem(x)
title('时域（fs=300Hz）')
xlabel('n')
ylabel('x(nT)')
subplot(212)
yl = -length(y)/2:length(y)/2-1;
stem(yl,Y)
title('频域（fs=300Hz）')
xlabel('f')
ylabel('X(f)')

 

 
 
 结果分析:
 这里对时长为1s的连续信号进行时域采样,通过fft函数得到频域采样的结果.
 由x[n]的代数式可知截止频率fm = 100Hz,则采样频率fs>=200Hz,则时域最大采样间隔为1/200s.这里取采样频率为fs=300Hz.
 信号基频为50π,最高频率为200π，为4倍基频，则频域采样的最少采样点为9个。 对比图一和图二,采样频率不同,即采样点数也不同,可以看到都能在对应的频率点有峰值,而其他点基本为0.可见当满足以上的采样规则时,采样点数对频谱基本没有影响.但是,当不满足以上的采样规则,如fs =100Hz,在图2-3可以看出不能正确显示频谱.
 
三、
 
1）取一个周期，计算其DFT；
 2）将1）中的 x[n]补零，使 0≤n ≤99，计算其DFT；
 3）对x[n]取0≤n ≤99，计算其DFT。
 实验要求：
 1、编程实现上述三个步骤；
 2、分析上述三个步骤得到的频谱有什么区别和联系。
 注：DFT点数取成等于时域信号的点数，即X=fft(x)
 
close all;
clear all;
 
n = 1:100;
nl = length(n);
x = cos(0.48*pi*n) + cos(0.52*pi*n);
%0.48pi/2 =0.24pi; 0.24pi = 24pi/100 = 12*2pi/50  ,0.52pi = 13*2pi/50 ;
T = 50; %周期
%(1)
x1 = x(1:T);
y1 = fft(x1,T);
Y1 = abs(fftshift(y1));
%(2)
y2 = fft(x1,nl);
Y2 = abs(fftshift(y2));
%(3)
y3 = fft(x,nl);
Y3 = abs(fftshift(y3));
 
figure
subplot(221)
stem(x)
title('时域')
xlabel('n')
ylabel('x(n)')
 
subplot(222)
l1 = (-length(Y1)/2:length(Y1)/2-1)/T;
stem(l1,Y1)
title('(1)')
xlabel('f')
ylabel('X(f)')
 
subplot(223)
l2 = (-length(Y2)/2:length(Y2)/2-1)/nl;
stem(l2,Y2)
title('(2)')
xlabel('f')
ylabel('X(f)')
 
subplot(224)
l3 = (-length(Y3)/2:length(Y3)/2-1)/nl;
stem(l3,Y3)
title('(3)')
xlabel('f')
ylabel('X(f)')
 

 结果分析:
 对比(1)和(3)发现频谱一致,说明了去时域的一个周期(50个采样点)或两个周期(100个采样点)都可以正确显示频谱图;
 但是对x[n]补0后，补0 的部分当做了信号的一个周期内的一部分，所以频谱发生了变化。

 
基于MATLAB仿真的连续信号频谱分析.doc.doc
 
实验十：基于MATLAB仿真的连续信号频谱分析
 
一、实验目的
 
1、学会运用MATLAB求连续时间信号的傅里叶变换；
 
2、学会运用MATLAB求连续时间信号的频谱图。
 
二、实验原理
 
1、周期信号的分解
 
根据傅里叶级数的原理，任何周期信号都可以表示为三角级数的组合称为f(t)的傅
 
合成波形所包含的谐波分量越多，除间断点附近外，它越接近于原波形，在间断点附近，即使合成的波形所含谐波次数足够多，也仍存在约9%的偏【示例1】用正弦信号的叠加近似合成一个频率为50Hz，幅度为3的MATLAB程序如下：
 
clear all;
 
fs=5000;
 
t=0:1/fs:0.1;
 
f0=50;sum=0;
 
subplot(211)
 
for n=1:2:9
 
plot(t,4/pi*1/n*sin(2*pi*n*f0*t),'k');
 
hold on;
 
end
 
title('信号叠加前');
 
subplot(212)
 
for n=1:2:9;
 
sum=sum+4/pi*1/n*sin(2*pi*n*f0*t);
 
end
 
plot(t,sum,'k')
 
10.1所示。
 
2、非周期信号的频谱分析
 
信号的傅里叶变换定义为
 
(.1)
 
傅里叶反变换定义为
 
(.2)
 
信号的傅里叶变换主要包括MATLAB符号运算和MATLAB数值分析两种方法。
 
(1) MATLAB符号运算求解法
 
MATLAB符号数学工具箱提供了直接求解傅里叶变换和傅里叶反变换的函数fourier( )及ifourier( )。傅里叶变换的语句格式分为三种。
 
F=fourier(f )；表示符号函数f的Fourier变换，默认返回是关于的函数。
 
F=fourier(f ,v)；表示返回函数F是关于符号对象v的函数，而不是默认的，即。
 
F=fourier(f ,u,v); 是对关于u的函数f进行变换，返回函数F是关于v的函数，即。
 
类似的，傅里叶反变换的语句格式也分为三种。
 
要注意的是，函数fourier( )及ifoureir( )都是接收由sym函数所定义的符号变量或者符号表达式。【2】用MATLAB符号运算求解法求单边指数信号的傅里叶变换。
 
解：MATLAB源程序为
 
>> ft=sym('exp(-2*t)*Heaviside(t)');
 
>> Fw=fourier(ft)
 
运行结果为
 
Fw =
 
1/(2+i*w)
 
(2) 连续时间信号的频谱图
 
信号f(t)的傅里叶变换表达了信号在处的频谱密度分布情况，这就是信号的傅里叶变换的物理意义。一般是复函数，可以表示为。我们把与曲线分别称为非周期信号的幅度频谱和相位频谱。
 
【3】用MATLAB命令绘制单边指数信号幅度谱和相位谱。
 
解：MATLAB源程序为
 
ft=sym('exp(-2*t)*Heaviside(t)');
 
Fw=fourier(ft);
 
subplot(211)
 
ezplot(abs(Fw),[-2*pi,2*pi]),grid on
 
title('幅度谱')
 
phase=atan(imag(Fw)/real(Fw));
 
subplot(212)
 
ezplot(phase,[-2*pi,2*pi]),grid on
 
title('相位谱')
 
程序运行结果如图.2所示。
 
(3)MATLAB数值计算求解法
 
Fourier( )和ifourier( )函数的一个局限性是，如果返回函数中有诸如狄拉克函数等项，则用ezplot( )函数无法作图。对某些信号求变换时，其返回函数可能包含一些不能直接用符号表达的式子，甚至可能出现提示“未被定义的函数或变量”，因而也不能对其返回函数作图。此外，在很多实际情况中，尽管信号是连续的，但经过抽样所获得的信号则是离散的数值量，因此无法表示成符号表达式，此时不能应用Fourier( )函数对进行处理，而只能用数值计算法来近似求解。下面介绍傅里叶变换的数值计算法。
 
从傅里叶变换定义出发，有
 
(3)
 
当足够小时，上式的近似情况可以满足实际需要。对于时限信号，或者在所研究的时间范围内衰减到足够小，上式可研究有限n的取值。假设是因果信号，则有
 
(4)
 
对(4)式的角频率进行离散化，假设离散化后得到N个样值，即
 
， (5)
 
因此有
 
采用行向量，用矩阵表示为
 
(6)
 
式(6)提供了用MATLAB实

 
 
 
傅立叶变换是一种线性的积分变换，常在将信号在时域（或空域）和频域之间变换时使用，在物理学和工程学中有许多应用。因其基本思想首先由法国学者约瑟夫·傅里叶系统地提出，所以以其名字来命名以示纪念。在不同的研究领域，傅立叶变换具有多种不同的变体形式，如连续傅立叶变换和离散傅立叶变换。最初傅立叶分析是作为热过程的解析分析的工具被提出的。
 
 
 
傅里叶变换经常用来计算存在噪声的时域信号的频谱，下面就采用FFT方法来分析存在噪声的时域信号的频谱
 
 
这里假设数据采样的频率为1000Hz，
一个信号包含频率为50Hz、振幅为0.7的正弦波和频率为120Hz、振幅为1的正弦波，
噪声为零平均值的随机噪声
 
% 时域信号的频谱分析

clear;
clc;
Fs = 1000; % 采样频率
T = 1 / Fs; % 采样时间
L = 1000; % 信号长度
t = (0:L-1)*T; % 时间向量
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));
figure(1);
plot(Fs*t(1:50), y(1:50));
title('零平均值噪声信号');
xlabel('time (milliseconds)');
NFFT = 2^nextpow2(L);
Y = fft(y, NFFT)/L;
f = Fs/2*linspace(0, 1, NFFT/2);
figure(2);
plot(f, 2*abs(Y(1:NFFT/2)));
title('y(t) 单边振幅频谱');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
 
计算结果的图像如下所示