数字信号处理翻转课堂笔记8

The Flipped Classroom8 of DSP

对应教材：《数字信号处理（第五版）》西安电子科技大学出版社，丁玉美、高西全著

一、要点

1、用DFT来计算序列的线性卷积时需要满足的条件（重点）；
2、重叠相加法计算线性卷积；
3、用DFT分析连续（模拟）信号的频谱的原理（难点），频率分辨率、频谱分析范围、采样点数、记录时间等之间的关系（重点）；
4、用DFT分析序列的频谱需要注意的问题；
5、用DFT分析信号频谱所面临的频谱混叠、栅栏效应和截断效应。

二、问题与解答

1、用DFT来计算序列的线性卷积有什么优势？其原理和步骤是什么（画出其实现原理框图，并予以说明）？根据实现原理框图，分析用DFT计算线性卷积的计算量，并与根据定义直接计算两个序列线性卷积的计算量相比，哪一个计算量更大？

2、

3、针对教材例3.4.1，不使用fftfilt函数，直接基于重叠相加法的原理，自行编程计算其卷积。据此总结重叠相加法的原理和步骤，并从理论上分析为什么卷积结果稳态输出是一个正弦序列。

4、针对一个周期序列，分别截取其1个周期、多个周期、非整数周期，用MATLAB求各截取序列的DFT，对结果进行对比分析，讨论其产生差异的原因。

5、构建一个包含多个谐波分量的连续周期信号，以满足采样定理的不同采样频率对其进行离散化，得到离散序列（注意要选择适当的采样频率，保证离散化之后的序列仍然具有周期性），然后分别截取不同数量的序列周期（注意确保所截取长度恰好为整数周期），基于MATLAB，用DFT对信号进行频谱分析（用stem画出离散幅度频谱图）。根据不同采样频率（至少两种不同的采样频率）、不同截取周期数（至少两种不同的截取周期数）条件下的频谱分析结果，以及模拟频率、数字频率之间的对应关系，指出各谐波分量在DFT分析结果中的位置，讨论频率分辨率、频谱分析范围与采样频率、截取长度之间的关系。再用一个对最高谐波频率分量不满足采样定理（对其他频率分量满足）的采样频率，对连续信号取样并进行DFT频谱分析，讨论其频谱混叠现象。

6、分别针对取样信号Sa(t)、余弦信号cos(t)，用适当的采样频率进行离散化，截取一定的长度的序列（余弦信号可取整数或者非整数周期），基于Matlab平台，用DFT分析这两个模拟信号的幅度频谱并利用频域内插公式（参见第7次课第8个讨论题），画出连续的频谱曲线（注意应该基于DFT结果绘制出对应模拟信号的幅度频谱，频率轴标注为对应的模拟频率。请回顾第2章有关内容，结合教材3.4.2节，根据所用的采样频率、模拟信号频谱与序列频谱的对应关系、模拟频率与数字频率的对应关系以及DFT的频谱分析范围，对频谱分析结果进行适当处理之后进行绘图）。将所得结果与取样信号、余弦信号的理论频谱进行比较，据此讨论时域无限长信号频谱分析的截断效应，研究如何抑制这种截断效应。同时分析：①为什么对截取的序列补0之后，再求DFT，可以使画出的频谱曲线更光滑（画出补零之后的DFT进行对比）。②改变采样率和截取长度，将如何影响频谱分析结果？

1、用DFT计算线性卷积的原理步骤与优势

用DFT来计算序列的线性卷积有什么优势？其原理和步骤是什么（画出其实现原理框图，并予以说明）？根据实现原理框图，分析用DFT计算线性卷积的计算量，并与根据定义直接计算两个序列线性卷积的计算量相比，哪一个计算量更大？

DFT计算序列的线性卷积优势有：
1、线性卷积可以完全使用DFT实现
2、序列的DFT为离散序列，便于计算机处理，DFT可以使用其快速算法FFT大大降低计算量。
原理和步骤：

其中，在x1(n)和x2(n)进行DFT前需要补零：
设x(k)的长度为N，h(k)的长度为M，以补零的方式延长到L=N+M-1点，将两个延长后的序列的DFT相乘得到线性卷积。
用定义直接计算线性卷积与利用DFT计算线性卷积相比，前者的计算量更大。

2、用DFT计算线性卷积与直接计算线性卷积比较

代码：

