快速傅立叶变换补0

简介

转载于:快速傅里叶变换(FFT)中为什么要“补零”?内积,卷积,傅里叶变换

说明一下,用MATLAB做FFT并不要求数据点个数必须为以2为基数的整数次方。之所以很多资料上说控制数据点数为以2为基数的整数次 方,是因为这样就能采用以2为基的FFT算法,提升运算性能。

如果数据点数不是以2为基数的整数次方,处理方法有两种,一种是在原始数据开头或末尾补零,即将数据补到以2为基数的整数次方,这是“补零”的一个用处;第二种是采用以任意数为基数的FFT算法。

MATLAB的FFT并非使用的是蝶形算法。对于非2的幂指数的长度,它使用Cooley-Tukey算法分解因数,再利用现有的库分别计算最后合成。

快速傅里叶变换FFT

比如,现在我有一个信号,这个信号中仅包含两个正(余)弦波,一个是 $1MHz$,一个是 $1.05MHz$,即 $x=\cos(2\pi\times1000000t)+\cos(2\pi\times1050000t)$。设定采样频率为 $F_s=100MHz$,如果采 $1000$ 个点,那么时域信号的时长就有 $10\mu s$。

图1. 1000个数据点

如果,直接对这1000个数据点其做快速傅里叶变换,将得到频谱图:

图2. 1000个数据点做FFT的频谱

可以发现,频谱点稀疏,在1MHz附近根本无法将1MHz 和1.05MHz 的两个频率分开。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
clear;clc
close all
%% FFT
Fs = 100e6; % 采样频率 Hz
T = 1/Fs; % 采样周期 s
L0 = 1000; % 信号长度
L = 1000; % 数据长度
t0 = (0:L0-1)*T; % 信号时间序列
x = cos(2*pi*1e6*t0) + cos(2*pi*1.05e6*t0); % 信号函数
t = (0:L-1)*T; % 数据时间序列
%% Plot
figure(1)
plot(t*1e6,x,'b-','linewidth',1.5)
title ('\fontsize{10}\fontname{Times New Roman}Time domain signal')
xlabel('\fontsize{10}\fontname{Times New Roman}\it t /\rm \mus')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itt\rm)')
grid on;
axis([0 10 -2 2])
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
%% FFT
Y = fft(x); % 快速傅里叶变换
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L0);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;% Rfft
figure(2)
plot(f, P1,'r-','Marker','.','markersize',10,'linewidth',1.5)
axis([0.5e6 1.5e6 0 1.5])
title('\fontsize{10}\fontname{Times New Roman}Power Spectrum')
xlabel('\fontsize{10}\fontname{Times New Roman}\it f /\rm Hz')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itf\rm)')
grid on;
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')

频率分辨率

发现频率成分无法被区分开来,第一反应应该就是:频率分辨率不够。那么,如何提高频率分辨率呢?首先要清楚,这里存在两种类型的频率分辨率。

一种叫波形分辨率,其由原始数据的时间长度决定:

另一种可以称之为视觉分辨率FFT分辨率,其由采样频率和参与FFT的数据点数决定:

之所以要区分,就是因为后面要进行“补零”的操作。如果不补零,直接对原始数据做FFT,那么这两种分辨率是相等的。

例如上面,有:

补零

那么,如果现在在原始数据点后补零会有什么效果呢?假设在这1000个原始数据点后面再补充零达到7000个点,那么数据变成了:

图3. 7000个补零后数据点

此时对其做快速傅里叶变换,结果如下:

图4. 7000个补零后数据点做FFT的频谱

可以发现,频谱点密集了不少,但是在$1MHz$附近依然无法将$1MHz$和$1.05MHz$的两个频率成分分开。这是因为从式(1)可以看出,波形分辨率只与原始数据的时长 $T$ 有关,而与参与FFT的数据点数无关。所以,虽然补了很多零,但波形分辨率依然为:$\Delta R_{w}=\frac{1}{10 \mu s}=100 k H z$,该分辨率大于$1MHz$和$1.05MHz$成分之间的距离 $50kHz$。这就好比用筛子分黄豆和大米,分辨率就好像是筛子上孔的大小,如果筛子的孔太大了,就没有办法把这两者分开。

而“时域补零相当于频域插值”,也就是说,补零操作增加了频域的插值点数,让频域曲线看起来更加光滑,也就是增加了FFT频率分辨率,注意式(2)所示,这是“补零”的另一个原因

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
clear;clc
close all
%% FFT
Fs = 100e6; % 采样频率 Hz
T = 1/Fs; % 采样周期 s
L0 = 1000; % 信号长度
L = 7000; % 数据长度
t0 = (0:L0-1)*T; % 信号时间序列
x = cos(2*pi*1e6*t0) + cos(2*pi*1.05e6*t0); % 信号函数
t = (0:L-1)*T; % 数据时间序列
y = zeros(1,L);
y(1:L0) = x;
%% Plot
figure(1)
plot(t*1e6,y,'b-','linewidth',1.5)
title('\fontsize{10}\fontname{Times New Roman}Time domain signal with Zero Padding')
xlabel('\fontsize{10}\fontname{Times New Roman}\it t /\rm \mus')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itt\rm)')
grid on;
axis([0 70 -2 2])
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
%% FFT
Y = fft(y); % 快速傅里叶变换
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L0);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;% Rfft
figure(2)
plot(f, P1,'r-','Marker','.','markersize',10,'linewidth',1.5)
axis([0.5e6 1.5e6 0 1.5])
title('\fontsize{10}\fontname{Times New Roman}Power Spectrum')
xlabel('\fontsize{10}\fontname{Times New Roman}\it f /\rm Hz')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itf\rm)')
grid on;
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')

