快速傅里叶变换FFT
fft在matlab中的函数就是fft,它式离散傅里叶变换的快速算法。fft的数学公式为:
X
[
k
]
=
1
N
∑
n
=
1
N
x
[
n
]
e
j
2
π
k
n
N
X[k] = \frac{1}{N}\sum_{n=1}^{N}x[n]e^{j\frac{2\pi kn}{N}}
X[k]=N1n=1∑Nx[n]ejN2πkn
原来的信号的序列长度为
N
N
N,则fft之后得到的结果依然是
N
N
N个点。
fft之后的
X
[
k
]
X[k]
X[k]表示了第
k
k
k个频率分量。那么如何将
k
k
k对应到真实频率上去呢?
真实频率
真实频率取决于两点:
- 得到原始序列 x [ n ] x[n] x[n]的采样频率;
- 采样点数;
例如采样点数是1000,采样频率是1000Hz,那么真实频率就是1,2,…,500Hz。
例如采样点数是2000,采样频率是1000Hz,那么真实频率就是0.5,1,1.5,…,500Hz。
由此可知道,采样率不变,采样点数越多,频率分辨率便越高。
fft之后得到的最高频率分量是采样频率的1/2,频率区分的刻度(即频率分辨力=采样频率/采样点数)。
假设真实频率为
f
(
k
)
f(k)
f(k),采样频率为
f
s
f_s
fs,采样点数为
N
N
N,真实频率为:
f
(
k
)
=
k
∗
(
f
s
N
)
f(k)=k*(\frac{f_s}{N})
f(k)=k∗(Nfs)
仿真实验
以下仿真程序,信号模型为多频点信号,采样率为1000,采样点数为2000个,因此频率分辨率为1000/2000=0.5Hz,频率值为0.5,1,1,5,…500。
代码如下:
%% 离散傅里叶变换
% 2022.5.26
% cleal all; close all; clc;
T = 2; % 总时间
fs = 1000; % 采样率
dt = 1/fs; % 采样间隔
t = 0:dt:T-dt; % 时间刻度向量
N = length(t); %序列长度
f1 = 100; % 频率分量1
f2 = 300; % 频率分量2
x = cos(2*pi*f1*t)+cos(2*pi*f2*t); % 基带信号
plot(t,x) % 绘图
X = fftshift(fft(x))*2/N; % fft
deltFrequency = fs/N; %频率分辨率
freqAxis = -fs/2:deltFrequency:fs/2-deltFrequency; % 频率正半轴
stem(freqAxis,abs(X));% 绘制幅度谱
axis([-501,501,0,1.3]); %控制坐标轴变量
幅度谱如下:
在进行fft并且绘图的时候,经常需要进行频率轴坐标的设置,有点繁琐,所以我编写了一个函数,可以将fft和绘制频谱图封装起来:
function [X] = fftPlot(x,fs)
% 该函数用于离散傅里叶变换,并且绘制频谱图,重点在于频谱图的下标显示为对应频率值
% Author:huasir 2023.11.22 @Beijing
% Input :
% * x: 输入信号
% * fs:采样率
% Output :
% * X:离散傅里叶变换的结果,复数
N = length(x);
X = fftshift(fft(x))*2/N; % fft
deltFrequency = fs/N; %频率分辨率
freqAxis = -fs/2:deltFrequency:fs/2-deltFrequency; % 频率正半轴
figure;
plot(freqAxis,abs(X));% 绘制幅度谱
end