%%      代码：
figure(1)
xn1=[1 1 1 1];                  %任取2个序列
xn2=[1 2 3 4];
subplot(2,2,1);                 %画图
stem(xn1);
title('x1(n)');
subplot(2,2,2);
stem(xn2);
title('x2(n)');Xn1=[1 1 1 1 0 0 0 0];          %取N>=N1+N2-1
Xn2=[1 2 3 4 0 0 0 0];          %利用DFT计算线性卷积
XK1=fft(Xn1)                    %快速傅里叶变换
XK2=fft(Xn2)
XK=XK1.*XK2;
xn=ifft (XK);                   %傅里叶反变换,至此便求得了线性卷积
subplot(2,2,3)                  %画图
stem(xn);
title('N1=8时用DFT计算线性卷积')Xn11=[1 1 1 1 0];               %取N1+N2-1>N>=max(N1,N2)
Xn22=[4 3 2 1 0];               %利用DFT计算线性卷积
XK11=fft (Xn11);                %快速傅里叶变换
XK22=fft (Xn22);
XK2=XK11.*XK22;
XN=ifft (XK2);                  %傅里叶反变换,至此便求得了线性卷积
subplot (2,2,4)                 %画图
stem(XN)
title('N2=5时用DFT计算线性卷积')figure(2)
y1=conv(Xn1,Xn2);
y2=conv(Xn11,Xn22);
subplot(2,1,1);
stem(y1);
title('N1=8时线性卷积')
subplot(2,1,2);
stem(y2);
title('N2=5时线性卷积')

运行结果：

分析：
x1(n)的长度为m，x2(n)的长度为k，线性卷积长度为N
则用DFT可以正确计算线性卷积的条件为：N>=m+k-1
如图，经过对比，N1满足条件，故可以还原，N2不满足条件，故不能还原，规律为将6、7点的数分别加到1、2点上。

3、用matlab用重叠相加法计算线性卷积

针对教材例3.4.1，不使用fftfilt函数，直接基于重叠相加法的原理，自行编程计算其卷积。据此总结重叠相加法的原理和步骤，并从理论上分析为什么卷积结果稳态输出是一个正弦序列。

代码：

%%      代码：
syms yn
n=0:40;
xn=cos(pi*n/10)+cos(2*pi*n/5);              %产生一个序列,参照教材例3.4.1
hn=[ones(1,5) zeros(1,36)];                 %产生1行41列的向量，前5列为1，其余为0
r1n=[ones(1,10) zeros(1,31)];               %先分段计算卷积，
r2n=[zeros(1,10) ones(1,10) zeros(1,21)];
r3n=[zeros(1,20) ones(1,10) zeros(1,11)];
r4n=[zeros(1,30) ones(1,10) 0];
x1n=xn.*r1n;
x2n=xn.*r2n;
x3n=xn.*r3n;
x4n=xn.*r4n;
y1n=conv(x1n,hn);
y2n=conv(x2n,hn);
y3n=conv(x3n,hn);
y4n=conv(x4n,hn);                           %至此，分段计算卷积结束
yn=y1n+y2n+y3n+y4n;                         %再把分段结果相加figure(1)                           %画图
subplot(4,1,1)                      %分段前x(n)
stem(x1n)
title('x1(n)-x4(n)')
subplot(4,1,2)
stem(x2n)
subplot(4,1,3)
stem(x3n)
subplot(4,1,4)
stem(x4n)figure(2)                           %分段前y(n)
subplot(4,1,1)
stem(y1n)
axis([0 45 -5 5])
title('y1(n)-y4(n)')
subplot(4,1,2)
stem(y2n)
axis([0 45 -5 5])
subplot(4,1,3)
stem(y3n)
axis([0 45 -5 5])
subplot(4,1,4)
stem(y4n)
axis([0 45 -5 5])figure(3)                           %合体后x(n)
stem(xn)
title('x(n)')figure(4)                           %合体后y(n)
stem(yn)
title('y(n)')
axis([0 45 -5 5])

运行结果：

重叠相加法就是把无限长序列x(n)等长分段，然后计算分段线性卷积，然后把分段卷积的结果叠加起来即可。因为输入序列是个余弦序列，所以输出也呈余弦变化。

课本溯源：
课本115页例3.4.1

代码：