频谱泄露

显然,根据上面的分析可知,在采样频率不变的情况下,要想将$1MHz$和$1.05MHz$这两个频率成分分析出来,光靠“补零”是不够的,必须要改变波形分辨率,也就是要延长原始数据的时长。现在以相同的采样频率对信号采7000个点作为原始信号:

图5. 7000个数据点

对其做快速傅里叶变换,结果如下:

图6. 7000个数据点做FFT的频谱

因为此时的波形分辨率为:$\Delta R_{w}=\frac{1}{70 \mu s} \approx 14 k H z$,小于$1MHz$和$1.05MHz$这两个频率成分之间的距离 $50kHz$ ,所以可以看出有两个明显的峰值。

但是会发现$1MHz$ 对应的幅值为1,与原始信号中该频率成分的幅值一致,但 $1.05MHz$对应的幅值明显低于1,但是其周边的点上却都有不小的幅值,这就是所谓的频谱泄露,因为数据点的个数影响,使得在$1MHz$ 处有谱线存在,但在$1.05MHz$ 处没有谱线存在,使测量结果偏离实际值 ,同时在实际频率点的能量分散到两侧的其它频率点上,并出现一些幅值较小的假谱。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
clear;clc
close all
%% FFT
Fs = 100e6; % 采样频率 Hz
T = 1/Fs; % 采样周期 s
L0 = 7000; % 信号长度
L = 7000; % 数据长度
t0 = (0:L0-1)*T; % 信号时间序列
x = cos(2*pi*1e6*t0) + cos(2*pi*1.05e6*t0); % 信号函数
t = (0:L-1)*T; % 数据时间序列
%% Plot
figure(1)
plot(t*1e6,x,'b-','linewidth',1.5)
title('\fontsize{10}\fontname{Times New Roman}Time domain signal')
xlabel('\fontsize{10}\fontname{Times New Roman}\it t /\rm \mus')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itt\rm)')
grid on;
axis([0 70 -2 2])
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
%% FFT
Y = fft(x); % 快速傅里叶变换
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L0);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs/2*linspace(0,1,L/2+1);% Rfft
figure(2)
plot(f, P1,'r-','Marker','.','markersize',10,'linewidth',1.5)
axis([0.5e6 1.5e6 0 1.5])
title('\fontsize{10}\fontname{Times New Roman}Power Spectrum')
xlabel('\fontsize{10}\fontname{Times New Roman}\it f /\rm Hz')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itf\rm)')
grid on;
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')

为了解决这个问题,可以设法使得谱线同时经过$1MHz$和$1.05MHz$这两个频率点,找到他们的公约数。

如果原始数据不变,在后面再补充1000个零点:

图7. 8000个补零后数据点

那么FFT分辨率就是 $12.5kHz$ ,是这两个频率的公约数 $1 M H z=80 \times 12.5 k H z$ ,$1.05 M H z=84 \times 12.5 k H z$,所以谱线同时经过 $1MHz$和$1.05MHz$ 这两个频率点。

对其做快速傅里叶变换,结果如下:

图8. 8000个补零后数据点做FFT的频谱

会发现$1MHz$和$1.05MHz$对应的幅值均为1,与原始信号一致。这也是一种补零操作带来的影响

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
clear;clc
close all
%% FFT
Fs = 100e6; % 采样频率 Hz
T = 1/Fs; % 采样周期 s
L0 = 7000; % 信号长度
L = 8000; % 数据长度
t0 = (0:L0-1)*T; % 信号时间序列
x = cos(2*pi*1e6*t0) + cos(2*pi*1.05e6*t0); % 信号函数
t = (0:L-1)*T; % 数据时间序列
y = zeros(1,L);
y(1:L0) = x;
%% Plot
figure(1)
plot(t*1e6,y,'b-','linewidth',1.5)
title('\fontsize{10}\fontname{Times New Roman}Time domain signal')
xlabel('\fontsize{10}\fontname{Times New Roman}\it t /\rm \mus')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itt\rm)')
grid on;
axis([0 80 -2 2])
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
%% FFT
Y = fft(y); % 快速傅里叶变换
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L0);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs/2*linspace(0,1,L/2+1);% Rfft
figure(2)
plot(f, P1,'r-','Marker','.','markersize',10,'linewidth',1.5)
axis([0.5e6 1.5e6 0 1.5])
title('\fontsize{10}\fontname{Times New Roman}Power Spectrum')
xlabel('\fontsize{10}\fontname{Times New Roman}\it f /\rm Hz')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itf\rm)')
grid on;
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')

图8中会有一些旁瓣出现,这是因为补零影响了原始信号,如果,直接采8000个点作为原始数据,即将程序中的L0改为8000,那么有:

图9. 8000个数据点

并对其做FFT,结果如下:

图10. 8000个数据点做FFT的频谱

这样就不存在补零带来的误差了。

FFT与卷积

时域卷积相当于FFT乘积。卷积的离散定义:

参考

  1. Zero Padding http://www.bitweenie.com/listings/fft-zero-padding/
  2. Zero Padding Theorem https://ccrma.stanford.edu/~jos/dft/Zero_Padding_Theorem_Spectral.html
一分一毛也是心意