数字信号处理2： 离散信号与系统的频谱分析

---

前言

  本人初次学习离散信号与系统的频谱分析。

一、实验目的

  1.掌握离散傅里叶变换（DFT）及快速傅里叶变换（FFT）的计算机实现方法。
  2.检验序列DFT的性质。
  3.掌握利用DFT（FFT）计算序列线性卷积的方法。
  4.学习用DFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差，以便在实际中正确应用DFT。
  5.了解采样频率对谱分析的影响。
  6.了解利用FFT进行语音信号分析的方法。

二、实验设备

  1.计算机
  2.Matlab软件2014a以上版本。

三、实验内容

  1.对不同序列进行离散傅里叶变换并进行分析；DFT共轭对称性质的应用（通过1次N点FFT计算2个N点实序列的DFT）。
  2.线性卷积及循环卷积的关系，以及利用DFT（FFT）进行线性卷积的方法。
  3.比较计算序列的DFT和FFT的运算时间。
  4.利用FFT实现带噪信号检测。
  5.利用FFT计算信号频谱。
  6.扩展部分主要是关于离散系统采样频率、时域持续时间、谱分辨率等参数之间的关系，频谱的内插恢复，对语音信号进行简单分析。

四、实验原理

  1.序列的离散傅里叶变换及性质
  2.利用DFT计算线性卷积
  3.利用DFT对信号进行谱分析
  4.DFT和FFT的运算量
  5.信噪比

五、实验步骤

1.序列的离散傅里叶变换及分析

  分别对复序列，实序列，实偶序列，实奇序列，虚奇序列进行离散傅里叶变换，得到实验结果并对其特点进行分析。实验所需序列自选。
解：A、设定复序列：x(n)= {1,1,1,1,1}+j*{1,2,3,4,5} ；实序列：x1(n)= {1,1,1,1,1}
程序：

x1=ones(1,5);
x2=1:5;
x=x1+j*x2;
y1=fft(x,10);
y2=fft(x1,10);
subplot(1,2,1)
stem(0:9,abs(y1),'filled');
title('复序列DFT的幅度谱');
subplot(1,2,2)
stem(0:9,abs(y2),'filled');
title('实序列DFT的幅度谱');

结果：

B、设定实偶序列：x(n)= {5,4,3,2,1,0,1,2,3,4}，实奇序列：x(n)= {0,4,3,2,1,0,-1,-2,-3,-4}，虚奇序列：x(n)= j*{0,4,3,2,1,0,-1,-2,-3,-4};
程序：

x=[5,4,3,2,1,0,1,2,3,4];
y=fft(x,10);
subplot(3,3,1)
stem(0:9,abs(y),'filled');
title('实偶序列DFT的幅度谱');
subplot(3,3,2)
stem(0:9,real(y),'filled');
title('实偶序列DFT的实部');
subplot(3,3,3)
stem(0:9,imag(y),'filled');
title('实偶序列DFT的虚部');
x=[0,4,3,2,1,0,-1,-2,-3,-4];
y=fft(x,10);
subplot(3,3,4)
stem(0:9,abs(y),'filled');
title('实奇序列DFT的幅度谱');
subplot(3,3,5)
stem(0:9,real(y),'filled');
title('实奇序列DFT的实部');
subplot(3,3,6)
stem(0:9,abs(imag(y)),'filled');
title('实奇序列DFT的虚部');
x=j*[0,4,3,2,1,0,-1,-2,-3,-4];
y=fft(x,10);
subplot(3,3,7)
stem(0:9,y,'filled');
title('虚奇序列DFT的幅度谱');
subplot(3,3,8)
stem(0:9,real(y),'filled');
title('虚奇序列DFT的实部');
subplot(3,3,9)
stem(0:9,imag(y),'filled');0
title('虚奇序列DFT的虚部');

结果：

分析：

2.利用共轭对称性，设计高效算法计算2个N点实序列的DFT。

   用一个N点FFT计算两个长度为N的实序列N点离散傅里叶变换，并将结果和直接使用两个N点DFT得到的结果进行比较。
