matlab编写dft函数_DFT之泄露

更多经验分享，关注 加油射频工程师 

(1)漫画

(2)引言

自从同事给推荐Understanding Digital Signal Processing这本书来，断断续续地一直再看。最近，终于看的稍微明白了一点。

这章共看了三遍：

(1) 看第一遍的时候，就是不管懂不懂，先往下看，所以看完后，感觉就像没看一样；

(2) 看第二遍的时候，就边看边翻译成中文，有不懂的，也放过；但是因为有形成文档笔记，所以本身心里还是比较有成就感的。这样就形成了正循环，看的时候，心情比较愉悦，所以理解的东西也会比较容易。

(3) 看第三遍的时候，就尝试着用自己的语言组织，用matlab去仿真验证，不懂的再去查其他资料帮助理解。

不得不说，写公众号确实是一个促进学习的好方式，以写促学。因为想在公众号里发表出来，就会去很认真地去写文档，然后就会把那些似懂非懂的问题去搞懂。

根据书本上的DFT上的内容，共整理了以下几篇文档：

DFT学习笔记，是关于DFT的变换公式以及其的一些特性

DFT之泄露，关于泄露的现象以及产生的原因；

DFT之窗函数，列举了几种窗函数，并进行了相关比较，以及窗函数的一些作用

DFT之扇贝损耗，以矩形函数为例，叙述了扇贝损耗的现象

DFT之零填充，指出进行零填充的原因，以及注意事项；

DFT之处理增益；

第一篇可以从往期文章看到，今天这篇主要是关于泄露。

(3)DFT之泄露

泄露的现象

DFT即是，对具有N个数的输入序列，进行一个N点的DFT变换，产生一个具有N个点的频域序列。而且，这些频域序列的间隔为：

前面的例子DFT学习笔记，之所以能够准确计算出频率，是因为输入信号的频率是fs/N的整数倍。

那如果不是整数倍，会怎么样呢？

还是上面的例子，把第二项的频率从2000改为1500，然后再运行一下程序，可以得到下述结果。并且将频率为2000Hz的计算结果列于左侧。

可以看到:

(1) 在m=1.5处未产生频谱，

(2) m=2,3,4,5处应该是0的，现在也有值了

(3) m=1,7处的值与频率为2000Hz时的结果相比，变小了。

从上面例子可以看出，当输入序列的频率不是fs/N的整数倍时，比如1.5fs/N时，其能量会被泄露到频域上的N个序列点上。

因此，在用DFT处理真实世界中的有限长的时间序列时，泄露是不可避免的。

那为什么会产生泄露呢？

这是因为，在使用DFT时，会不可避免的使用窗口函数。

以下列函数为例：

为一连续时间函数。

为了计算其DFT，我们对其以fs=8000对其采样，得到xin[n]，即

但我们进行DFT运算时，只是使用了其前8个点，记为x[n]。所以可以认为x[n]=xin[n]w[n]。

绘制上图的matlab程序：

clc;

clear all;

ts=1/8000;

N=64;

n=0:1:(N-1);

m=0:1:7;

Xn=sin(2*pi*1500*n*ts);

Xm=sin(2*pi*1500*m*ts);

subplot(3,1,1)

%画原始的连续函数

plot(n,Xn)

title('原始连续函数')

ylabel('xin(t)')

%绘制采样后的函数

subplot(3,1,2)

stem(n,Xn,'k')

hold on;

stem(m,Xm,'r--')

title('黑色为采样后的点，红色为DFT变换所用的点')

ylabel('xin(n)')

%绘制窗函数

subplot(3,1,3)

stem(m,rectwin(8),'b')

hold on;

k=8:1:70

stem(k,zeros(63),'b')

title('矩形窗函数')

ylabel('w(n)')

也就是说，当我们想通过DFT换算计算xin(t)频谱的时候，我们其实是做了这些步骤：

(1) 对连续信号xin(t)进行采样，得到离散信号xin[n]

(2) 截取离散信号xin[n]中的一段，即与窗函数w[n]相乘，得到x[n]。

所以，我们其实计算的是x[n]=x_in[n]•w[n]的频谱，即对xin[n]进行矩形加窗后的函数的频谱。

问题是，加窗后的函数的频谱与原来的频谱有什么区别呢？

要说明DFT计算的频谱的问题，先了解一下DFT与DTFT的关系，以及DTFT变换。

DTFT与DFT之间的关系：

DTFT的变换公式为：

 

DFT的变换公式为：

 

所以，如果x[n]的长度为L，且L≤N，则其N点DFT为DTFT的采样，即：

 

而x[n]正好满足这一要求，所以可以对X(w)做以上操作，以得到其DFT。

 时域乘积等于频域卷积。

(1)先计算窗函数的DTFT：

窗函数

w[n]的DTFT的计算结果为：

 

当K=8时，其幅度函数如下图所示：

 

(2)计算x[n]的DTFT

l 对xin[n]做如下变换：

 

DTFT具有下列性质：

所以

用matlab画出x[n]的DTFT的幅度函数：

以上两副图的matlab实现程序：

clc;

clear all;

 K=8

 %计算矩形窗函数DTFT变换的幅度函数

w=-pi:0.01:pi;

Ww=exp(-j*w*(K-1)/2).*sin(w*K/2)./sin(w/2)

figure (1);

plot(w,abs(Ww))

set(gca,'xtick',[-pi:pi/4:pi])

set(gca,'XTickLabel',{'-pi';'-3pi/4';'-pi/2';'-pi/4';'0';'pi/4';'pi/2';'3pi/4';'pi'})

%计算x[n]的DTFT变换

n=0:1:(K-1);

w0=3*pi/8;

xn=sin(w0*n);

ww=0:0.01:2*pi;

w1=ww-w0;

w2=ww+w0;

Ww1=exp(-j*w1*(K-1)/2).*sin(w1*K/2)./sin(w1/2);

Ww2=exp(-j*w2*(K-1)/2).*sin(w2*K/2)./sin(w2/2);

xw=(Ww1-Ww2)/(2*j);

figure (2);

plot(ww,abs(xw));

set(gca,'xtick',[0:pi/4:2*pi])

set(gca,'XTickLabel',{'0';'1pi/4';'pi/2';'3pi/4';'pi';'5pi/4';'3pi/2';'7pi/4';'2pi'})

%计算x[n]的FFT变换，并与x[n]的DTFT变换比较

hold on;

Y=fft(xn,K);

stem(2*pi/K*n,abs(Y))

如果，则可以得到x[n]的DTFT与FFT之间的关系为：

(3)总结说明

  

从上面两幅图，我们可以看到：

(1) 加窗函数x[n]计算出来的DTFT是连续频谱，也就是说，其包含了所有的频率分量。(2) DFT是在DTFT上采样的结果，如果输入信号的频率正好等于采样点，则计算出来的频率是精确的。

(3) 如果输入信号的频率与采样频率点不等，则计算出来的频率则不等于输入信号频率，而且会有很多泄露。实际情况，经常是这样，输入信号频率≠采样频率点。

所以说，因为窗函数的存在，其旁瓣导致有泄露；所以可以用泄露小的窗函数来减小其影响。

参考文献：

[1]RICHARD G. LYONS Understanding digital signal processing

[2]J. Fessler, chapter 5, The Discrete Fourier Transform

部分图片及文字来源于网络，如有侵权，麻烦后台留言，立马删除，谢谢！

长按图片关注微信号