%%      代码：
clear all
L=41;
N=5;
M=10;
hn=ones(1,N);
hn1=[hn zeros(1,L-N)];              %产生h(n)，补零是为了绘图好看
n=0:L-1;
xn=cos(pi*n/10)+cos(2*pi*n/5);        %产生x(n)的L个样值
yn=fftfilt(hn,xn,M);                %调用fftfilt计算卷积%以下为绘图部分
subplot(3,1,1);
stem(n,hn1,'.');
ylabel('h(n)');
axis([0,45,0,1.2])subplot(3,1,2);
stem(n,xn,'.');
ylabel('x(n)')
xlabel('n');
line([0,L+N-1],[0,0])subplot(3,1,3);
m=0:length(yn)-1;
stem(m,yn,'.');
xlabel('n');
ylabel('y(n)');
line([0,L+N-1],[0,0])

运行结果：

4、用DFT分析周期序列频谱

针对一个周期序列，分别截取其1个周期、多个周期、非整数周期，用MATLAB求各截取序列的DFT，对结果进行对比分析，讨论其产生差异的原因。
代码：

%%      代码：
n=0:79;
w=20;
xn=cos(2*pi/w*n);
L1=w;
L2=1.5*w;
L3=3*w;
r1=[ones(1,L1) zeros(1,length(xn)-L1)];  %一个周期
r2=[ones(1,L2) zeros(1,length(xn)-L2)];  %非整数周期
r3=[ones(1,L3) zeros(1,length(xn)-L3)];  %三个周期
y1=xn.*r1;
y2=xn.*r2;
y3=xn.*r3;                              %取xn的不同个周期
ft1=fft(y1,L1);
ft2=fft(y2,L2);
ft3=fft(y3,L3);                         %调用fft函数求各截取序列的DFT
figure                                   %绘制截取部分图像
subplot(3,1,1)
stem(n,y1)
subplot(3,1,2)
stem(n,y2)
subplot(3,1,3)
stem(n,y3)figure                                 %绘制截取部分幅频
subplot(3,1,1)
stem(abs(ft1))
title('abs1')
subplot(3,1,2)
stem(abs(ft2))
title('abs1.5')
subplot(3,1,3)
stem(abs(ft3))
title('abs3')figure                                %绘制截取部分相频
subplot(3,1,1)
stem(angle(ft1))
title('angle1')
subplot(3,1,2)
stem(angle(ft2))
title('angle1.5')
subplot(3,1,3)
stem(angle(ft3))
title('angle3')

运行结果：

当截取的长度为序列的周期以及倍数时，此时DFT为周期序列DFS的主值序列，取样点也为原来的倍数。

5、用DFT分析连续信号的频谱

构建一个包含多个谐波分量的连续周期信号，以满足采样定理的不同采样频率对其进行离散化，得到离散序列（注意要选择适当的采样频率，保证离散化之后的序列仍然具有周期性），然后分别截取不同数量的序列周期（注意确保所截取长度恰好为整数周期），基于MATLAB，用DFT对信号进行频谱分析（用stem画出离散幅度频谱图）。根据不同采样频率（至少两种不同的采样频率）、不同截取周期数（至少两种不同的截取周期数）条件下的频谱分析结果，以及模拟频率、数字频率之间的对应关系，指出各谐波分量在DFT分析结果中的位置，讨论频率分辨率、频谱分析范围与采样频率、截取长度之间的关系。再用一个对最高谐波频率分量不满足采样定理（对其他频率分量满足）的采样频率，对连续信号取样并进行DFT频谱分析，讨论其频谱混叠现象。

代码：

%%      代码：
clear;
f=600;
N=120;
T=1/f;
n=0:N-1;
xn=cos(2*pi*100*T*n)+2*cos(2*pi*200*T*n)+2*cos(2*pi*120*T*n);
figure(1)
stem(n,xn);hold on
t=0:0.0001:N/f;
xt=cos(2*pi*100*t)+2*cos(2*pi*200*t)+2*cos(2*pi*120*t);
plot(t*f,xt);
hold off
%N=120;
X=fft(xn,N);
XX=fftshift(X);
figure(2)
subplot(2,1,1)
stem((-N/2:N/2-1)/N/T,abs(XX)/N);
xlabel("f(Hz)");
title("离散化的截取模拟信号幅度频谱")w=-pi:0.001:pi;
Xw=0;
%下面基于频域采样内插公式，求出更多点的傅里叶变换，以画出连续的频谱曲线
for k=1:NXw=Xw+1/N*X(k)*sin((w-2*pi*(k-1)/N)/2*N)./sin((w-2*pi*(k-1)/N)/2).*exp(-j*(w-2*pi*(k-1)/N)*(N-1)/2);
end;
subplot(2,1,2)
plot(w/2/pi/T,abs(Xw)/N);
xlabel("f(Hz)");
title("截取模拟信号的幅度频谱(由离散频谱恢复得到的连续频谱)")