解：设定$x_1(n)$={1,2,3,4,5,6}，0≤n≤5; $x_2(n)$={1,2,1,3,1,4}，0≤n≤5;
程序：

 
x1=[1,2,3,4,5,6];
x2=[1,2,1,3,1,4];
y1=fft(x1);
y2=fft(x2);
x=x1+j*x2;
y=fft(x,6);
[xec,xoc]=circevod(y);
subplot(2,2,1)
stem(0:5,abs(y1),'filled');
title('x_1(n)的直接六点DFT的幅度谱');
grid on;
subplot(2,2,2)
stem(0:5,abs(y2),'filled');
title('x_2(n)的直接六点DFT的幅度谱');
grid on;
subplot(2,2,3)
stem(0:5,abs(xec),'filled');
title('x_1(n)高效算法DFT的幅度谱');
grid on;
subplot(2,2,4)
stem(0:5,abs(xoc),'filled');
title('x_2(n)高效算法DFT的幅度谱');
grid on;

结果：

分析：

3.线性卷积及循环卷积的实现及二者关系分析

  计算两序列的线性卷积及循环卷积，循环卷积采用2种计算方法（时域、频域方法）。设序列x1长度为M，序列x2长度为N，循环卷积长度为L，分别计算L大于、等于、小于（M+N-1）时的循环卷积。序列x1、x2、L自选，得到实验结果并对线性卷积及循环卷积的关系进行分析。
解：设定$x_1(n)$={1,2,3,4}，0≤n≤5; $x_2(n)$={1,2,1,3}，0≤n≤5;其中M=N=4，分别采用2种方法（时域、频域）计算L=4/7/12的循环卷积
程序：

figure(1)
x1=[1,2,3,4];
x2=[1,2,1,3];
subplot(3,2,1)
y1=circonv(x1,x2,4);
stem(0:3,abs(y1),'filled');
title('从时域计算4点循环卷积');
grid on;
subplot(3,2,3)
y1=circonv(x1,x2,7);
stem(0:6,abs(y1),'filled');
title('从时域计算7点循环卷积');
grid on;
subplot(3,2,5)
y1=circonv(x1,x2,12);
stem(0:11,abs(y1),'filled');
title('从时域计算12点循环卷积');
grid on;
subplot(3,2,2)
y1=fft(x1,4);
y2=fft(x2,4);
x=ifft(y1.*y2)
stem(0:3,abs(x),'filled');
title('从频域计算4点循环卷积');
grid on;
subplot(3,2,4)
y1=fft(x1,7);
y2=fft(x2,7);
x=ifft(y1.*y2)
stem(0:6,abs(x),'filled');
title('从频域计算7点循环卷积');
grid on;
subplot(3,2,6)
y1=fft(x1,12);
y2=fft(x2,12);
x=ifft(y1.*y2)
stem(0:11,abs(x),'filled');
title('从频域计算12点循环卷积');
grid on;figure(2)
x1=[1,2,3,4];
x2=[1,2,1,3];
y=conv(x1,x2);
stem(0:6,y,'filled')
grid on;
title('线性卷积的结果');

结果：

分析：

4.比较DFT和FFT的运算时间

  （1）自行选择进行计算的序列（或产生随机序列）。
  （2）DFT点数N分别取64、256、1024、2048、……
  （3）分别采用直接计算（按DFT定义计算）和fft函数计算，利用计时函数tic、toc，分别统计直接计算和fft计算所用的时间，并对结果进行比较（含直接计算和fft计算的比较、不同点数N所需计算时间和节约时间的比较等）。
解：设定序列为x=[1,2,3,4,ones(1,60)]
程序：

x=[1,2,3,4,ones(1,60)];
tic
X1=dft(x,64);
toc
tic
X2=fft(x,64);
toc
时间已过 0.003188 秒。
时间已过 0.001967 秒。
x=[1,2,3,4,ones(1,252)];
tic
X1=dft(x,256);
toc
tic
X2=fft(x,256);
toc
时间已过 0.012424 秒。
时间已过 0.002118 秒。
x=[1,2,3,4,ones(1,1020)];
tic
X1=dft(x,1024);
toc
tic
X2=fft(x,1024);
toc
时间已过 0.249468 秒。
时间已过 0.008443 秒。
x=[1,2,3,4,ones(1,2044)];
tic
X1=dft(x,2048);
toc
tic
X2=fft(x,2048);
toc
时间已过 1.078018 秒。
时间已过 0.021545 秒 
结果：
点数	DFT	FFT
64	0.003188 秒	0.001967 秒
256	0.012424 秒	0.002118 秒
1024	0.249468 秒	0.008443 秒
2048	1.078018 秒	0.021545 秒

