前半部分内容转自：http://www.ilovematlab.cn/thread-119939-1-1.html

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/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的结果的模值如图所示。



                      图1 FFT结果
    从图中我们可以看到，在第1点、第51点、和第76点附近有比较大的值。我们分别将这三个点附近的数据拿上来细看：
1点： 512+0i
2点： -2.6195E-14 - 1.4162E-13i 
3点： -2.8586E-14 - 1.1898E-13i

50点：-6.2076E-13 - 2.1713E-12i
51点：332.55 - 192i
52点：-1.6707E-12 - 1.5241E-12i

75点：-2.2199E-13 -1.0076E-12i
76点：3.4315E-12 + 192i
77点：-3.0263E-14 +7.5609E-13i
   
 很明显，1点、51点、76点的值都比较大，它附近的点值都很小，可以认为是0，即在那些频率点上的信号幅度为0。接着，我们来计算各点的幅度值。分别计算这三个点的模值，结果如下：
1点： 512
51点：384
76点：192
 

按照公式，可以计算出

直流分量为：512/N=512/256=2；

50Hz信号的幅度为：384/(N/2)=384/(256/2)=3；

75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5。

可见，从频谱分析出来的幅度是正确的。（需要注意的是，由于该信号的频率分辨率为Fs/N=1Hz，也即fft能分辨的最大的信号细节为1Hz，若模拟的cos三角信号中有小数频段的信号时，则无法正确计算出的该频段的幅度，同时，若该信号使用fft的阶数不为信号的采样频率的整数倍，也会出现此类情况。这也是fft的栅栏效应导致的）

然后再来计算相位信息。直流信号没有相位可言，不用管它。

先计算50Hz信号的相位，atan2(-192, 332.55)=-0.5236,结果是弧度，换算为角度就是180*(-0.5236)/pi=-30.0001。

再计算75Hz信号的相位，atan2(192, 3.4315E-12)=1.5708弧度，换算成角度就是180*1.5708/pi=90.0002。可见，相位也是对的。
根据FFT结果以及上面的分析计算，我们就可以写出信号的表达式了，它就是我们开始提供的信号。

总结：

假设采样频率为Fs，采样点数为N，做FFT之后，某一点n（n从1开始）表示的频率为：Fn=(n-1)*Fs/N；

n点的模值除以N/2就是对应该频率下的信号的幅度（对于直流信号是除以N）；

n点的相位即是对应该频率下的信号的相位。相位的计算可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值，范围从-pi到pi。要精确到xHz，则需要采样长度为1/x秒的信号，并做FFT。

要提高频率分辨率，就需要增加采样点数，这在一些实际的应用中是不现实的，需要在较短的时间内完成分析。

解决这个问题的方法有频率细分法，比较简单的方法是采样比较短时间的信号，然后在后面补充一定数量的0，使其长度达到需要的点数，再做FFT，这在一定程度上能够提高频率分辨力。具体的频率细分法可参考相关文献。

[附录：本测试数据使用的matlab程序]
close all; %先关闭所有图片
Adc=2;  %直流分量幅度
A1=3;   %频率F1信号的幅度
A2=1.5; %频率F2信号的幅度
F1=50;  %信号1频率(Hz)
F2=75;  %信号2频率(Hz)
Fs=256; %采样频率(Hz)
P1=-30; %信号1相位(度)
P2=90;  %信号相位(度)
N=256;  %采样点数
t=[0:1/Fs:N/Fs]; %采样时刻

%信号
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
%显示原始信号
plot(S);
title('原始信号');

figure;
Y = fft(S,N); %做FFT变换
Ayy = (abs(Y)); %取模
plot(Ayy(1:N)); %显示原始的FFT模值结果
title('FFT 模值');

figure;
Ayy=Ayy/(N/2);   %换算成实际的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; %换算成实际的频率值
plot(F(1:N/2),Ayy(1:N/2));   %显示换算后的FFT模值结果
title('幅度-频率曲线图');

figure;
Pyy=[1:N/2];
for i="1:N/2"
Pyy(i)=phase(Y(i)); %计算相位
Pyy(i)=Pyy(i)*180/pi; %换算为角度
end;
plot(F(1:N/2),Pyy(1:N/2));   %显示相位图
title('相位-频率曲线图');

泄露与窗函数