运行结果：

6、用matlab分析截断效应（以取样信号和余弦信号为例）

分别针对取样信号Sa(t)、余弦信号cos(t)，用适当的采样频率进行离散化，截取一定的长度的序列（余弦信号可取整数或者非整数周期），基于Matlab平台，用DFT分析这两个模拟信号的幅度频谱并利用频域内插公式（参见第7次课第8个讨论题），画出连续的频谱曲线（注意应该基于DFT结果绘制出对应模拟信号的幅度频谱，频率轴标注为对应的模拟频率。请回顾第2章有关内容，结合教材3.4.2节，根据所用的采样频率、模拟信号频谱与序列频谱的对应关系、模拟频率与数字频率的对应关系以及DFT的频谱分析范围，对频谱分析结果进行适当处理之后进行绘图）。将所得结果与取样信号、余弦信号的理论频谱进行比较，据此讨论时域无限长信号频谱分析的截断效应，研究如何抑制这种截断效应。同时分析：①为什么对截取的序列补0之后，再求DFT，可以使画出的频谱曲线更光滑（画出补零之后的DFT进行对比）。②改变采样率和截取长度，将如何影响频谱分析结果？

%%      代码：
clear;
T=1;
TL=40;                      %截断
n=-TL:TL;
N=length(n);
xn=cos(T*n);
yn=sinc(T/pi*n);
figure(1)                   %画Sa(t)函数和cos函数图
subplot(2,1,1)
stem(xn)
subplot(2,1,2)
stem(yn)figure(2)                   %画出Sa(t)函数和cos函数截断后的DFT
X=fft(xn,N);
subplot(2,1,1)
stem((0:N/2-1)/N*2*pi,abs(X(1:floor(N/2))))
Y=fft(yn,length(yn));
subplot(2,1,2)
stem((0:N/2-1)/N*2*pi,abs(Y(1:floor(N/2))))w=0:0.001:pi;
Xw=0;
Yw=0;
%下面基于频域采样内插公式，求出更多点的傅里叶变换，以画出连续的频谱曲线进行对比
for k=1:NXw=Xw+1/N*X(k)*sin((w-2*pi*(k-1)/N)/2*N)./sin((w-2*pi*(k-1)/N)/2).*exp(-j*(w-2*pi*(k-1)/N)*(N-1)/2);Yw=Yw+1/N*Y(k)*sin((w-2*pi*(k-1)/N)/2*N)./sin((w-2*pi*(k-1)/N)/2).*exp(-j*(w-2*pi*(k-1)/N)*(N-1)/2);
end;
figure(3)                           %将离散谱和连续谱放在同一张图中对比
subplot(2,1,1)
plot(w,abs(Xw));
hold on
stem((0:N/2-1)/N*2*pi,abs(X(1:floor(N/2))));
hold off
subplot(2,1,2)
plot(w,abs(Yw));
hold on
stem((0:N/2-1)/N*2*pi,abs(Y(1:floor(N/2))))
hold offN=N+200;                            %原来N=50左右，现在N又加了200，补了200个0，再进行DFT
figure(4)
X1=fft(xn,N);
subplot(2,1,1)
stem((0:N/2-1)/N*2*pi,abs(X1(1:floor(N/2))))
Y1=fft(yn,N);
subplot(2,1,2)
stem((0:N/2-1)/N*2*pi,abs(Y1(1:floor(N/2))))

运行结果:

仅改变代码第三行的TL：从40->100

TL=100;                      %截断

运行结果：

可见，频率分辨率明显变优。

分析：
1、由于截断，会导致频谱泄露，比如cos函数的频谱，不来应该是两根“线“”段，却成了一座“山”，Sa函数的频谱本来应该是一个平的“倒台阶”，结果截断后不平了，导致了频谱的不准；
2、DFT截断后补零，就是让描述频谱的点数增加了而已，仅有的效果就是“让一段区间内不准的频谱用更多点描述”，所以看起来频谱变光滑了；
3、虽然看起来变光滑了，但是由于频谱本身就是“不准的”（这个“不准”由截断引起，不由DFT的点数引起），所以虽然增加了DFT的点数，该不准还不准；
4、一旦增加了截断的长度（仅改变代码第三行的TL：从40->100），频谱瞬间变准了，“山变瘦了”，“台阶变平了”。

结论：
1、截断后补零的方法治标不治本，不能提高频率分辨率
2、要提高频率分辨率，就必须对原始信号截取的长度加长

三、反思总结

未理清第5题的逻辑