结果：

点数	DFT	FFT
64	0.003188 秒	0.001967 秒
256	0.012424 秒	0.002118 秒
1024	0.249468 秒	0.008443 秒
2048	1.078018 秒	0.021545 秒

分析：

5.利用FFT求信号频谱及分析采样频率、噪声对频谱的影响

  （1）设模拟信号$x(t)=3cos(8\pi t)+6cos(20\pi t)$，以$t=0.01n,0\le n\le N-1$进行采样离散化，画出所得序列$x_N(n)$的N点DFT幅值谱（N分别取50，256），利用FFT实现。分析信号频谱所对应频率轴的数字频率和频率之间的关系。在点数N不变的条件下降低采样频率（需满足采样定理），观察信号频谱的变化，分析产生变化的原因。
解：A、采样频率为100hz，分别画出N=50，256点DFT幅度谱如下所示：
程序：

N=50;
n=0:N-1;
T=0.01;
x=3*cos(8*pi*n*T)+6*cos(20*pi*n*T);
subplot(2,2,1);
stem(n,x,'.');
title('f_s=100hz,N=50时的信号');
Y=fft(x,N);
k1=0:N-1;w1=2*pi/N*k1;
subplot(2,2,2);
stem(n,abs(Y),'.');
title('f_s=100hz,N=50时信号的频谱');
N=256;
n=0:N-1;
T=0.01;
x=3*cos(8*pi*n*T)+6*cos(20*pi*n*T);
subplot(2,2,3);
stem(n,x,'.');
title('f_s=100hz,N=256时的信号');
Y=fft(x,N);
k1=0:N-1;w1=2*pi/N*k1;
subplot(2,2,4);
stem(n,abs(Y),'.');
title('f_s=100hz,N=256时信号的频谱') 

结果：

B、采样频率为50hz，分别画出N=50，256点DFT幅度谱如下所示：
程序:

N=50;
n=0:N-1;
T=0.02;
x=3*cos(8*pi*n*T)+6*cos(20*pi*n*T);
subplot(2,2,1);
stem(n,x,'.');
title('f_s=50hz,N=50时的信号');
Y=fft(x,N);
k1=0:N-1;w1=2*pi/N*k1;
subplot(2,2,2);
stem(n,abs(Y),'.');
title('f_s=50hz,N=50时信号的频谱');
N=256
n=0:N-1;
T=0.02;
x=3*cos(8*pi*n*T)+6*cos(20*pi*n*T);
subplot(2,2,3);
stem(n,x,'.');
title('f_s=50hz,N=256时的信号');
Y=fft(x,N);
k1=0:N-1;w1=2*pi/N*k1;
subplot(2,2,4);
stem(n,abs(Y),'.');
title('f_s=50hz,N=256时信号的频谱');

结果：

分析：

  （2）对序列$x_N(n)$（N任取50或者256，需表示为行向量）加入噪声$（0.2\ast randn(1,N)）$，比较有无噪声时的信号谱，加大噪声到$2\ast randn(1,N)$和$10\ast randn(1,N)$，计算信噪比，画出并比较不同噪声下时域波形和幅度频谱，讨论噪声对信号分析的影响。
答：程序：