每次FFT变换只能对有限长度的时域数据进行变换，因此，需要对时域信号进行信号截断。信号截断有两种，一种是周期截断，一种是非周期截断，哪怕原始信号是周期信号。若周期截断，则FFT频谱为单一谱线，得到的频率成分为原始信号的真实频率，并且幅值与原始信号的幅值相等。
    若为非周期截断，截断后的信号起始时刻和结束时刻的幅值不等，将这个信号再进行重构，在连接处信号的幅值不连续，出现跳跃；对截断后的信号做FFT，频谱出现拖尾，峰值处的频率与原始信号的频率相近，但并不相等。另一方面，峰值处的幅值已不再等于原始信号的幅值，为原始信号幅值的64%（矩形窗的影响）。而幅值的其他部分（36%幅值）则分布在整个频带的其他谱线上。拖尾现象这种非常严重的误差，称为泄漏，是数字信号处理所遭遇的最严重误差。现实世界中，在做FFT分析时，很难保证截断的信号为周期信号，因此，泄漏不可避免。为了将这个泄漏误差减少到最小程度（注意是减少，而不是消除），我们需要使用加权函数，也叫窗函数。加窗主要是为了使时域信号似乎更好地满足FFT处理的周期性要求，减少泄漏。非周期截断的信号与窗函数相乘得到的信号起始点与最末点达到相同（比如都为0），变成一个类似周期截断的信号。窗函数只能减少泄漏，不能消除泄漏。

功率谱估计

功率谱密度（PSD），单位为：unit^2/Hz代表单位频率上信号的能量，所以是密度谱，幅值代表频段内的有效值平方。

分段加窗函数的含义。分段之后，需要补0至和原先信号一样长，这就相当于每小段加矩形窗，矩形窗不仅有1值还有0值。在连续的世界里非常理所当然的事情，在离散的世界里就发生了变化。如果不补0到原来长度，频率分辨率就大大降低了。

将一个信号从时域通过FFT变换到频域，得到的直接结果就是所谓的频谱，频谱是复数形式，有幅值和相位。由于频谱是复数形式，包含相位信息，当信号中包含不相关的噪声成分时，由于噪声成分的相位是杂乱无序的，那么多次线性平均之后，可以将不相关的噪声平均掉。

功率谱密度PSD表征的是单位频率上的能量分布。它等于自功率谱除以频率分辨率，因此，它的单位为（信号单位^2/Hz）。由于自谱是实数，因此，功率谱也是实数，可进行线性平均。它只有RMS格式。

自谱或称为自功率谱，本质是由频谱计算得到的，它是复数频谱乘以它的共轭。因此，自谱是实数，没有相位信息。由于它是实数，因此可以进行线性平均。

由于它是复数频谱与它的共轭的乘积，因此自谱有平方形式，平方形式的自谱称为自功率谱Power。对平方形式的自谱再求平方根，对应为线性形式，称为线性自功率谱AutoPower Linear。当用线性自谱来表示一个随机信号的结果时，由于信号的总能量是一定的，当采用不同的频率分辨率会导致幅值大小不一样。因为，频率分辨率高，则谱线越密集，因而分配到每条谱线上的能量就少，对应的幅值就低。当用PSD（PSD=自功率谱除以频率分辨率）来表示随机信号的结果时，由于PSD是实数，所以没有相位信号，只有幅值信息。

PSD计算时的步骤为
1 对每一分段数据进行FFT变换，并求它的幅值谱
2 对幅值谱进行平方
3 将双边谱转化为单边谱？？
4 除以频率分辨率delt(f)=1/T=fs/nfft

举个例子：
幅值为1，频率为16Hz的正弦信号，使用1024Hz采样，2048点进行功率谱密度计算，频率分辨率为1024/2048=0.5Hz，求出的功率谱单边谱在第32根谱线处的值为1(P1(33)=1)，解释为：信号FFT变换后得到的双边谱，幅值分别为0.5，平方后为0.25，转化为单边乘2为0.5，再除以频率分辨率为1。将1乘以0.5（频率分辨率），正好为该信号有效值0.707的平方。

   
fs=1024;
second=2;
N = second*fs;
order=2048;

dt = 1/fs;
t = 0:dt:second-dt;%时间

fw=16;
Cz = sin(2*pi*fw*t);

figure(1)
subplot(211);
plot(t, Cz)
y=Cz;

