深入浅出解释FFT（四）——fft分析信号频率和相位_深入浅出解析fft

很基础的问题往往很重要，做仿真时候有一个点的差错都会导致结果的错误。在网上找了前人写的东东，总结下希望对大家有帮助，让大家少走一些弯路。

1.  信号的时域采样点N和频域采样点数相同

%##################################################################

clear all; close all;
Adc =1.25;                  %直流分量幅度
A1 =1;                      %频率F1信号的幅度
A2 =0.25;                    %频率F2信号的幅度
F1 =100;                   %信号1频率(Hz)
F2 =1000;                   %信号2频率(Hz)
Fs =5120;                  %采样频率(Hz)

Ts =1/Fs;                     %时域采样间隔，采样周期
P1 =-30;                   %信号1相位(度)
P2 =90;                    %信号2相位(度)
N =256;                     %采样点数，请注意时域采样N点，FFT时也是N点！
t = [0 : 1/Fs : N/Fs];        %采样时刻，注意不是序列是真正的时间！t = [0:N]*Ts

% t = [0 : N]//Fs;  %共采样257个点，为了好看起见多采样了最后一个点，最后一个点在实际中应该是下一个采样周期的第一个点

%生成信号
%S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180); %考虑相位
signal1=A1*sin(2*pi*F1*t); % 周期T1 = 2*pi/(2*pi*F1) = 0.01s，采样周期Ts = 0.00019531 < (T1)/2

signal2=A2*sin(2*pi*F2*t); % 周期T2 = 2*pi/(2*pi*F2) = 0.001s，采样周期Ts = 0.00019531 < (T2)/2
S=A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t);                      %无直流偏置
S=Adc+A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t);                  %直流偏置
%S=Adc+[A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t)]/2;

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 模值');

F = ([1:N] - 1) * Fs /N;              %换算成实际的频率值， N*Fs/2 对应着w = pi
% 因为MATLAB中FFT的变换矩阵不是一个酉矩阵(Unitary Matrix)，该阵除以1/sqrt(N)就是个酉矩阵。故经过变换后对信号有放大作用，

% 所以要在fft处理后结果除以N/2来，修正此“放大”作用。但是结果在直流的那一点是错误的，实际上直流应该除以N修正。

Ytemp = 2/N * fft(S,N); % 做FFT变换后除以N/2修正

figure; plot(F, abs(Ytemp));title('实际幅度-频率曲线图 --- abs前修正');

figure;
Ayy = Ayy / (N /2);                   %换算成实际的幅度
Ayy(1) = Ayy(1) /2;                   %能看出来为什么直流信号除以N了吗？我开始没看出来啊

%假设采样频率为Fs，采样点数为N，做FFT之后，某一点n（n从1开始）表示的频率为：Fn=(n-1)*Fs/N；该点的模值除以N/2就是对应

%该频率下的信号的幅度（对于直流信号是除以N）

 

 

根据公式验证 --- 某点n所表示的频率为： Fn=(n-1)*Fs/N，因为没有fftshift

                                    f1=（6-1）*5120/256=100 （Hz） ---- 验证结果正确

                                    f2=（51-1）*5120/256=1000（Hz）---- 验证结果正确

 根据公式验证 --- 对于n=1点的信号，是直流分量，幅度即为A1/N   --- A1=320/256=1.25 --- 验证结果正确

                对于n点（n≠1，且n<=N/2）  幅度A --- A=An/(N/2)=128/（256/2）=1 --- 验证结果正确

plot(F(1:N/2),Ayy(1:N/2));            %显示换算后的FFT模值结果，只看0~N*Fs/2正频谱部分
title('实际幅度-频率曲线图 --- abs后修正');

figure;

% 相位的计算可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值，范围从-pi到pi。

% atan2(500, 148)=x,结果是弧度，换算为角度就是180*(-x)/pi=相位值。

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结果以及上面的分析计算，--- 频率,幅度.相位....我们就可以写出信号的表达式了，它就是我们开始提供的信号。------ 信号重建

2.  信号的时域采样点N和频域采样点数NFFT不同