figure(1)
N=50;
n=0:N-1;
t=0.01*n;
q=n*2/N;
x1=3*cos(8*pi*t)+6*cos(20*pi*t);
subplot(2,2,1)
plot(n,x1);title('Pure signal')
y=fft(x1,N);
subplot(2,2,2)
plot(q,abs(y));title('FFT of pure signal(N=50)')
x2=x1+0.2*randn(1,N);
SNR1=snr(x1,x2)
subplot(2,2,3)
plot(n,x2);title('noisy signal')
y=fft(x2,N);
subplot(2,2,4)
plot(q,abs(y));title('FFT of signal with noise (N=50)')figure(2)
N=50;
n=0:N-1;
t=0.01*n;
q=n*2/N;
x1=3*cos(8*pi*t)+6*cos(20*pi*t);
x2=x1+2*randn(1,N);
x3=x1+10*randn(1,N);
SNR2=snr(x1,x2)
SNR3=snr(x1,x3)
subplot(2,2,1)
plot(n,x2);title('noisy signal（2*randn(1,N)）')
y2=fft(x2,N);
subplot(2,2,2)
plot(q,abs(y2));title('FFT of signal with more noise (N=50)')
subplot(2,2,3)
plot(n,x3);title('noisy signal（10*randn(1,N)）')
y3=fft(x3,N);
subplot(2,2,4)
plot(q,abs(y3));title('FFT of signal with more 
noise(N=50)') SNR1 =0.0783
SNR2 =-0.6612
SNR3 =-6.4308

图形结果：

分析：

6.创新训练拓展内容

6.1 信号持续时间、频谱分析范围、采样点数和谱分辨率的关系。

  （1）已知模拟信号$x(t)=cos(200\pi t)+cos(800\pi t)+cos(816\pi t)+cos(1000\pi t)$，欲采用数字方法分析该信号的频谱，要求能明确区分该信号的所有频率分量（两个相邻频谱分量之间至少间隔一个离散频率点）。为了便于采用基2-FFT来分析该信号的频谱并避免频谱失真，希望截取该信号的整数个周期进行频谱分析，并且采样点数为2M（M为正整数），请确定适当的采样频率（不发生频谱混叠且满足前述要求的最低采样频率）、采样点数（时域最小记录时间），并利用所确定的参数编写程序，对该信号进行频谱分析，画出模拟信号的幅度频谱，分析其频谱分析范围和谱分辨率。
解：$f_{smin}=1024hz，N_{min}=256$（具体过程在下面分析里）
程序

T=1/1024;
n=0:255;
t=n*T;
x=cos(200*pi*t)+cos(800*pi*t)+cos(816*pi*t)+cos(1000*pi*t);
subplot(2,1,1);
stem(n,x,'.') 
title('采样1024hz的信号');
axis([0,256,0,4])
subplot(2,1,2);
XK2=fft(x);
w=[0:5000]*pi/5000;
N=256;k=n;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi*512,abs(X2))
xlabel('hz');  
title('采样1024hz的幅度谱（N=256）');
axis([0,512,0,150]) 

结果：

分析：

  （2）保持采样点数不变，采样频率增大1倍，重新分析信号的幅度频谱。比较信号幅度频谱、谱分辨率、频谱分析范围的变化并分析其原因。
解：程序：

T=1/2048;
n=0:255;
t=n*T;
x=cos(200*pi*t)+cos(800*pi*t)+cos(816*pi*t)+cos(1000*pi*t);
subplot(2,1,1);
stem(n,x,'.') 
title('采样2048hz的信号');
axis([0,256,0,4])
subplot(2,1,2);
XK2=fft(x);
w=[0:5000]*pi/5000;
N=256;k=n;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi*1024,abs(X2))
xlabel('hz');  
title('采样2048hz的幅度谱（N=256）');
axis([0,1024,0,180]) 

结果：

分析：

  （3）保持采样频率不变，分别将采样点数增大和减少1倍，重新分析信号的幅度频谱。比较信号幅度频谱、谱分辨率、频谱分析范围的变化并分析其原因。
解：程序：

T=1/1024;
n=0:127;
t=n*T;
x=cos(200*pi*t)+cos(800*pi*t)+cos(816*pi*t)+cos(1000*pi*t);
subplot(2,2,1);
stem(n,x,'.') 
title('采样1024hz的信号');
axis([0,128,0,4])
subplot(2,2,2);
XK2=fft(x);
w=[0:5000]*pi/5000;
N=128;k=n;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi*512,abs(X2))
xlabel('hz');  
title('采样1024hz的幅度谱（N=128）');
axis([0,512,0,85]) 
T=1/1024;
n=0:511;
t=n*T;
x=cos(200*pi*t)+cos(800*pi*t)+cos(816*pi*t)+cos(1000*pi*t);
subplot(2,2,3);
stem(n,x,'.') 
title('采样1024hz的信号');
axis([0,512,0,4])
subplot(2,2,4);
XK2=fft(x);
w=[0:5000]*pi/5000;
N=512;k=n;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi*512,abs(X2))
xlabel('hz');  
title('采样1024hz的幅度谱（N=512）');
axis([0,512,0,300]) 

