matlab fft 函数 全面详解 横坐标定标 幅度恢复等

相关解释都在注释里了，写的很详细。

% matlab version 2020b
% time 20211108
% zzx fft for newer
%%
clc;clear;close all;
%% 以下步骤是循序渐进的，是按照推理过程编写的：
% -----------------step1 确定信号幅度
A1=7;
A2=3; 
% -----------------step2 确定信号的频率
f0=20;%单位：Hz or kHz or其他，这个不影响，为了描述方便，后续都不再带单位。
f1=5;
% -----------------step3 根据信号频率以及奈奎斯特采样定理，确定采样频率最小值
fs = 100;   % fs>min_fs=2*max(f0,f1),即，信号的采样频率须大于信号频率最大值的两倍
            % 此处如果采样频率过小，无法得出真实频域信号。
ts=1/fs;% 确定采样周期
% -----------------step4 对连续信号进行“采样”
t = 0:ts:200;%注意，时间本来是连续的，此处用ts进行量化，相当于对连续信号进行了采样
S1 = A1*cos(2*pi*f0*t);
S2 = A2*cos(2*pi*f1*t);
S = S1 + S2;
% -----------------step5 添加噪声
N = randn(1,length(S));
X = S + N;
% -----------------step6 fft
Y = fft(X);

%到上面这一步没什么问题，使人容易困惑的地方在于横坐标 纵坐标的意义
% -----------------step7 确定横坐标分辨率&绘图
len = length(Y);
x_axis_error = 0:len-1;% 【error】这样定横坐标是错误的，这只是表示了样点的索引，每个样点没有物理意义。
delta_f = fs/len;%每个样点的 频率分辨率（注意理解除法的意义）
x_axis = (0:len-1) * delta_f;%乘法左边是分配样点索引，乘法右边是说明每个样点的 频率分辨率。
figure(1);
plot(x_axis,abs(Y));
xlabel('Frequency'); ylabel('Amplitude'); title('FFT Result'); grid minor;
%绘出此图你会发现三个问题：
% 1、前两个峰值横坐标不完全=f0 f1。
%    这是因为频率分辨率的问题，横坐标其实是离散的，按照delta_f递增，所以无法确保峰值正好在f0 f1;
% 2、应该有两个峰值，却出现了四个，而且是中间轴对称的?
% 3、纵坐标为毛≠A1 A2?

% 问题2 直接给出数学推导解释(其实就是简单的三角变换，看不懂的记得问候你数学老师)。
temp1 = (A1/2) * (exp(1i*2*pi*f0*t) + exp(1i*2*pi*-f0*t));
bool_result = (S1 == temp1);%运行后发现两者相等。
% fft绘制的结果，其实是按照temp1的形式绘制的，有负频率。这里只给出了S1的推导，S2同理，所以可以看到四个峰值。并且：
% 直接做fft的结果，信号的前半部分对应频率[0,fs/2]，后半部分对应[-fs/2,0]。为了将零频点移到频谱中间，需要使用fftshift函数
% -----------------step8 fftshift
Y_shift = fftshift(Y);
x_axis_shift = (-len/2:len/2-1) * delta_f;
figure(2);
plot(x_axis_shift,abs(Y_shift));
xlabel('Frequency'); ylabel('Amplitude'); title('FFT Shift Result'); grid minor;

%问题3 
% -----------------step9 峰值恢复
% fft后，幅值增加了len倍（DFT公式），并且由于按照temp1绘制，峰值又分成了两份，因此真实的幅值为：
Y_real = fftshift(2*Y/len);
figure(3);
plot(x_axis_shift,abs(Y_real));%此时，横坐标为真实频率，纵坐标为真实幅度
xlabel('Frequency'); ylabel('Real Amplitude'); title('FFT Real Result'); grid minor;