%如果采样的时间周期是信号周期的倍数，可能泄露就会避免。频谱泄露是由于非整周期采样导致的。

% 我以前有个错误的观念：信号的时域采样点N必须和FFT的计算点数NFFT相同，才会给处理和解释带来便利。

% 原来是模型的不同产生的影响，采用模拟频率 f 建模不会产生末尾补零使得FFT频率不一致的问题，看下面的程序。

clear; close all;

% f1 = 20Hz, f2 = 40Hz %若 f1 = 15 可以观察频谱泄露情况

%  （1）数据个数Ndata=32，FFT所用的采样点数NFFT=32
fs=240; %采样频率
Ndata=32; %数据长度
N=32; %FFT的数据长度
n=0:Ndata-1;t=n/fs;   %数据对应的时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);   %时间域信号，x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); 可以发现频谱泄露
y=fft(x,N);   %信号的Fourier变换
mag=abs(y);    %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=32');grid on;

% （2）数据个数Ndata=32，FFT所用的采样点数NFFT=126
Ndata=32;   %数据个数
N=128;     %FFT采用的数据长度
n=0:Ndata-1;t=n/fs;   %时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=128');grid on;

% （3）数据个数Ndata = 136，FFT所用的采样点数NFFT=128
Ndata=136;   %数据个数
N=128;     %FFT采用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真实频率
subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=128');grid on;

% （4）数据个数Ndata=136，FFT所用的采样点数NFFT=512
Ndata=136;    %数据个数
N=512;    %FFT所用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真实频率
subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;
% （5）数据个数Ndata=136，FFT所用的采样点数NFFT=4096
Ndata=136;    %数据个数
N=4096;    %FFT所用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真实频率

figure,plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;

3.  信号的时域采样点N和频域采样点数NFFT不同 -- 比较两种不同的信号模型

clear; close all;

% 该信号的数字周期 N = 8，模拟周期T = N*Ts = 0.008s ，实际频率 f = 125 Hz

w = pi/4;
N = 2*pi/w; % N = 8 
n = 0:N-1; % n = 0:N 可能会更好看一些，但是要清楚第N+1点可是下一个采样周期的第一个点
x = sin(n*w); % 不含有采样信息的x，补零的FFT对频谱有影响
h=plot(n, x,'-o'); %注意n没有定标，没有物理含义！
Fs = 1000; %采样频率为1000Hz
Ts = 1/Fs; 
t = n*Ts; % t 时间序列的给法永远只有这么一种，请铭记！
T = N*Ts; %模拟周期T
f = w*Fs/(2*pi); % 信号的真实频率f
figure; plot(t, x)
% freq = n*Fs/N - Fs/2; %频率序列的定标，注意此式子是本来面目，不要给此式子穿个马甲就不认识了，下面是此式子的变形
% freq = (n/N - 1/2)*Fs;        freq = (n -N/2)/N * Fs; % 注意0~N-1时第N 点对应着 Fs， N/2 点对应着 Fs/2， Fs/N 就是频率分辨率
% freq = (n - N/2)/N;  就是归一化频率，从-0.5 到 0.5 - 1/N , 看看不是 -0.5 ~ 0.5 哦，实际结合硬件编程就要这样干

Nfft = 1024; % 此处 Nffft/N = 128，所以要对结果的频谱进行128点的抽取，才可以得到正确的频率定标！
X_Nfft = fftshift(abs(fft(x,Nfft)));
freqNormalized = ((0:Nfft-1) - Nfft/2)/Nfft;
freq_ = freqNormalized*Fs;
freq = freqNormalized(1:128:end)*Fs; % 结果和上面的结果完全相同 -500  -375 -250  -125     0   125  250   375

X_ = X_Nfft(1:128:end);

figure; plot(freq_, X_Nfft);grid

figure; plot(freq , X_);grid

%-------------------------------------------------------------------------
Ndata=136;    %数据个数
N=4096;    %FFT所用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x =0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t); %含有采样信息，补零后FFT对信号的频谱没影响
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真实频率

figure,plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;