yfft = fft(y(1:N),order);
P2 = abs(yfft/N);
P1 = P2(1:N/2+1);%单边谱
P1(2:end-1) = 2*P1(2:end-1);%由于P1(1)是直流

因为实数信号的双边幅值谱都是对称的，双边谱中包含负频率，在物理系统中是没有的，因此用单边谱就够了，这时候把负频率成分附加到相应的正频率成分，也就是在双边谱的基础上乘以２。

1、matlab中fft函数之后，信号幅度值变大，要得到真实幅度值大小，要将变换后结果除以N再乘以2（直流即零频率处除以N），除以N得到双边频谱，乘以2得到单边频谱。如果原始信号就是复数信号，则fft得到的就是单边频谱不用乘以2。

2、对fs采样频率信号，使用fft函数之后，最高频率为采样频率一半，fft之后的值关于半采样率共轭对称。

3、对时域信号做傅里叶变换后对其值求模，得到的谱叫做频谱的幅度谱，它的纵坐标没有物理量纲，这个幅度谱与相位谱在一起称为频谱。
     再对傅里叶变换后的值做进一步处理，具体叫什么谱视纵坐标而定，如果纵坐标是能量，就叫能量频谱即频域能量谱，这里区分时域的能量，时域能量是时间对应的能量，能量频谱则是频率对应的能量；如果是功率，就叫功率频谱即频域功率谱。而能量频谱与功率频谱之间是一个线性关系，他们差了一个信号时长的倍数。
     再说功率谱密度，它和功率谱本质上是一个东西，功率谱指的是整个一个谱，是一个频率序列，功率谱密度这是指某个具体频率值对应的功率值。人们习惯上将二者不加区分，经常混用，都是指代频率点处的功率值。

clc;
clear;

fs=1000;%采样率
sec=100;%100秒的数据
t = 0:1/fs:sec-1/fs;

y = 4*sin(2*pi*50*t)+1*sin(2*pi*20*t)+ 3*sin(2*pi*80*t);
% y = 2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*70*t+pi*90/180);

nfft=2000;
y1 = y(1:nfft);
N=length(y1);
xaxis=fs* (0: N/2-1)/N;

figure(1)

subplot(4,1,1);
plot( y1);
title('原始信号');

%fft后，此时求出的Xk除以N得到双边频谱（幅度谱和相位谱都在这里），乘以2得到单边频谱。如果原始信号就是复数信号，则fft得到的就是单边频谱不用乘以2。
Xk=fft(y1,nfft);%若y1的长度小于nfft，则fft函数会自动补0至长度为nfft。若长度大于nfft，则截断，只对前nfft个数据进行fft处理
biXk=abs(Xk)./nfft;%双边幅度谱
oneXk(1) = biXk(1);
oneXk(2:nfft/2) = biXk(2:nfft/2).*2;%单边幅度谱 需要注意的是，fft的栅栏效应导致的幅度低于实际值的问题
subplot(4,1,2);
plot(xaxis, oneXk);
 title('单边幅度谱');
 
fEnergy = Xk.*conj(Xk)/nfft;%双边能量谱，此算式即Xk^2/nfft   会随着nfft的增大而增大
fEnergy2 = fEnergy(1: nfft/2);
fEnergy2(2:nfft/2) = fEnergy2(2:nfft/2).*2;%单边能量谱
 subplot(4,1,3);
 plot(xaxis, fEnergy2);
 title('单边能量谱');
 
fxP2 = (Xk./ nfft).* (conj(Xk)./ nfft)%双边功率谱  同能量谱差了一个信号时长的倍数
fxP1 = fxP2(1: nfft/2);
fxP1(2:end) = 2* fxP1(2:end);%单边功率谱
subplot(4,1,4);
plot(xaxis, fxP1);
title('单边功率谱');

curhz=50;%当前频率的功率谱密度
a=fxP1(curhz/fs*nfft+1);

相关链接：

PSD的单位及计算方法[转]_weixin_34273046的博客-CSDN博客

PSD的单位及计算方法[转] - 马语者 - 博客园

FFT 频谱分辨率_不积跬步，无以至千里-CSDN博客_频谱分辨率

10-11:功率谱估计及频率分辨率 - 简书

说说脑电信号的频谱、能量、功率谱、功率谱密-燕山徐一 | 个人博客