结果：

分析：

6.2频谱的内插函数恢复（频域采样恢复）

  基于频域采样X(k)及内插公式（教材式(3.3.7)和式(3.3.8)），用Matlab编程实现离散谱（DFT）到连续谱的转换：采用N=4的矩形序列$x(n)=R_4(n)$，分别求出其4点和8点DFT，然后利用内插公式（为保证恢复得到的频谱曲线足够光滑，在Matlab程序中对连续变量ω的离散化，应取200个离散点以上），分别由其4点和8点DFT恢复其连续频谱$X(e^{j\omega })$，画出$X(e^{j\omega })$的幅度频谱图（应为光滑连续曲线），并比较不同DFT点数恢复的结果以及$R_4(n)$的理论频谱（参见教材图2.2.1）。

解：程序：

subplot(2,2,1)
x=[1,1,1,1];
XK1=fft(x,4);
k=0:3;N=4;
w=[-2000:2000]*pi/1000;
X1=XK1/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi,abs(X1))
title('N=4的DFT内插公式恢复结果');
xlabel('w/π');
subplot(2,2,2)
n=0:3; 
x=(1).^n;  
w=[-500:500]*pi/250; 
X=x*exp(-j*n'*w);   
plot(w/pi,abs(X))
title('N=4的DFT的理论频谱');
xlabel('w/π');
subplot(2,2,3)
x=[1,1,1,1];
XK2=fft(x,8);
k=0:7;N=8;
w=[-2000:2000]*pi/1000;
X2=XK2/N*(sin(w*N/2-pi*k')./sin(w/2-pi*k'/N)).*exp(-j*w*(N-1)/2);
plot(w/pi,abs(X2))
title('N=8的DFT内插公式恢复结果');
xlabel('w/π');
subplot(2,2,4)
n=0:3; 
x=(1).^n 
w=[-500:500]*pi/250; 
X=x*exp(-j*(0:7)'*w);   
plot(w/pi,abs(X))
title('N=8的DFT的理论频谱');
xlabel('w/π');

结果：

分析：

6.3对语音信号进行简单分析。

  基于MATLAB平台，利用DFT对以适当的采样频率进行采集的不同人的单音节语音信号进行频谱分析（注意首先要对采集到的语音信号进行适当的截取），画出语音信号的波形及频谱。说明该频谱所对应的模拟频率，分析语音信号的频率分布特点，对比不同人说同一个音节的频谱的差异（最好分别用男声和女声进行分析）。对语音信号波形、单音节语音信号的截取方式、DFT频谱、频谱分析结果等进行讨论。
解：程序：

figure(1)
recObj = audiorecorder;
disp('Start speaking.')
recordblocking(recObj, 2);
disp('End of Recording.');
play(recObj);
y = getaudiodata(recObj);
subplot(1,2,1);
stem(y,'.');  
xlabel('时间');
ylabel('幅度');
title('信号波形');  
subplot(1,2,2); 
stem(abs(fft(y)),'.');
title('信号频谱');
xlabel('频率');
ylabel('幅度');figure(2)
sec1=5 
sec2 =12 
y1=y(((1000*sec1+1):1000*sec2),:);
subplot(1,2,1);
stem(y1,'.');   
xlabel('时间');
ylabel('幅度');
title('适当截取后信号的波形');  
subplot(1,2,2); 
stem(abs(fft(y1)),'.');
title('适当截取后信号的频谱');
xlabel('频率');
ylabel('幅度'); 

结果：

分析：

六、实验结论与心得体会（手写）

七、实验参考资料

  1.高西全,丁玉美.数字信号处理[M].西安：西安电子科技大学出版社,2008
  2.张德丰.详解MATLAB 数字信号处理[M].北京：电子工业出版社,2010
  3.王月明,张宝华.MATLAB基础与应用教程[M].北京:北京大学出版社,2012
  4.陈怀琛.数字信号处理教程：MATLAB释义与实现（第3版）[M].北京:电子工业出版社,2013
  5.宋知用.MATLAB数字信号处理85个实用案例精讲：入门到进阶[M].北京:北京航空航天大学出版社,2016