DFT
1. 频域采样和时域混叠¶
对于 \(x\left\lbrack n\right\rbrack =a^n u\left\lbrack n\right\rbrack ,\;\;\;0<a<1\)
\(X\left(\omega \right)\) 的频域如下图
a=0.8;
omega = pi*(0:0.01:2);
X = 1./(1-a*exp(-1i*omega));
plot(omega/pi,abs(X))
xlabel("freqs/(w/\pi)")
ylabel("maginitude")
采样N个点 \(X\left(\frac{2\pi k}{N}\right),k=0,1,\ldotp \ldotp \ldotp ,N-1\) ,重建得到的序列为
很显然，N越大，混叠的影响就越小
%比较重建和原来的时域
clf
N=0:1:100;
x = power(a,N);
hold on
plot(N,x)
n1=0:1:4;
xp1=power(a,n1)./(1-power(a,5));
plot(n1,xp1);
xlim([0 100]);
n2=0:1:50;
xp2= power(a,n2)./(1-power(a,50));
plot(n2,xp2);
xlim([0 100]);
hold off
2. DFT¶
下面注意定义，L是信号的采样点数，N是做DFT的点数，N>=L
%考虑下面的DTFT（和滑动平均滤波器是一样的）
L=10;
x = ones(1,L);
w = 0:0.01:10;
X = sin(w*L/2)./sin(w/2).*exp(-1i*w*(L-1)/2);
plot(w/pi,(abs(X)));
xlim([0,2]);
通过上面的例子，知道了DFT长度和补0会发生什么了吗？
首先，对于一个有限长L的序列x[n]来说，L点DFT是和时域唯一对应的，这没有错。
但是点数太少的时候，频域上信息太少了，什么都看不出来。
所以我们要增加DFT的长度N，也就是补零。补零的过程中相当于插值，增加了分辨率 \(\omega =\frac{2\pi }{N}\) （数值越小表示分辨率越大）我们得到了更多的频域的信息。换句话说，DFT的点数N就是在DTFT频域 \(\left\lbrack 0,2\pi \right\rbrack\) 上采样的点数。
注意：此处的分辨率（resolution）用于描述频域采样的间隔，实际上描述的是清晰度，和后面窗函数的分辨能力是不同的。
DFT进行频域分析¶
我们采样得到的信号的长度永远是有限的，也就是说不自觉的加上了矩形窗，窗会使功率泄露。
我们拿最简单的余弦信号进行举例（proakis的例子）
omega0 = 0.2*pi;
omega1 = 0.22*pi;
omega2 = 0.6*pi;
n = 0:1:24;
x = cos(omega0*n)+cos(omega1*n)+cos(omega2*n);
X1 = fft(x,2048);
n = 0:1:49;
x = cos(omega0*n)+cos(omega1*n)+cos(omega2*n);
X2 = fft(x,2048);
n = 0:1:1999;
x = cos(omega0*n)+cos(omega1*n)+cos(omega2*n);
X3 = fft(x,2048);
plot(abs(X1));
可以看到，截取 \(x\left\lbrack n\right\rbrack\) 的25个和50个点做fft都会吃掉0.22pi的部分，只有200个点的时候才能够分辨出来。这是为什么呢？
来看一下矩形窗的fourier transform（就是复制了开始的滑动平均滤波器）
L=10;
x = ones(1,L);
w = 0:0.01:10;
X = sin(w*L/2)./sin(w/2).*exp(-1i*w*(L-1)/2);
plot(w/pi,(abs(X)));
xlim([0,2]) %0--2\pi figure
主瓣的第一个零点是 \(\frac{2\pi }{L}\) ,这里L=10，所以零点是0.2π。
只有当 \(|\omega_1 -\omega_2 |>\frac{2\pi }{L}\) 的时候才能分辨出来。（画图看）（这严谨吗）
所以采样L个点， \(\frac{2\pi }{L}\) 就是分辨率。
或者说，主瓣的宽度表征了分辨率。宽度越窄，分辨率越高，能分辨差距更小的频率。
拿书本里的结论来说：分辨率受主瓣（main lobe）宽度的影响，泄露的程度取决于主瓣和旁瓣的相对幅度。
The resolution is influenced primarily by the width of the main lobe of W(ejω), whereas the degree of leakage depends on the relative amplitude of the main lobe to the side lobes of W(ejω).
但是上面的例子怎么看出泄露？正弦波原来只有一个频率，加窗后增加了很多频率的分量。旁瓣的幅度越大，自然说明泄露的越多咯。
矩形波的分辨率最好，但与此同时，它泄露的也是最多的。
（这里也强调一遍，和padding zero说的分辨率不是同一个分辨率）
汉宁窗¶
前面看到矩形窗泄露的频率到旁瓣上了，我们可以选择汉宁窗来降低旁瓣（代价是增加主瓣的宽度）
频率特性如下图
W=fftshift(fft(w,1000));
omega = 0:2*pi/1000:2*pi-2*pi/1000;
plot(omega/pi, abs(W));
title("freqency domain")
圆周卷积¶
若序列1的长度是L，序列2的长度是P，只有结果序列的长度是L+P-1时，它才会和线性卷积的结果一致。即做圆周卷积前要补零到合适的长度（L+P-1），否则就会混叠。
举例：
4点圆周卷积发生了混叠，7点的是和线性卷积答案一致。
看了上面也可以知道，padding-zero只会在数值上影响DFT的结果，但是不会影响卷积（滤波）的结果。
长数据滤波¶
数据很长时，常常要分块滤波。这里就衍生出两种方法：overlap-add和overlap-save。
overlap-save¶
8p和6p的结果是完全一致的，说明补零后fft可以用于线性滤波。
至于4点fft则出现了混叠：x[4]叠加到x[0]上，x[5]叠加到x[1]上。而后面是正确的。
即：若输入数据块的大小为N，做N点DFT，与FIR长度为M进行滤波，那么输出的前面M-1个点会出现混叠，后面都正确。
（更加广义的是做N1点DFT，则有N+M-1-N1个点出现混叠）
所以overlap-save的原理就是对N个点做DFT，舍弃前面M-1个点，保留后面的L个点。对每一段都这样。
overlap-add¶
输入padding几个0，输出末尾就有几个0
先padding 0，到L+M-1，得到的输出y也是L+M-1，每一段前面M-1个点和前一段后面M-1个点相加即可。
教材的图已经很形象。
ch10¶
第十章主要讲窗函数，和前面类似，但是更加详细
注意：我们考虑的都是窗函数对无限长信号处理后造成的影响。
example10.3¶
加窗对fourier分析的影响
采样率是10kHz，矩形窗长度为64；
矩形窗形状：
w = ones(1,64);
W = fftshift(fft(w,1024));
%fftshift以原点为中心了就
w = -pi:2*pi/1024:pi-2*pi/1024;
plot(w,abs(W));
xlim([-pi,pi])
fs = 10e3;
%x[n] = x(nT) = x(n/fs)
n=0:63;
omega1 = 2/6*pi*fs;
omega2 = 2/3*pi*fs;
y = cos(omega1*n/fs) + 0.75*cos(omega2*n/fs);
Y = fftshift(fft(y,1024));
%矩形窗只有63个点，补大量0可以增加更多细节。
w = -pi:2*pi/1024:pi-2*pi/1024;
plot(w,abs(Y))
xlim([-pi,pi])
omega1 = 2/14*pi*fs;
omega2 = 4/15*pi*fs;
y = cos(omega1*n/fs) + 0.75*cos(omega2*n/fs);
Y = fftshift(fft(y,1024));
w = -pi:2*pi/1024:pi-2*pi/1024;
plot(w,abs(Y))
xlim([-pi,pi])
omega1 = 2/12*pi*fs;
omega2 = 2/14*pi*fs;
y = cos(omega1*n/fs) + 0.75*cos(omega2*n/fs);
Y = fftshift(fft(y,1024));
w = -pi:2*pi/1024:pi-2*pi/1024;
plot(w,abs(Y))
xlim([-pi,pi])
上图中此时频率泄露最为明显，两个峰的峰值都受到对方的影响，远远小于32和24了，只是勉强可以分辨
omega1 = 2/15*pi*fs;
omega2 = 2/14*pi*fs;
y = cos(omega1*n/fs) + 0.75*cos(omega2*n/fs);
Y = fftshift(fft(y,1024));
w = -pi:2*pi/1024:pi-2*pi/1024;
plot(w,abs(Y))
xlim([-pi,pi])
长度L是64所以分辨率是2pi/64. 2pi/15-2pi/14<<2pi/64了所以当然不能分辨。
kaiser window¶
书本后面的例子会用到很多kaiser窗。
得益于MATLAB函数我们直接调用即可。
w = kaiser(20,3);
%200points and beta=3 kaiser window.
W = mag2db(abs(fft(w,4096)));
f = 0:2/4096:2-2/4096;
plot(f,W)
xlim([0,1])
xlabel("normed freq *\pi rad/sample")
title("freq domain")
绘制不同的Kaiser窗，我们可以比较一下窗函数的性质：
主瓣宽度和窗的长度成反比，相对旁瓣高度和窗的长度关系不大。
例10.6,7,8¶
用Kaiser窗对正弦信号DFT分析
omega1 = 2*pi/14;
omega2 = 4*pi/15;
%L=64, beta=5.48
w = kaiser(64,5.48);
n=0:63;
x = cos(omega1*n)+0.75*cos(omega2*n);
v = w' .* x;
X = fftshift(fft(x,1024));
V = fftshift(fft(v,1024));
f = -pi:2*pi/1024:pi-2*pi/1024;
plot(f,abs(X))
title("64 points square window")
n1=0:31;
w1 = kaiser(32,5.48);
x1 = cos(omega1*n1)+0.75*cos(omega2*n1);
v1 = w1' .* x1;
V1 = fftshift(fft(v1,1024));
f1 = -pi:2*pi/1024:pi-2*pi/1024;
plot(f1,abs(V1))
title("32 points kaiser window")
看上面3张图片，你就知道为什么说补零不能提高分辨率了（PS，这个分辨率是指the ability to resolve close frequencies, 分辨不同频率的能力）。补零增加的能力是频域上的清晰度，即在已经有的DTFT上采样更多的样本。而DTFT的形状（即分辨不同频率能力的本质）是由信号和窗函数的长度和形状决定的。
我们用窗函数是截取无限长的信号，才有这种真不真实的问题；有限长的信号本身如此，就没有这种问题了。
窗函数的长度形状和分辨能力的关系:
良好的形状可以减少相对旁瓣高度，减少频率泄露。
而窗函数的长度越长，主瓣宽度就越小，分辨能力越强。反过来，窗函数越短，数据DFT就越容易反映出窗函数的频谱，离我们想要的就越远。
下面看一下这道题：
proakis e7.4.1¶
%模拟信号e^(-t)的频谱
omega = -20:0.1:20;
X = 1./(1+1i*2*pi*omega);
plot(omega, abs(X));
title("original freq")
%窗是100 points时的频域
n=0:1:99;
x=power(0.95,n);
X1=fftshift(fft(x, 200));
plot(abs(X1))
title("100 points window")
%窗是20 points时的频域
n=0:1:19;
x=power(0.95,n);
X2=fftshift(fft(x,200));
plot(abs(X2));
title("20 points window")
%在只有20个点的时候，宽频谱矩形窗函数的效果就显现出来了。
%比如主峰和旁瓣，信号也泄露不少
% 信号点数越少 窗函数就越明显。
%这道题也可以说明DFT补零纯粹是为了增加清晰度，
% 和真实情况没有半毛钱关系
N1 = fftshift(fft(ones(1,1000),2048));
plot(abs(N1))
xlim([0,2000])
title("1000 points square window function")
