一、DTFT、DFT、FFT、DFS的关系
1、DTFT与DFT的关系
在实际应用中,计算机只能处理离散信号,所以对连续信号
x
(
t
)
x(t)
x(t)进行时域采样,得到一组离散样本
x
(
n
)
x(n)
x(n),对它进行傅里叶变换得到
X
(
w
)
=
∑
n
=
−
∞
∞
x
(
n
)
e
−
j
w
n
X(w)=\sum_{n=-∞}^{∞} x(n)e^{-jwn}
X(w)=n=−∞∑∞x(n)e−jwn
上式即为离散时间傅里叶变换(DTFT),由于变换后得到的频域值仍然是连续的,继续对频域进行采样,得到
X
(
k
)
=
∑
n
=
0
N
−
1
x
(
n
)
e
−
j
2
π
k
n
N
X(k)=\sum_{n=0}^{N-1} x(n)e^{-j\frac{2\pi k n}{N}}
X(k)=n=0∑N−1x(n)e−jN2πkn
上式就是离散傅里叶变换(DFT)。目前计算机中常用的快速傅里叶变换(FFT),就是DFT的快速算法。
也就是说:DFT是DTFT的离散化,DTFT在频域是连续的,在归一化ω轴上是连续的,以2π为周期。
2、DFS与DFT的关系
离散傅里叶级数(DFS)是针对周期序列的,即无限长序列;
而将DFS取出一个周期就是DFT,即DFT是有限长序列;
实际上,DFT就是DFS在一个周期内的取值,一个序列的DFT实际上就是这个序列的频谱在一个周期内等间隔采样的样点值。
二、DFT的实质
离散傅里叶变换:DFT(Discrete Fourier Transform)
其实质是有限长序列傅里叶变换的有限点离散采样,从而实现了频域离散化,使数字信号处理可以在频域采用数值运算的方法进行,这样就大大增加了数字信号处理的灵活性。
更重要的是,DFT有多种快速算法,统称为快速傅里叶变换(Fast Fourier Transform,FFT),从而使信号的实时处理和设备的简化得以实现。
三、DFT的定义
设
x
(
n
)
x(n)
x(n)是一个长度为
M
M
M的有限长序列,则定义
x
(
n
)
x(n)
x(n)的
N
N
N点离散傅里叶变换为
X
(
k
)
=
D
F
T
[
x
(
n
)
]
=
∑
n
=
0
N
−
1
x
(
n
)
W
N
k
n
k
=
0
,
1
,
⋯
,
N
−
1
X(k)=DFT[x(n)]=\sum_{n=0}^{N-1}x(n)W_{N}^{kn} \qquad k=0,1,\cdots,N-1
X(k)=DFT[x(n)]=n=0∑N−1x(n)WNknk=0,1,⋯,N−1
X
(
k
)
X(k)
X(k)的离散傅里叶逆变换(Inverse Discrete Fourier Transfrom,IDFT)为
x
(
n
)
=
I
D
F
T
[
X
(
k
)
]
=
1
N
∑
k
=
0
N
−
1
X
(
k
)
W
N
−
k
n
n
=
0
,
1
,
⋯
,
N
−
1
x(n)=IDFT[X(k)]={1\over N}\sum_{k=0}^{N-1}X(k)W_{N}^{-kn} \qquad n=0,1,\cdots,N-1
x(n)=IDFT[X(k)]=N1k=0∑N−1X(k)WN−knn=0,1,⋯,N−1
式中,
W
N
=
e
−
j
2
π
N
W_{N}=e^{-j \frac{2\pi}{N}}
WN=e−jN2π,
N
N
N称为
D
F
T
DFT
DFT变换区间长度,
N
≥
M
N\ge M
N≥M。
四、DFT的物理意义
DFT与傅里叶变换和Z变换的关系:
设序列
x
(
n
)
x(n)
x(n)的长度为
M
M
M,其
Z
Z
Z变换和
N
(
N
≥
M
)
N(N\ge M)
N(N≥M)点DFT分别为
X
(
z
)
=
Z
T
[
x
(
n
)
]
=
∑
n
=
0
M
−
1
x
(
n
)
z
−
n
X(z)=ZT[x(n)]=\sum_{n=0}^{M-1}x(n)z^{-n}
X(z)=ZT[x(n)]=n=0∑M−1x(n)z−n
X
(
k
)
=
D
F
T
[
x
(
n
)
]
N
=
∑
n
=
0
M
−
1
x
(
n
)
W
N
k
n
k
=
0
,
1
,
⋯
,
N
−
1
X(k)=DFT[x(n)]_N=\sum_{n=0}^{M-1}x(n)W_{N}^{kn} \qquad k=0,1,\cdots,N-1
X(k)=DFT[x(n)]N=n=0∑M−1x(n)WNknk=0,1,⋯,N−1
比较上面二式可得关系式
X
(
k
)
=
X
(
z
)
∣
z
=
e
j
2
π
N
k
k
=
0
,
1
,
⋯
,
N
−
1
(1)
X(k)=X(z)\Big|_{z=e^{j\frac{2\pi}{N}k}} \qquad k=0,1,\cdots,N-1 \tag{1}
X(k)=X(z)∣∣∣z=ejN2πkk=0,1,⋯,N−1(1)
或
X
(
k
)
=
X
(
e
j
ω
)
∣
ω
=
2
π
N
k
k
=
0
,
1
,
⋯
,
N
−
1
(2)
X(k)=X(e^{j\omega})\Big|_{\omega=\frac{2\pi}{N}k} \qquad k=0,1,\cdots,N-1 \tag{2}
X(k)=X(ejω)∣∣∣ω=N2πkk=0,1,⋯,N−1(2)
(1)式表明序列
x
(
n
)
x(n)
x(n)的
N
N
N点
D
F
T
DFT
DFT是
x
(
n
)
x(n)
x(n)的
Z
Z
Z变换在单位圆上的
N
N
N点等间隔采样。(2)式则说明
X
(
k
)
X(k)
X(k)为
x
(
n
)
x(n)
x(n)的傅里叶变换
X
(
e
j
ω
)
X(e^{j\omega})
X(ejω)在区间
[
0
,
2
π
]
[0,2\pi]
[0,2π]上的
N
N
N点等间隔采样。这就是DFT的物理意义。
由此可见,DFT的变化区间长度N不同,表示对 X ( e j ω ) X(e^{j\omega}) X(ejω)在区间 [ 0 , 2 π ] [0,2\pi] [0,2π]上的采样间隔和采样点数不同,所以DFT的变换结果不同。
xn=boxcar(4);
M=1024; % DFT变换区间长度
Xw=fft(xn,M); % fft函数在xn后面补零,返回xn的M点DFT变换结果向量
wk=2*(0:M-1)/M;
X8k=fft(xn,8); % 计算xn的8点DFT
X16k=fft(xn,16); % 计算xn的16点DFT
subplot(3,2,1);
plot(wk,abs(Xw));
title('(a)x(n)的幅频特性曲线 ');xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');axis([0,2,0,4.5]);
subplot(3,2,3);
k=0:7; stem(k,abs(X8k),'.');
title('(b)x(n)的8点DFT');xlabel('k');ylabel('|X(k)|');axis([0,8,0,4.5]);
subplot(3,2,5);
k=0:15; stem(k,abs(X16k),'.')
title('(c)x(n)的16点DFT');xlabel('k');ylabel('|X(k)|');axis([0,16,0,4.5]);
%%%%%%%%%% DFT的MATLAB计算示例 %%%%%%%%%%
xn=[1 1 1 1]; %输入长度M=4的时域序列向量
Xk16=fft(xn,16); %计算xn的16点DFT
Xk32=fft(xn,32); %计算xn的32点DFT
k=0:15;
wk=2*k/16; %产生16点DFT对应的采样点频率(关于π归一化值)
subplot(2,2,1);
stem(wk,abs(Xk16),'.'); %绘制16点DFT的幅频特性图
title('(a)16点DFT的幅频特性图');
xlabel('ω/π');
ylabel('幅度');
subplot(2,2,3);
stem(wk,angle(Xk16),'.'); %绘制16点DFT的相频特性图
line([0,2],[0,0]); %绘制一条从(0,0)到(2,0)的线条
title('(b)16点DFT的相频特性图')
xlabel('ω/π');
ylabel('相位');
axis([0,2,-3.5,3.5]); %更改坐标轴范围,使 x 轴的范围从0到2,y轴的范围从-3.5 到3.5
k=0:31;
wk=2*k/32; %产生32点DFT对应的采样点频率(关于π归一化值)
subplot(2,2,2);
stem(wk,abs(Xk32),'.'); %绘制32点DFT的幅频特性图
title('(c)32点DFT的幅频特性图');
xlabel('ω/π');
ylabel('幅度');
subplot(2,2,4);
stem(wk,angle(Xk32),'.'); %绘制32点DFT的相频特性图
line([0,2],[0,0]);
title('(d)32点DFT的相频特性图');
xlabel('ω/π');
ylabel('相位');
axis([0,2,-3.5,3.5]);
五、用DFT对信号进行谱分析
所谓信号的谱分析,就是计算信号的傅里叶变换。连续信号与系统的傅里叶分析不便于直接用计算机进行计算,使其应用受到限制。而DFT是一种时域和频域均离散化的变换,适合数值运算,成为用计算机分析离散信号和系统的有力工具。对连续信号和系统,可以通过时域采样,应用DFT进行近似谱分析。
1、用DFT对连续信号进行谱分析
1.1 实际应用DFT所面临的问题
工程实际中,经常遇到连续信号 x a ( t ) x_a(t) xa(t),其频谱函数 X a ( j Ω ) X_a(jΩ) Xa(jΩ)也是连续函数。为了利用DFT对 x a ( t ) x_a(t) xa(t)进行频谱分析,先对 x a ( t ) x_a(t) xa(t)进行时域采样,得到 x ( n ) = x a ( n T ) x(n)=x_a(nT) x(n)=xa(nT),再对 x ( n ) x(n) x(n)进行DFT,得到的 X ( k ) X(k) X(k)则是 x ( n ) x(n) x(n)的傅里叶变换 X ( e j ω ) X(e^{j\omega}) X(ejω)在频率区间 [ 0 , 2 π ] [0,2\pi] [0,2π]上的 N N N点等间隔采样。这里 x ( n ) x(n) x(n)和 X ( k ) X(k) X(k)均为有限长序列。
然而,由傅里叶变换理论知道,若信号持续时间有限长,则其频谱无限宽;若信号的频谱有限宽,则其持续时间必然为无限长。所以严格地讲,持续时间有限的带限信号是不存在的。因此,按采样定理采样时,上述两种情况下的采样序列 x ( n ) = x a ( n T ) x(n)=x_a(nT) x(n)=xa(nT)均应为无限长,不满足DFT的变换条件。
实际上对频谱很宽的信号,为防止时域采样后产生频谱混叠失真,可用预滤波器滤除幅度较小的高频成分,使连续信号的带宽小于折叠频率。对于持续时间很长的信号,采样点数太多,以致无法存储和计算,只好截取有限点进行DFT。
由上述可见,用DFT对连续信号进行频谱分析必然是近似的,其近似程度与信号带宽、采样频率和截取长度有关。实际上从工程角度看,滤除幅度很小的高频成分和截取幅度很小的部分时间信号是允许的。因此,在下面分析中,假设 x a ( t ) x_a(t) xa(t)是经过预滤波和截取处理的有限长带限信号。
1.2 对实际连续信号的预处理及假设
设连续信号 x a ( t ) x_a(t) xa(t)持续时间为 T p T_p Tp,最高频率为 f c f_c fc。 x a ( t ) x_a(t) xa(t)的傅里叶变换为 X a ( j Ω ) X_a(jΩ) Xa(jΩ),对 x a ( t ) x_a(t) xa(t)进行时域采样得到 x ( n ) = x a ( n T ) x(n)=x_a(nT) x(n)=xa(nT), x ( n ) x(n) x(n)的傅里叶变换为 X ( e j ω ) X(e^{j\omega}) X(ejω)。
由假设条件可知
x
(
n
)
x(n)
x(n)的长度为
N
=
T
p
T
=
T
p
F
s
N=\frac{T_p}{T}=T_p F_s
N=TTp=TpFs
式中,
T
T
T为采样间隔,
F
s
=
1
T
F_s=\frac{1}{T}
Fs=T1为采样频率。
用 X ( k ) X(k) X(k)表示 x ( n ) x(n) x(n)的N点DFT,下面推导出 X ( k ) X(k) X(k)与 X a ( j Ω ) X_a(jΩ) Xa(jΩ)的关系,最后由此关系归纳出用 X ( k ) X(k) X(k)表示 X a ( j Ω ) X_a(jΩ) Xa(jΩ)的方法,即用DFT对连续信号进行谱分析的方法。
1.3 用DFT表示连续信号傅里叶变换的推导过程
时域离散信号的傅里叶变换和模拟信号傅里叶变换之间的关系为:
X
(
e
j
ω
)
=
X
^
(
j
Ω
)
∣
Ω
=
ω
T
=
1
T
∑
k
=
−
∞
∞
X
a
(
j
ω
−
2
π
k
T
)
X(e^{j\omega})=\hat X(j\Omega)\Big|_{\Omega=\frac{\omega}{T}}=\frac{1}{T}\sum_{k=-∞}^{∞}X_a(j\frac{\omega-2\pi k}{T})
X(ejω)=X^(jΩ)∣∣∣Ω=Tω=T1k=−∞∑∞Xa(jTω−2πk)
由上式知道,
x
(
n
)
x(n)
x(n)的傅里叶变换
X
(
e
j
ω
)
X(e^{j\omega})
X(ejω)与
x
a
(
t
)
x_a(t)
xa(t)的傅里叶变换
X
a
(
j
Ω
)
X_a(jΩ)
Xa(jΩ)满足如下关系:
X
(
e
j
ω
)
=
1
T
∑
m
=
−
∞
∞
X
a
[
j
(
ω
T
−
2
π
T
m
)
]
X(e^{j\omega})=\frac{1}{T}\sum_{m=-∞}^{∞}X_a\Big[j({\omega \over T}-\frac{2\pi}{T}m)\Big]
X(ejω)=T1m=−∞∑∞Xa[j(Tω−T2πm)]
将
ω
=
Ω
T
\omega=\Omega T
ω=ΩT代入上式,得到:
X
(
e
j
Ω
T
)
=
1
T
∑
m
=
−
∞
∞
X
a
[
j
(
Ω
−
2
π
T
m
)
]
=
1
T
X
~
a
(
j
Ω
)
(3)
X(e^{j\Omega T})=\frac{1}{T}\sum_{m=-∞}^{∞}X_a\Big[j(\Omega-\frac{2\pi}{T}m)\Big]=\frac{1}{T}\tilde X_a(j\Omega) \tag{3}
X(ejΩT)=T1m=−∞∑∞Xa[j(Ω−T2πm)]=T1X~a(jΩ)(3)
式中
X
~
a
(
j
Ω
)
=
∑
m
=
−
∞
∞
X
a
[
j
(
Ω
−
2
π
T
m
)
]
\tilde X_a(j\Omega)=\sum_{m=-∞}^{∞}X_a\Big[j(\Omega-\frac{2\pi}{T}m)\Big]
X~a(jΩ)=m=−∞∑∞Xa[j(Ω−T2πm)]
表示模拟信号频谱
X
a
(
j
Ω
)
X_a(j\Omega)
Xa(jΩ)的周期延拓函数。
由
x
(
n
)
x(n)
x(n)的N点DFT的定义有
X
(
k
)
=
D
F
T
[
x
(
n
)
]
N
=
X
(
e
j
ω
)
∣
ω
=
2
π
N
k
k
=
0
,
1
,
⋯
,
N
−
1
(4)
X(k)=DFT[x(n)]_N=X(e^{j\omega})\Big|_{\omega=\frac{2\pi}{N}k} \qquad k=0,1,\cdots,N-1 \tag{4}
X(k)=DFT[x(n)]N=X(ejω)∣∣∣ω=N2πkk=0,1,⋯,N−1(4)
将(4)式代入(3)式,得到:
X
(
k
)
=
X
(
e
j
2
π
N
k
)
=
1
T
X
~
a
(
j
2
π
N
T
k
)
=
1
T
X
~
a
(
j
2
π
T
p
k
)
k
=
0
,
1
,
⋯
,
N
−
1
(5)
X(k)=X(e^{j\frac{2\pi}{N}k})=\frac{1}{T}\tilde X_a(j\frac{2\pi}{NT}k)=\frac{1}{T}\tilde X_a(j\frac{2\pi}{T_p}k) \qquad k=0,1,\cdots,N-1 \tag{5}
X(k)=X(ejN2πk)=T1X~a(jNT2πk)=T1X~a(jTp2πk)k=0,1,⋯,N−1(5)
(5)式说明了
X
(
k
)
X(k)
X(k)与
X
a
(
j
Ω
)
X_a(jΩ)
Xa(jΩ)的关系。
1.4 将以自变量ω表示的关系转换为以频率f为自变量
为了符合一般的频谱描述习惯,以频率
f
f
f为自变量,整理(5)式。
令
{
X
a
′
(
f
)
=
X
a
(
j
Ω
)
∣
Ω
=
2
π
f
=
X
a
(
j
2
π
f
)
X
~
a
′
(
f
)
=
X
~
a
(
Ω
)
∣
Ω
=
2
π
f
=
X
~
a
(
2
π
f
)
\begin{cases} X'_a(f)=X_a(j \Omega)\Big|_{ \Omega=2\pi f}=X_a(j2\pi f )\\ \tilde X'_a(f)=\tilde X_a(\Omega)\Big|_{ \Omega=2\pi f}=\tilde X_a(2\pi f ) \end{cases}
⎩⎪⎨⎪⎧Xa′(f)=Xa(jΩ)∣∣∣Ω=2πf=Xa(j2πf)X~a′(f)=X~a(Ω)∣∣∣Ω=2πf=X~a(2πf)
则(5)式变为
X
(
k
)
=
1
T
X
~
a
′
(
f
)
∣
f
=
k
T
p
=
1
T
X
~
a
′
(
k
F
)
k
=
0
,
1
,
⋯
,
N
−
1
(5)
X(k)=\frac{1}{T}\tilde X'_a(f)\Big|_{f=\frac{k}{T_p}}=\frac{1}{T}\tilde X'_a(kF) \qquad k=0,1,\cdots,N-1 \tag{5}
X(k)=T1X~a′(f)∣∣∣f=Tpk=T1X~a′(kF)k=0,1,⋯,N−1(5)
由此可得
X
~
a
′
(
k
F
)
=
T
X
(
k
)
=
T
⋅
D
F
T
[
x
(
n
)
]
N
k
=
0
,
1
,
⋯
,
N
−
1
(6)
\tilde X'_a(kF)=TX(k)=T\cdot DFT[x(n)]_N \qquad k=0,1,\cdots,N-1 \tag{6}
X~a′(kF)=TX(k)=T⋅DFT[x(n)]Nk=0,1,⋯,N−1(6)
式中,F表示对模拟信号频谱的采样间隔,所以称之为频率分辨率,
T
p
=
N
T
T_p=NT
Tp=NT为截断时间长度。
F
=
1
T
p
=
1
N
T
=
F
s
N
F=\frac{1}{T_p}=\frac{1}{NT}=\frac{F_s}{N}
F=Tp1=NT1=NFs
1.5 对推导结果的分析
(6)式表明可以通过对连续时间采样并进行DFT再乘以T,得到模拟信号频谱的周期延拓函数在第一个周期
[
0
,
F
s
]
[0,F_s]
[0,Fs]上的N点等间隔采样
X
~
a
′
(
k
F
)
\tilde X'_a(kF)
X~a′(kF)。
如图所示,对满足假设的持续时间有限的带限信号,在满足时域采样定理时,
X
~
a
′
(
k
F
)
\tilde X'_a(kF)
X~a′(kF)包含了模拟信号频谱的全部信息(
k
=
0
,
1
,
2
,
⋯
,
N
/
2
k=0,1,2,\cdots,N/2
k=0,1,2,⋯,N/2,表示正频率频谱采样;
k
=
N
/
2
+
1
,
N
/
2
+
2
,
⋯
,
N
−
1
k=N/2+1,N/2+2,\cdots,N-1
k=N/2+1,N/2+2,⋯,N−1,表示负频率频谱采样)。所以上述分析方法不丢失信息,即可由
X
(
k
)
X(k)
X(k)恢复
X
a
(
j
Ω
)
X_a(jΩ)
Xa(jΩ)或
x
a
(
t
)
x_a(t)
xa(t)。
但直接由分析结果
X
(
k
)
X(k)
X(k)看不到
X
a
(
j
Ω
)
X_a(jΩ)
Xa(jΩ)的全部频谱特性,而只能看到N个离散采样点的谱线,这就是所谓的栅栏效应。对实信号,其频谱函数具有共轭对称性,所以分析正频率频谱就足够了。不存在频谱混叠失真时,正频率
[
0
,
F
s
/
2
]
[0,F_s/2]
[0,Fs/2]的频谱采样为
X
~
a
′
(
k
F
)
=
T
X
(
k
)
=
T
⋅
D
F
T
[
x
(
n
)
]
N
k
=
0
,
1
,
⋯
,
N
/
2
\tilde X'_a(kF)=TX(k)=T\cdot DFT[x(n)]_N \qquad k=0,1,\cdots,N/2
X~a′(kF)=TX(k)=T⋅DFT[x(n)]Nk=0,1,⋯,N/2
如果 x a ( t ) x_a(t) xa(t)持续时间无限长,上述分析中要进行截断处理,所以会产生所谓的截断效应,从而使谱分析产生误差。
1.6 举例说明截断效应
理想低通滤波器的单位冲激响应 h a ( t ) h_a(t) ha(t)及其频响函数 H a ( f ) H_a(f) Ha(f)如下图:
图中
h
a
(
t
)
=
s
i
n
[
π
(
t
−
α
)
]
π
(
t
−
α
)
α
=
T
p
2
h_a(t)=\frac{sin[\pi(t-\alpha)]}{\pi(t-\alpha)} \qquad \alpha=\frac{T_p}{2}
ha(t)=π(t−α)sin[π(t−α)]α=2Tp
现在用DFT来分析
h
a
(
t
)
h_a(t)
ha(t)的频率响应特性。由于
h
a
(
t
)
h_a(t)
ha(t)的持续时间为无限长,因此要截取一段
T
p
T_p
Tp,假设
T
p
=
8
T_p=8
Tp=8s,采样间隔
T
=
0.25
T=0.25
T=0.25s(即采样频率
F
s
F_s
Fs=4Hz),采样点数
N
=
T
p
/
T
=
32
N=T_p/T=32
N=Tp/T=32,频率采样间隔即频率分辨率
F
=
1
/
T
p
=
0.125
F=1/T_p=0.125
F=1/Tp=0.125Hz;由于
h
a
(
t
)
h_a(t)
ha(t)为实信号,因此仅取正频率
[
0
,
F
s
/
2
]
[0,F_s/2]
[0,Fs/2]频谱采样:
H
(
k
F
)
=
T
⋅
D
F
T
[
h
(
n
)
]
0
≤
k
≤
16
H(kF)=T\cdot DFT[h(n)] \qquad 0\leq k \leq 16
H(kF)=T⋅DFT[h(n)]0≤k≤16
其中,
h
(
n
)
=
h
a
(
n
T
)
R
32
(
n
)
h(n)=h_a(nT)R_{32}(n)
h(n)=ha(nT)R32(n)。
H ( k F ) H(kF) H(kF)如上图(c)黑点所示。由图可见,低频部分近似理想低通频响特性,而高频误差较大,且整个频响都有波动。这些误差就是由于对 h a ( t ) h_a(t) ha(t)截断所产生的,所以通常称之为截断效应。为减少这种截断误差,可适当加长 T p T_p Tp,增加采样点数N或用窗函数处理后再进行DFT。
1.7 谱分析范围及频率采样间隔
在对连续信号进行谱分析时,主要关心两个问题,即谱分析范围和频率分辨率。
谱分析范围为 [ 0 , F s / 2 ] [0,F_s/2] [0,Fs/2],直接受采样频率 F s F_s Fs的限制。为了不产生频谱混叠失真,通常要求信号的最高频率 f c ≤ F s / 2 f_c\leq F_s/2 fc≤Fs/2。频率分辨率用频率采样间隔 F F F描述, F F F表示谱分析中能够分辨的两个频率分量的最小间隔。显然, F F F越小,谱分析就越接近 X a ( j f ) X_a(jf) Xa(jf),所以 F F F越小,频率分辨率越高。
1.8 DFT对连续信号谱分析的参数选择原则
在已知信号的最高频率
f
c
f_c
fc(即谱分析范围)时,为了避免频率混叠现象,要求采样速率
F
s
F_s
Fs满足下式:
F
s
>
2
f
c
F_s>2f_c
Fs>2fc
谱分辨率
F
=
F
s
/
N
F=F_s/N
F=Fs/N,如果保持采样点数
N
N
N不变,要提高频谱分辨率(减小
F
F
F)就必须降低采样频率,采样频率的降低会引起谱分析范围变窄和频谱混叠失真。如维持
F
s
F_s
Fs不变,为提高频谱分辨率可以增加采样点数
N
N
N,因为
N
T
=
T
p
NT=T_p
NT=Tp,
T
=
1
/
F
s
T=1/F_s
T=1/Fs,只有增加对信号的观察时间
T
p
T_p
Tp,才能增加
N
N
N。
T
p
T_p
Tp和
N
N
N的选择:
N
>
2
f
c
F
N>\frac{2f_c}{F}
N>F2fc
T
p
≥
1
F
T_p\ge \frac{1}{F}
Tp≥F1
2、用DFT对离散信号(序列)进行谱分析
单位圆上的Z变换就是序列的傅里叶变换,即
X
(
e
j
w
)
=
X
(
z
)
∣
z
=
e
j
w
X(e^{jw})=X(z)\Big|_{z=e^{jw}}
X(ejw)=X(z)∣∣∣z=ejw
X
(
e
j
w
)
X(e^{jw})
X(ejw)是
w
w
w的连续周期函数。如果对序列
x
(
n
)
x(n)
x(n)进行
N
N
N点DFT得到
X
(
k
)
X(k)
X(k),则
X
(
k
)
X(k)
X(k)是在区间
[
0
,
2
π
]
[0,2\pi]
[0,2π]上对
X
(
e
j
w
)
X(e^{jw})
X(ejw)的
N
N
N点等间隔采样,频谱分辨率就是采样间隔
2
π
/
N
2\pi/N
2π/N。因此序列的傅里叶变换可利用DFT(即FFT)来计算。
对周期为N的周期序列
x
~
(
n
)
\tilde x(n)
x~(n),根据周期序列的傅里叶变换表达式知,其频谱函数为
X
(
e
j
ω
)
=
F
T
[
x
~
(
n
)
]
=
2
π
N
∑
k
=
−
∞
∞
X
~
(
k
)
δ
(
ω
−
2
π
N
k
)
X(e^{j\omega})=FT[\tilde x(n)]=\frac{2\pi}{N}\sum_{k=-∞}^{∞}\tilde X(k)\delta (\omega-\frac{2\pi}{N}k)
X(ejω)=FT[x~(n)]=N2πk=−∞∑∞X~(k)δ(ω−N2πk)
其中
X
~
(
k
)
=
D
F
S
[
x
~
(
n
)
]
=
∑
n
=
0
N
−
1
x
~
(
n
)
e
−
j
2
π
N
k
n
\tilde X(k)=DFS[\tilde x(n)]=\sum_{n=0}^{N-1}\tilde x(n)e^{-j\frac{2\pi}{N}kn}
X~(k)=DFS[x~(n)]=n=0∑N−1x~(n)e−jN2πkn
由于
X
~
(
k
)
\tilde X(k)
X~(k)以N为周期,因而
X
(
e
j
ω
)
X(e^{j\omega})
X(ejω)也是以
2
π
2\pi
2π为周期的离散谱,每个周期有N条谱线,第k条谱线位于
ω
=
(
2
π
/
N
)
k
\omega=(2\pi/N)k
ω=(2π/N)k处,代表
x
~
(
n
)
\tilde x(n)
x~(n)的k次谐波分量。而且,谱线的相对大小与
X
~
(
k
)
\tilde X(k)
X~(k)成正比。
由此可见,周期序列的频谱结构可用其离散傅里叶级数系数
X
~
(
k
)
\tilde X(k)
X~(k)表示。由DFT的隐含周期性知道,截取
x
~
(
n
)
\tilde x(n)
x~(n)的主值序列
x
(
n
)
=
x
~
(
n
)
R
N
(
n
)
x(n)=\tilde x(n)R_N(n)
x(n)=x~(n)RN(n),并进行N点DFT,得到:
X
(
k
)
=
D
F
T
[
x
(
n
)
]
N
=
D
F
T
[
x
~
(
n
)
R
N
(
n
)
]
=
X
~
(
k
)
R
N
(
k
)
X(k)=DFT[x(n)]_N=DFT[ \tilde x(n)R_N(n) ]=\tilde X(k)R_N(k)
X(k)=DFT[x(n)]N=DFT[x~(n)RN(n)]=X~(k)RN(k)
所以可用
X
(
k
)
X(k)
X(k)表示
x
~
(
n
)
\tilde x(n)
x~(n)的频谱结构。
3、用DFT进行谱分析的误差问题
DFT(实际中用FFT计算)可用来对连续信号和数字信号进行谱分析。在实际分析过程中,要对连续信号采样和截断,有些非时限数据序列也要截断,由此可能引起分析误差。
3.1 混叠现象
- 对连续信号进行谱分析时,首先要对其采样,变成时域离散信号后才能用DFT(FFT)进行谱分析。
- 采样速率 F s F_s Fs必须满足采样定理,否则会在 w = π w=\pi w=π(对应模拟信号 f = F s / 2 f=F_s/2 f=Fs/2)附近产生频谱混叠现象。这时用DFT分析的结果必然在 f = F s / 2 f=F_s/2 f=Fs/2附近产生较大误差。因此,理论上必须满足 F s ≥ 2 f c F_s≥2f_c Fs≥2fc( f c f_c fc为连续信号的最高频率)。
- 对 F s F_s Fs确定的情况,一般在采样前进行预滤波,滤除高于折叠频率 F s / 2 F_s/2 Fs/2的频率成分,以免发生频率混叠现象。
3.2 栅栏效应
N点DFT是在频率区间 [ 0 , 2 π ] [0,2\pi] [0,2π]上对时域离散信号的频谱进行N点等间隔采样,而采样点之间的频谱是看不到的。这就好像从N个栅栏缝隙中观看信号的频谱情况,仅得到N个缝隙中看到的频谱函数值。因此称这种现象为栅栏效应。
由于栅栏效应,有可能漏掉(挡住)大的频谱分量。为了把原来被“栅栏”挡住的频谱分量检测出来,就必须提高频率分辨率。
(1)对有限长序列,可以在原序列尾部补零。
(2)对无限长序列,可以增大截取长度及DFT变换区间长度,从而使频域采样间隔变小,增加频域采样点数和采样点位置,使原来漏掉的某些频谱分量被检测出来。
(3)对连续信号的谱分析,只要采样率
F
s
F_s
Fs足够高,且采样点数满足频率分辨率要求(前述参数选择原则),就可以认为DFT后得到的离散谱的包络近似代表原信号的频谱。
3.3 截断效应
实际中遇到的序列 x ( n ) x(n) x(n)可能是无限长的,用DFT对其进行谱分析时,必须将其截短,形成有限长序列 y ( n ) = x ( n ) w ( n ) y(n)=x(n)w(n) y(n)=x(n)w(n), w ( n ) w(n) w(n)称为窗函数,长度为N。 w ( n ) = R N ( n ) w(n)=R_N(n) w(n)=RN(n),称为矩形窗函数。
根据傅里叶变换的频域卷积原理,有
Y
(
e
j
w
)
=
F
T
[
y
(
n
)
]
=
1
2
π
X
(
e
j
w
)
∗
W
(
e
j
w
)
=
1
2
π
∫
−
π
π
X
(
e
j
θ
)
W
(
e
j
(
w
−
θ
)
)
d
θ
Y(e^{jw})=FT[y(n)]=\frac{1}{2\pi}X(e^{jw})*W(e^{jw})=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\theta})W(e^{j(w-\theta)})\,d\theta
Y(ejw)=FT[y(n)]=2π1X(ejw)∗W(ejw)=2π1∫−ππX(ejθ)W(ej(w−θ))dθ
其中
X
(
e
j
w
)
=
F
T
[
x
(
n
)
]
W
(
e
j
w
)
=
F
T
[
w
(
n
)
]
X(e^{jw})=FT[x(n)] \qquad W(e^{jw})=FT[w(n)]
X(ejw)=FT[x(n)]W(ejw)=FT[w(n)]
对矩形窗函数
w
(
n
)
=
R
N
(
n
)
w(n)=R_N(n)
w(n)=RN(n),有
W
(
e
j
w
)
=
F
T
[
w
(
n
)
]
=
e
−
j
w
N
−
1
2
s
i
n
(
w
N
/
2
)
s
i
n
(
w
/
2
)
=
W
g
(
w
)
e
j
φ
(
w
)
W(e^{jw})=FT[w(n)]=e^{-jw\frac{N-1}{2}}\frac{sin(wN/2)}{sin(w/2)}=W_g(w)e^{j\varphi(w)}
W(ejw)=FT[w(n)]=e−jw2N−1sin(w/2)sin(wN/2)=Wg(w)ejφ(w)
W
g
(
w
)
=
s
i
n
(
w
N
/
2
)
s
i
n
(
w
/
2
)
W_g(w)=\frac{sin(wN/2)}{sin(w/2)}
Wg(w)=sin(w/2)sin(wN/2)
例如,
x
(
n
)
=
c
o
s
(
π
4
n
)
x(n)=cos(\frac{\pi}{4}n)
x(n)=cos(4πn),其频谱为
X
(
e
j
w
)
=
π
∑
l
=
−
∞
∞
[
δ
(
w
−
π
4
−
2
π
l
)
+
δ
(
w
+
π
4
−
2
π
l
)
]
X(e^{jw})=\pi \sum_{l=-∞}^{∞}\Big[\delta(w-\frac{\pi}{4}-2\pi l) +\delta(w+\frac{\pi}{4}-2\pi l) \Big]
X(ejw)=πl=−∞∑∞[δ(w−4π−2πl)+δ(w+4π−2πl)]
将
x
(
n
)
x(n)
x(n)截断后,
y
(
n
)
=
x
(
n
)
R
N
(
n
)
y(n)=x(n)R_N(n)
y(n)=x(n)RN(n)的频谱为在
δ
\delta
δ谱线处
W
g
(
w
)
W_g(w)
Wg(w)图像的展宽。
由此可见,截断后序列的频谱 Y ( e j w ) Y(e^{jw}) Y(ejw)与原序列 X ( e j w ) X(e^{jw}) X(ejw)必然有差别,这种差别对谱分析的影响主要表现在如下两个方面:
3.3.1 泄露
原来序列 x ( n ) x(n) x(n)的频谱是离散谱线,经截断后,使原来的离散谱线向附近展宽,通常称这种展宽为泄露。显然,泄露使频谱变模糊,使谱分辨率降低。
从上图中可以看出,频谱泄露程度与窗函数幅度谱的主瓣宽度直接相关,而在所有的窗函数中,矩形窗的主瓣是最窄的,但其旁瓣的幅度也最大。所以,在窗函数长度N相同时,用矩形窗截取,产生的泄露最小。
3.3.2 谱间干扰
在主谱线两边形成很多旁瓣,引起不同频率分量间的干扰(简称谱间干扰),特别是强信号谱的旁瓣可能湮没弱信号的主谱线,或者把强信号谱的旁瓣误以为是另一频率的信号的谱线,从而造成假信号,这样就会使谱分析产生较大偏差。由于矩形窗的旁瓣最大,所以,用矩形窗截取时,产生的谱间干扰最大。
泄露与谱间干扰
由于泄露和谱间干扰是由信号截断引起的,因此称之为截断效应。矩形窗 ∣ ω ∣ < 2 π N |\omega|<\frac{2\pi}{N} ∣ω∣<N2π的部分为主瓣,增加N可使 W g ( ω ) W_g(\omega) Wg(ω)的主瓣变窄,减小泄露,提高频率分辨率,但旁瓣的相对幅度并不减小。为了减小谱间干扰,应用其它形状的窗函数 w ( n ) w(n) w(n)代替矩形窗。但在N一定时,旁瓣幅度越小的窗函数,其主瓣就越宽。所以,在DFT变换区间(即截取长度)N一定时,只能以降低谱分析分辨率为代价,换取谱间干扰的减小。
栅栏效应与频率分辨率
栅栏效应与频率分辨率是不同的两个概念。如果截取长度为N的一段数据序列,则可以在其后面补N个零,再进行2N点DFT,使栅栏宽度减半,从而减轻了栅栏效应。
但是,这种截短后补零的方法不能提高频率分辨率。因为截短已经使频谱变模糊,补零后仅使采样间隔变小,但得到的频谱采样仍是已经变模糊的频谱,所以频率分辨率没有提高。因此,要提高频率分辨率,就必须对原始信号截取的长度加长(对模拟信号,就是增加采样时间 T p T_p Tp的长度)。