-
前几天上数字信号处理(本以为第二次上这个课只是简单地重复过去学习过的内容,但是这次有了
很多新
的发现),书上说对时域信号补零之后再作
DFT
并不能提高频谱的频率分辨率,提高采样频率也不能提高
DFT
谱的频率分辨率。这个很新鲜,以前上课时没有考虑过这个问题,以前的课本好像也没有开辟专
门的
章节论述这个问题。
提高采样频率不能提高频率分辨率的原因其实很简单,因为提高了采样频率,虽然在相同
的观察时长
那的点数增多了,但与此同时采样频率也变大了,点数增加几倍采样频率增加
几倍,所以不改变观察时长
而仅仅提高采样频率并不能提高
DF
T
谱的频率分辨率。
<
/p>
但是时域补零呢?采样频率没有变化,而点数增加无疑会减小
DF
T
谱的相邻谱线间隔,相邻谱线间
隔的缩小为什么不能提高频率
的分辨率呢?书上是这样写的:
“
错把
?
计算分辨率
?
当成了
?
物理分辨率
?……
补
零没有对原信号增加任何新的信息,因此不可能提高分辨率。但补零
…
…
补零还可以对原
X(k)
做插值。<
/p>
”
(《数字信号处理
——
理论、算法与实现(第二版)》清华大学出版社,胡广书)
为了更好地理解这个问题,我又一次借用
MATLAB
的强大力量,写了一个简单的程序如下:
< br>
%
点数
n=0:127;
%
频率
f=0.1;
%
信号,正弦叠加矩形
y1=sin(2*pi*f*n);
y1(1:16)=y1(1:16)+1;
%
绘制
y1
的
fft
谱幅度
%
谱线较多,直接画的包络
figure;
plot(abs(fft(y1)));
%
对信号进行截短
y2s=y1(1:32);
<
/p>
%
绘制
y1
截断
后没有补零的
fft
谱幅度
figure;
fy2s=abs(fft(y2s));
stem(fy2s);
%
然后补零使
y1
和
y
2
一样长
y2=[y2s zeros(1,128-32)];
%
打开一个绘图窗口
figure;
%
绘制
y1
的
fft
谱幅度
%
谱线较多,直接画包络
plot(abs(fft(y1)));
%
在同一个
figure
中继续绘图
hold on;
%
绘制
y2
的
fft
谱幅度
(
红色
)
%
谱线较多,直接画包络
plot(abs(fft(y2)),
'r'
);
%
绘制
y2s
的
fft
幅度谱
stem(
1:4:128,fy2s,
'k'
);
hold off;
程序生成的图像如下:
图
1
p>
(原信号的谱,因为点数较多,绘制的是包络)
图
2
p>
(截短信号的谱)
图
3
(比较
原信号【蓝】,截短信号【黑】,补零信号【红】三者谱的关系,补零信号的谱由于点数较多,
< br>绘制的是包络)
程序的大意是:首先生成了一个长度为
128
< br>点的信号,绘制了它的
DFT
谱,然后将该信号截短,求
其
DFT
谱,然后对截短的信号补零,
使其长度为
128
点再求
DFT
谱,并将原信号的谱与补零信号的谱进
行比较,结果一目了然。
从图三中我们可以发现,红色
的补零信号的谱,仅仅是对黑色的截短信号的谱的插值,也就是说补零
信号的谱,是通过
截短信号的谱进行了推测(插值算法)得来的,它并不能反映原信号的谱(因为原信号
在
截短的过程中部分信息丢失了,而补零并没有将这些丢失的信息找回来),所以虽然补零信号的谱线间
隔变小了,
但是除了从截短信号的谱中取出来的
3
2
根谱线以外,
其余的
96
根谱线都是无效的。
去掉这些
无效的谱线,采样频
率不变,有效的谱线数不变,所以其物理频率分辨率自然没有改变。
在时域补零起到了对频谱的插值作用,考虑到傅立
叶变换的对称性,考虑到变换与反变换的表达式除
了系数和符号的差别外并没有根本的不
同,推断出这样一个结果:在频域补零也会造成时域的插值(当然
在频域进行操作一定要
注意谱的对称性,否则反变换回来得到的将是一个复信号,至于为什么会这样参看
《
p>
xialulee
来说明什么是负频率》
),于是我又写了如下的程序:
f=0.1;
n=0:31;
%
产生正弦信号
y1=sin(2*pi*f*n);
%y1
的
f
ft
谱
f1=fft(y1);
%
给
y1
的谱补零
f2=[f1(1:
end
/2+1)
zeros(1,31) f1(
end
/2+1:
end
)];
stem(abs(f2));
%
对补零后的谱作反变换
y2=ifft(f2);
%
打开绘图窗口
figure;
%
绘制多幅图形
hold on;
%
绘制
y2
stem(y2);
%
绘制原信号(除以
2
以比较插值效果)
stem(1:2:64,y1/2,
'r'
);
hold off;
程序的输出图形如下:
图
4
仿佛证明了这个推测。这是不是就是所谓的傅立叶插值
关于
DFT
计算正弦没有频谱泄漏的问
题
本来是昨天写的文章,但是本着
一天发表不超过一篇的原则,就拖到今天了。
星期二上数字信号处理,
老师讲解了
“
关于正弦信号抽样的讨论
”
一节。
其中有一小节是
“
对正弦信号截
短的原则
”
,提出:如果抽样频率大于正弦频率的两倍,且是正弦信号频率的整数倍;如果抽样点数包含整
< br>数个正弦周期,则对应的
DFT
谱没有泄漏。
符合上述要求的正弦信号,仍然是
截短的正弦信号,仍然相当于乘了一个矩形窗,但是
DFT
谱没
有
出现泄漏的现象
——
书中配合一幅插
图讲解
——
是因为虽然
DTFT
谱出现了泄漏,
但
DFT
是对
DTFT
的抽
样,而此时,除
了两个点采到了波峰以外,其余采样点都采到了
DTFT
谱曲线
的过零点上。
有插图就
能很形象地说明这个问题,为了绘制与书上类似的插图,我编写了如下的
Matlab<
/p>
程序
(
Matlab
版本
7.6
):
文件名
dtft.m
用于
DTFT
运算
function
X=dtft(x,range,w)
X=x*exp(-j*(range(1):range(2)).'*w);
p>
文件名
testdtft.m
产生信号并运
算并绘图
%
信号
x(n)
的
n
n=0:15;
%dtft
的点数,虽然
dtft
是连续谱,但连续量无法用数字计算机计算
nw=513;
%
产生信号,一个周期采样
16
个点
x=cos(2*pi/16*n);
%
产生用于
dtft
的频率向量
w=linspace(0,2*pi,nw);
%
计算
x
的
dtft
谱
Xdtft=dtft(x,[n(1),n(
end
< br>)],w);
%
用
fft
计算
x
的
dft
谱
Xdft=fft(x);
%
打开画图窗口
figure;
%
绘制多幅图形
hold on;
%
绘制
dtft
谱
plot(abs(Xdtft));
%
绘制
df
t
谱
st
em(1:((nw-1)/numel(n)):nw-1,abs(Xdft),
'
r'
);
运行
testdtft.m
,得到了如下的图形:
从上图中可以看出,
除了两个点以外,
其余都采到了过零点上,
于是
DFT
谱
“
好似
”
没有
发生频谱泄漏。
(
DFT
谱上非零的两
个点没有采到
DTFT
谱的波峰上,不知是我的程序出了问题还
是什么
“
有限字长效应
”
抑或是什么我不太了解的原因,但是此图仅仅用来说明一下情况,所以目的已经达到,没有深究。我想
最
可能的原因是两个
sinc
的旁瓣影
响了对方的主瓣,从而使得主瓣的最大值点发生了偏移)
从采样采到过零点这一角度可以说明:为什么对于有限长序列
D
FT
仍然有可能做到没有频谱泄漏。
不过我觉得这个问题还可以
从另一个方面进行说明:
DFT
是由
DFS
导出来的,
DFT
是将时域有限长的序列进行拓展,变成一个无限长的周期信号,
之后
对这个无限长的周期信号计算
DFS
,
DFS
同样是一个无限长的周期序列,再截取
DFS
的主周期,得到的
就是所谓
< br>DFT
谱了。
<
/p>
因此,对于满足了那两个条件的正弦信号(如果抽样频率大于正弦频率的两倍,且是正弦信
号的整数
倍;如果抽样点数包含整数个正弦周期),对其进行拓展,则得到了一个完整的
无限长的正弦序列。因此,
频谱泄漏在此处被消除,于是后面的
DFT
谱中看不到泄漏的现象了。
所以,总结一下,我觉得为什么
DFT
谱对于上述的信号没有出现频谱泄漏可以从两个不同的角度进
行说明:
< br>
角度一:
DFT
是对
DTFT
的采样,
DTFT
的确存在着由于信号截短而产生的泄漏,但是
DFT
除了两
个点外,其余采样点都采到了
DTFT
谱曲线的过零点了,所以在
DFT
中看不到频谱泄漏;
< br>角度二:
DFT
由
DFS
导出,隐含着周期性。对时域有限长序列进行拓展,相当于取消了矩形窗,导
< br>致频谱泄漏现象的消失。
最后说一点题外话,有一种免费的类似
Matlab
的
软件叫做
Scilab
。现在推出了
5
.0
的版本,我试验
了一下它的
Mat
lab
代码自动转换成
Scilab
代
码的功能,发现有些地方还是得手动修改才能运行,代码如
下:
文件名
function [X] = dtft(x,range,w)
// Ouput variables
initialisation (not found in input variables)
X=[];
// Display
mode
mode(0);
//
Display warning for floating point exception
ieee(1);
X =
x*exp(-(%i*(mtlb_imp(mtlb_double(mtlb_e(ra
nge,1)),mtlb_double(mtlb_e(range,2))).'))*
< br>w);
endfunction
文件名
// Display mode
mode(0);
// Display warning for
floating point exception
ieee(1);
//
信号
x
(n)
的
n
n = 0:15; <
/p>
//dtft
的点数,虽然
dtft
p>
是连续谱,但连续量无法用数字计算机计算
nw = 513;
//
产生信号,
一个周期采样
16
个点
x = cos(((2*%pi)/16)*n);
//<
/p>
产生用于
dtft
的频率向量
w = linspace(0,2*%pi,nw);
//
计算
x
的
dtft
谱
Xdtft
= dtft(x,[n(1),n($$)],w);
//
用
fft
计算
x
的
dft
谱
Xdft = fft(x,-1);
//
打开画图窗口
// !! L.14: Matlab function figure not
yet converted.
// !! L.14: Matlab
function figure not yet converted, original
calling sequence used.
figure;
//
绘制多幅图形
< br>set(gca(),
//
绘制
dtft
谱
plot(abs(Xdtft));
//
< br>绘制
dft
谱
// !! L.20: Matlab function numel not
yet converted, original calling sequence used.
// !! L.20: Matlab function stem not
yet converted, original calling sequence used.
plot2d3(mtlb_imp(1,(nw-1)/length(n),nw-1),
abs(Xdft));
上面的代码是自动转换之后我又手动修改的,有两处,一个是
Matlab
的
numel
函数,在
Scilab
中
应换成
length<
/p>
吧。另一个是
stem
,改成了
plot2d3
,去掉了第三个用于设置线条属性的字符串参数。
由于没有编写代码上色程序,所以只好黑白贴上来了。
关于
FFT
的相位谱
昨天收到一条消息,发送消息的人表示出对
FFT
相位谱的忧虑,他给出了一段
Matla
b
程序,显示了
FFT
的相位谱并没有
真实反映正弦信号的初相。
今天我研究了一下这个问题,得到了一些比较有启发的结论。本来想直接以消息的形式对问题进行答
复,但是考虑到消息中不太好控制图片的显示方式,而且这个问题如果能和更多的人进行讨论,想
必会得
到更多的启发,所以直接发到了这里。
先看一下我收到的程序,作为研究对象的信号是这样产生的:
T=128;
N=128;
dt=T/N;
t=dt*(1:N);
x=2*cos(2*t-pi/4);
...
(我觉得这个信号存在一点问题,因为
t
是从
1
开始的,所以它
的初相应该和
-pi/4
有点差别吧。)
为什么进行
FFT<
/p>
,用
angle
得到相位-频率特性却不
能反映这个信号的初始相位?
< br>胡广书的
《数字信号处理-理论、
算法与实现
(第二版)
》
第三章第八节
《关于正弦信号抽样的讨论》
,
得出了关于正弦信号抽样的六
个结论,最后总结了一个总的原则:抽样频率应为信号频率的整数倍,抽样
点数应包含整
周期。
书中的结论五与
采样频率和抽样点数有很大的关联。结论五主要说只有满足了上面的那个总的原则,
频谱
泄漏才不会发生。
我想不光是幅度谱的频谱泄漏现象,
抽样频率
和抽样点数同样会对相位谱产生影响。
考虑一个无限长的正弦信号(相当于初相为-
90°
),如果我们截取它的整数个周期,然后对截短的
信号进行周期延拓,则得到的延拓
的信号与原无限长正弦没有区别。
现在我们再次对这个无限长的正弦进行截短,长度为
1.5
< br>个周期,然后对截短信号进行周期延拓,看