一、傅里叶变换
傅里叶变换是一种全局的变换,时域信号经过傅里叶变换后,就变成了频域信号,所以从傅里叶变换频谱上是无法看出时域信息的。
傅里叶变换只适合处理平稳信号,对于非平稳信号,由于频率会随时间变化,为了捕获这种时变特性,需要对信号进行时频分析,比如:短时傅里叶变换、小波变换、希尔伯特黄变换等。短时傅里叶变换详见:【20211206】【信号处理】时频分析 —— 短时傅里叶变换(STFT)
二、举个栗子
使用 Matlab 生成一个单频信号,并进行 FFT,作图观察。
%% 参数设置
fs = 50; % 采样频率
L = 10; % 时间长度
t = 0 : 1/fs : L; % 时间坐标
%% 生成信号
fsignal = 10; % 信号频率
A = 10; % 信号幅值
s = A * sin(2 * pi * fsignal * t);
%% FFT
numfft = 256; % fft点数
s_fft = fft(s, numfft); % fft
p = abs(s_fft(1:numfft/2)) .^ 2; % 能量(只取正频部分)
f = (0 : numfft/2-1) / numfft * fs; % 频率(只取正频部分)
%% 作图
figure(1);
subplot(2, 1, 1);
plot(s); axis('tight'); title('信号s');
subplot(2, 1, 2);
plot(f, p); axis('tight'); title('信号s的傅里叶变换');
对运行结果的解读:生成一个频率为 10Hz 的正弦信号,对信号的采样频率为 50Hz,即 1s 的时间段内抽 50 个采样点,对该时域信号进行 FFT,只取正频部分的信号并做图。使用 abs 的平方计算能量,使用采样频率 fs 计算 FFT 图的横坐标轴坐标。得到的图形在 10Hz 处有一个尖峰,表明该处能量很高,这是因为该信号是一个 10Hz 的单频信号,所以频谱只在此处有值。
三、含有直流分量的信号做FFT问题
%%
clear; clc; close all; warning off;
%% 参数设置
fs = 50; % 采样频率
L = 10; % 时间长度
t = 0 : 1/fs : L; % 时间坐标
%% 生成信号
fsignal = 10; % 信号频率
s1 = 10 * sin(2 * pi * fsignal * t);
s2 = 10 * ones(1, length(t));
s = s1 + s2;
sDC = s - mean(s); % 信号s去直流
%% FFT
numfft = 2048; % fft点数
s_fft = fft(s, numfft); % fft
p = abs(s_fft(1:numfft/2)) .^ 2; % 能量(只取正频部分)
f = (0 : numfft/2-1) / numfft * fs; % 频率(只取正频部分)
sDC_fft = fft(sDC, numfft);
p2 = abs(sDC_fft(1:numfft/2)) .^ 2; % 能量(只取正频部分)
f2 = (0 : numfft/2-1) / numfft * fs; % 频率(只取正频部分)
%% 作图
figure(1);
subplot(3, 1, 1);
plot(s, 'linewidth', 1.5); axis('tight'); title('信号s');
subplot(3, 1, 2);
plot(f, p, 'linewidth', 1.5); axis('tight'); title('信号s的傅里叶变换');
subplot(3, 1, 3);
plot(f2, p2, 'linewidth', 1.5); axis('tight'); title('去直流信号sDC的傅里叶变换');
可以发现,如果信号含有直流分量,那么体现在频谱上,在 f=0 处会有一个能量尖峰,这是因为直流分量就是频率为 0 的信号。如果不进行去直流处理,那么会影响根据 FFT 频谱提取信号基频,因此要进行去直流处理!
四、知识点
1. 区分采样频率和信号频率
采样频率:是指 1s 内采多少个采样点,是抽样的概念。
信号频率:是指信号周期的倒数,是信号的概念。
进行 FFT 变换时,计算横坐标频率坐标值时使用的是采样频率 fs,从频谱分析到的能量尖峰对应的横坐标位置是信号频率。
2. 区分时间分辨率和频率分辨率
时间分辨率和频率分辨率时时频分析中很重要的两个指标,若有采样频率为 ,FFT 点数为 ,则:
时间分辨率:
频率分辨率:
(参考:关于频谱分析中两个重要指标:频率分辨率和时间分辨率的理解及计算)
3. FFT 点数的选取
FFT 点数的选取和要求分析的频率分辨率有关系,做完 FFT 之后纵坐标是 采样频率/2 被平均划分到 N/2 个点,所以 N 越大,FFT 之后的频谱频率分辨率越高。
频率分辨率 = 采样频率 / FFT点数
(1)例一
%% 参数设置
fs = 50; % 采样频率
L = 10; % 时间长度
t = 0 : 1/fs : L; % 时间坐标
%% 生成信号
fsignal = 10; % 信号频率
A = 10; % 信号幅值
s = A * sin(2 * pi * fsignal * t);
%% FFT
numfft = 64; % fft点数
s_fft = fft(s, numfft); % fft
p = abs(s_fft(1:numfft/2)) .^ 2; % 能量(只取正频部分)
f = (0 : numfft/2-1) / numfft * fs; % 频率(只取正频部分)
numfft2 = 1024; % fft点数
s_fft2 = fft(s, numfft2); % fft
p2 = abs(s_fft2(1:numfft2/2)) .^ 2; % 能量(只取正频部分)
f2 = (0 : numfft2/2-1) / numfft2 * fs; % 频率(只取正频部分)
%% 作图
figure(1);
subplot(3, 1, 1);
plot(s); axis('tight'); title('信号s');
subplot(3, 1, 2);
plot(f, p); axis('tight'); title('64点傅里叶变换');
subplot(3, 1, 3);
plot(f2, p2); axis('tight'); title('1024点傅里叶变换');
(2)例二
%% 观察多种采样点前提下的频谱变化
clear; clc; close all; warning off;
%% 参数设置
fs = 100; % 采样频率
L = 10; % 时间长度
t = 0 : 1/fs : L; % 时间坐标
%% 生成信号
fsignal = 20; % 信号频率
s1 = 10 * sin(2 * pi * fsignal * t);
s2 = 10 * ones(1, length(t));
s = s1 + s2;
sDC = s - mean(s); % 信号s去直流
%% FFT
numfft = 32;
numfft1 = 64; % fft点数
numfft2 = 128;
numfft3 = 256;
numfft4 = 512;
sDC_fft = fft(sDC, numfft);
p = abs(sDC_fft(1:numfft/2)) .^ 2; % 能量(只取正频部分)
f = (0 : numfft/2-1) / numfft * fs; % 频率(只取正频部分)
sDC_fft1 = fft(sDC, numfft1);
p1 = abs(sDC_fft1(1:numfft1/2)) .^ 2; % 能量(只取正频部分)
f1 = (0 : numfft1/2-1) / numfft1 * fs; % 频率(只取正频部分)
sDC_fft2 = fft(sDC, numfft2);
p2 = abs(sDC_fft2(1:numfft2/2)) .^ 2; % 能量(只取正频部分)
f2 = (0 : numfft2/2-1) / numfft2 * fs; % 频率(只取正频部分)
sDC_fft3 = fft(sDC, numfft3);
p3 = abs(sDC_fft3(1:numfft3/2)) .^ 2; % 能量(只取正频部分)
f3 = (0 : numfft3/2-1) / numfft3 * fs; % 频率(只取正频部分)
sDC_fft4 = fft(sDC, numfft4);
p4 = abs(sDC_fft4(1:numfft4/2)) .^ 2; % 能量(只取正频部分)
f4 = (0 : numfft4/2-1) / numfft4 * fs; % 频率(只取正频部分)
%% 作图
figure(1);
subplot(5, 1, 1);
plot(f, p, 'linewidth', 1.5); axis('tight'); title('32点傅里叶变换');
subplot(5, 1, 2);
plot(f1, p1, 'linewidth', 1.5); axis('tight'); title('64点傅里叶变换');
subplot(5, 1, 3);
plot(f2, p2, 'linewidth', 1.5); axis('tight'); title('128点傅里叶变换');
subplot(5, 1, 4);
plot(f3, p3, 'linewidth', 1.5); axis('tight'); title('256点傅里叶变换');
subplot(5, 1, 5);
plot(f4, p4, 'linewidth', 1.5); axis('tight'); title('512点傅里叶变换');
(参考:数字信号处理】Matlab做fft时点数N怎么选取)
(参考:如何决定要使用多少点来做fft)
(参考:用MATLAB演示采样频率和点数对FFT的影响-好东西一定要及时保存)
(参考:快速傅里叶变换FFT点数和采样率)